Skip to content

GitLab

Commit b4b6bb7e authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Function to report residual error for each unknown in model; changes to MFD... 

Function to report residual error for each unknown in model; changes to MFD scheme; tweaks to SuperLU download
parent f67c2fde
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ b4b6bb7e
......@@ -371,20 +371,21 @@ if(USE_SOLVER_SUPERLU)
CMAKE_ARGS "-DTPL_PARMETIS_INCLUDE_DIRS=${PARMETIS_INCLUDE_DIR} "
"-DTPL_PARMETIS_LIBRARIES=${PASS_PARMETIS_LIBRARIES} "
"-DTPL_LAPACK_LIBRARIES=${PASS_LAPACK_LIBRARIES} "
"-DBLAS_LIBRARIES=${PASS_BLAS_LIBRARIES} " #should come from find_pacakge(LAPACK)
"-DBLAS_LIBRARIES=${PASS_BLAS_LIBRARIES} "
"-Denable_examples=OFF "
"-Denable_tests=OFF "
"-Denable_doc=OFF "
"-DCMAKE_INSTALL_PREFIX=${LIB_DOWNLOAD_PATH} "
"-DUSE_XSDK_DEFAULTS=TRUE "
"-DXSDK_ENABLE_Fortran=FALSE "
"-DNOFORTRAN=TRUE "
"-DCMAKE_C_FLAGS='${C_STANDARD_FLAG} -DMETIS_EXPORT=' "
"-DUSE_XSDK_DEFAULTS=OFF "
"-DUSE_XSDK_DEFAULTS_DEFAULT=OFF "
"-DXSDK_ENABLE_Fortran=OFF "
"-DNOFORTRAN=ON "
"-DCMAKE_C_FLAGS='-DMETIS_EXPORT=' "
"-DCMAKE_CXX_FLAGS='-DMETIS_EXPORT=' "
"-DCMAKE_BUILD_TYPE=Release "
"-DMPI_C_COMPILER=${MPI_C_COMPILER} "
"-DMPI_CXX_COMPILER=${MPI_CXX_COMPILER} "
"-DCMAKE_DISABLE_FIND_PACKAGE_OpenMP=TRUE "
"-DCMAKE_DISABLE_FIND_PACKAGE_OpenMP=ON "
"-DCMAKE_POSITION_INDEPENDENT_CODE=ON "
)
......
Examples/ADMFD/diffusion.cpp
View file @ b4b6bb7e
......@@ -91,7 +91,7 @@ int main(int argc,char ** argv)
table[NORMAL] = FACE; //Compute normals
table[ORIENTATION] = FACE; //Check and fix normal orientation
table[MEASURE] = CELL | FACE; //Compute volumes and areas
//table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass
table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass
m->PrepareGeometricData(table); //Ask to precompute the data
BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
......@@ -177,9 +177,9 @@ int main(int argc,char ** argv)
#pragma omp parallel
#endif
{
rMatrix NK, R, Areas;
rMatrix NK, L, R, Areas;
rMatrix x(1,3), xf(1,3), n(1,3);
double area; //area of the face
double area, dist; //area of the face
double volume; //volume of the cell
#if defined(USE_OMP)
#pragma omp for
......@@ -192,29 +192,38 @@ int main(int argc,char ** argv)
int NF = (int)faces.size(); //number of faces;
rMatrix W(NF,NF);
volume = cell->Volume(); //volume of the cell
cell->Centroid(x.data());
cell->Barycenter(x.data());
//get permeability for the cell
rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),
cell->RealArrayDF(tag_K).size());
NK.Resize(NF,3); //co-normals
R.Resize(NF,3); //directions
L.Resize(NF,NF);
Areas.Resize(NF,NF); //areas
Areas.Zero();
L.Zero();
for(int k = 0; k < NF; ++k) //loop over faces
{
area = faces[k].Area();
faces[k].Centroid(xf.data());
faces[k].Barycenter(xf.data());
faces[k].OrientedUnitNormal(cell->self(),n.data());
dist = n.DotProduct(xf-x);
// assemble matrix of directions
R(k,k+1,0,3) = (xf-x)*area;
R(k,k+1,0,3) = (xf-x);
// assemble matrix of co-normals
NK(k,k+1,0,3) = n*K;
NK(k,k+1,0,3) = n*K*area;
L(k,k) = n.DotProduct(n*K)*area/dist;
Areas(k,k) = area;
} //end of loop over faces
W = NK*(NK.Transpose()*R).PseudoInvert(1.0e-12)*NK.Transpose(); //stability part
W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())*
(1.0/(static_cast<real>(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace());
W = Areas*W*Areas;
//~ W = NK*(NK.Transpose()*R).PseudoInvert(1.0e-12)*NK.Transpose(); //stability part
//~ W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())*
//~ (2.0/(static_cast<real>(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace());
double ca = 0, ba = ca;
W = (NK+ca*L*R)*((NK+ca*L*R).Transpose()*R).Invert()*(NK+ca*L*R).Transpose();
W+= L - (1+ba)*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
//~ W = Areas*W*Areas;
//access data structure for gradient matrix in mesh
real_array store_W = cell->RealArrayDV(tag_W);
//resize the structure
......@@ -223,7 +232,7 @@ int main(int argc,char ** argv)
std::copy(W.data(),W.data()+NF*NF,store_W.data());
tag_WG[cell].resize(3*NF);
tag_WG(cell,3,NF) = (NK.Transpose()*R).PseudoInvert(1.0e-12)*NK.Transpose()*Areas;
tag_WG(cell,3,NF) = (NK.Transpose()*R+ca*L*R).PseudoInvert(1.0e-12)*(NK+ca*L*R).Transpose();
} //end of loop over cells
}
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
......@@ -286,6 +295,7 @@ int main(int argc,char ** argv)
{
vMatrix pF; //vector of pressure differences on faces
vMatrix FLUX; //computed flux on faces
rMatrix n(3,1);
#if defined(USE_OMP)
#pragma omp for
#endif
......@@ -317,7 +327,15 @@ int main(int argc,char ** argv)
if( faces[k].Boundary() )
{
double a = 1, b = 0, c = 0;
double a = 0, b = 1, c = 0;
//~ faces[k].UnitNormal(n.data());
//~ double a = 1, b = 0, c = 0;
//~ if( fabs(n(2,0)-1) < 1.0e-4 )
//~ {
//~ a = 0;
//~ b = 1;
//~ }
if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
{
real_array BC = faces[k].RealArray(tag_BC);
......@@ -354,8 +372,8 @@ int main(int argc,char ** argv)
if( R.Norm() < 1.0e-4 ) break;
//Solver S(Solver::INNER_ILU2);
Solver S(Solver::INNER_MPTILUC);
//~ Solver S(Solver::K3BIILU2);
//~ Solver S(Solver::INNER_MPTILUC);
Solver S(Solver::K3BIILU2);
//Solver S("superlu");
S.SetParameter("verbosity","1");
S.SetParameter("relative_tolerance", "1.0e-14");
......
Examples/TestSuite/ns_turek.cpp
View file @ b4b6bb7e
......@@ -120,28 +120,31 @@ int main(int argc, char ** argv)
if( it->nbAdjElements(FACE,cylinder) ) it->SetMarker(cylinder);
std::cout << "project boundary nodes onto cylinder " << std::endl;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) if( it->GetMarker(cylinder) )
if( fix_cylinder != 3 )
{
double x = it->Coords()[0], y = it->Coords()[1], dx, dy;
double r = sqrt((x-0.5)*(x-0.5) + (y-0.2)*(y-0.2));
if( r )
std::cout << "project boundary nodes onto cylinder " << std::endl;
for(Mesh::iteratorNode it = m->BeginNode(); it != m->EndNode(); ++it) if( it->GetMarker(cylinder) )
{
dx = (0.05/r-1.0)*(x-0.5);
dy = (0.05/r-1.0)*(y-0.2);
//std::cout << "at " << x << "," << y << " r " << r << " dx " << dx << " dy " << dy << std::endl;
it->Coords()[0] += dx;
it->Coords()[1] += dy;
double x = it->Coords()[0], y = it->Coords()[1], dx, dy;
double r = sqrt((x-0.5)*(x-0.5) + (y-0.2)*(y-0.2));
if( r )
{
dx = (0.05/r-1.0)*(x-0.5);
dy = (0.05/r-1.0)*(y-0.2);
//std::cout << "at " << x << "," << y << " r " << r << " dx " << dx << " dy " << dy << std::endl;
it->Coords()[0] += dx;
it->Coords()[1] += dy;
}
else std::cout << "node at center of cylinder: " << x << "," << y << std::endl;
}
else std::cout << "node at center of cylinder: " << x << "," << y << std::endl;
}
else fix_cylinder = 2;
if( fix_cylinder == 2 )
{
TagRealArray vec_t = m->CreateTag("vec_t",DATA_REAL,FACE,FACE,2);
std::cout << "project centers of boundary faces onto cylinder " << std::endl;
int iter = 0;
while(iter < 1000)
while(iter < 5000)
{
double A = 0, err = 0, fA;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->GetMarker(cylinder) )
......
Source/Autodiff/model.cpp
View file @ b4b6bb7e
......@@ -197,9 +197,17 @@ namespace INMOST
bool success = true;
R.Clear();
Automatizator::MakeCurrent(&aut);
double total_time = Timer(), model_time = 0;
for(std::vector< std::pair<std::string, AbstractSubModel *> >::const_iterator it = SubModels.begin();
it != SubModels.end(); ++it)
{
model_time = Timer();
success &= it->second->FillResidual(R);
model_time = Timer() - model_time;
//~ std::cout << it->first << ":" << model_time << " ";
}
total_time = Timer() - total_time;
//~ std::cout << "total:" << total_time << std::endl;
Automatizator::RemoveCurrent();
return success;
}
......@@ -287,5 +295,63 @@ namespace INMOST
it != SubModels.end(); ++it)
it->second->Adaptation(m,tree);
}
void Model::ReportErrors(const Residual & R) const
{
for(std::vector< std::pair< std::string, AbstractEntry *> >::const_iterator it = Entries.begin(); it != Entries.end(); ++it)
{
ElementType etypes = it->second->GetElementType();
Mesh * m = it->second->GetMeshLink();
//if( m->GetProcessorRank() == 0 )
// std::cout << "Errors for entry " << it->first << std::endl;
for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype) ) if ( etype & etypes )
{
if( !m->HaveGlobalID(etype) ) m->AssignGlobalID(etype); //to report element number
//account for variable size!
std::vector<double> err_comp(it->second->Size(),0.0);
//std::vector<Element> err_elem(it->second->Size());
std::vector<double> err_int(it->second->Size(),0.0);
double max_err = 0, int_err = 0;
Element e;
for(Mesh::iteratorElement jt = m->BeginElement(etype); jt != m->EndElement(); ++jt) if( jt->GetStatus() != Element::Ghost )
{
rpMatrix err = R.Value(it->second->Index(jt->self()));
double block_err = err.FrobeniusNorm();
if( block_err > max_err )
{
e = jt->self();
max_err = block_err;
}
int_err += block_err*block_err;
int N = err.Rows()*err.Cols();
if( N > it->second->Size() ) continue; //No account for variable size
//~ std::cout << jt->LocalID() << " block err " << block_err << " comp err ";
for(int k = 0; k < N; ++k)
{
//~ std::cout << err.data()[k] << " ";
err_int[k] += err.data()[k]*err.data()[k];
if( fabs(err.data()[k]) > err_comp[k] )
{
err_comp[k] = fabs(err.data()[k]);
//err_elem[k] = jt->self();
}
}
//~ std::cout << std::endl;
}
m->Integrate(&err_int[0],it->second->Size());
m->AggregateMax(&err_comp[0],it->second->Size());
int_err = m->Integrate(int_err);
max_err = m->AggregateMax(max_err); //we need MPI_Allreduce to zero proc to know index of element with largest error
if( m->GetProcessorRank() == 0 )
{
std::cout << std::setw(30) << it->first << " element type " << ElementTypeName(etype) << " error integral " << std::setw(12) << sqrt(int_err) << " maximal " << std::setw(12) << max_err << std::endl;// " on element " << e.GlobalID() << std::endl;
for(int k = 0; k < (int)err_int.size(); ++k)
{
std::cout << "\t\t" << std::setw(2) << k << " error integral " << std::setw(12) << sqrt(err_int[k]) << " maximal " << std::setw(12) << err_comp[k] << std::endl;
}
}
}
}
}
}
#endif
Source/Autodiff/residual.cpp
View file @ b4b6bb7e
......@@ -109,6 +109,15 @@ namespace INMOST
new (&ret(i,j)) multivar_expression_reference(residual[rows(i,j)],&jacobian[rows(i,j)]);
return ret;
}
rpMatrix Residual::Value(const AbstractMatrix<INMOST_DATA_INTEGER_TYPE> & rows) const
{
rpMatrix ret(rows.Rows(),rows.Cols());
for(INMOST_DATA_ENUM_TYPE i = 0; i < rows.Rows(); ++i)
for(INMOST_DATA_ENUM_TYPE j = 0; j < rows.Cols(); ++j)
ret(i,j) = residual[rows(i,j)];
return ret;
}
#endif //USE_SOLVER
} //namespace INMOST
......
Source/Headers/inmost_model.h
View file @ b4b6bb7e
......@@ -147,6 +147,8 @@ namespace INMOST
/// Adapt the data of the model after the mesh refinement/coarsement.
/// Those model that use the adapted mesh should update their data
virtual void Adaptation(Mesh & m) const;
void ReportErrors(const Residual & R) const;
};
}
......
Source/Headers/inmost_residual.h
View file @ b4b6bb7e
......@@ -50,10 +50,18 @@ namespace INMOST
/// @return A structure that can be used in or assigned an automatic differentiation expression.
__INLINE multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row)
{return multivar_expression_reference(residual[row],&jacobian[row]);}
/// Retrive a residual value corresponding to certain equation.
/// @param row Equation number.
/// @return A structure that can be used in or assigned an automatic differentiation expression.
__INLINE double Value(INMOST_DATA_ENUM_TYPE row) const {return residual[row];}
/// Retrive a vector of entries in residual, corresponding to a set of equations.
/// @param rows A row-vector of equation numbers.
/// @param A structure that can be used in or assigned an automatic differentiation matrix expression.
Matrix<multivar_expression_reference> operator [](const AbstractMatrix<INMOST_DATA_INTEGER_TYPE> & rows);
/// Retrive a vector of entries in residual, corresponding to a set of equations.
/// @param rows A row-vector of equation numbers.
/// @param A structure that can be used in or assigned an automatic differentiation matrix expression.
rpMatrix Value(const AbstractMatrix<INMOST_DATA_INTEGER_TYPE> & rows) const;
/// Retrive hessian matrix. Use in nonlinear solver.
Sparse::HessianMatrix & GetHessian() {return hessian;}
/// Retrive hessian matrix without right of modificaiton.
......
Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp
View file @ b4b6bb7e
......@@ -19,7 +19,7 @@ using namespace INMOST;
#define REORDER_RCM
//#define REORDER_NNZ
#if defined(USE_SOLVER_METIS)
//~ #define REORDER_METIS_ND
#define REORDER_METIS_ND
#endif
#if defined(USE_SOLVER_MONDRIAAN)
//#define REORDER_MONDRIAAN
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment