v0.14.0
Loading...
Searching...
No Matches
operators_tests.cpp

Test operators in forms integrators.

Test operators in forms integrators

Date
2022-12-11

TODO: Add more operators.

/**
* @file operators_tests.cpp
* @example operators_tests.cpp
* @brief Test operators in forms integrators
* @date 2022-12-11
*
* @copyright Copyright (c) 2022
*
* TODO: Add more operators.
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
using DomainEle = VolumeElementForcesAndSourcesCore;
};
//! [Define dimension]
constexpr int SPACE_DIM = 3;
constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType I =
IntegrationType::GAUSS; //< selected integration type
// Specializations for tested operators
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
template <int FIELD_DIM>
constexpr bool debug = false;
int main(int argc, char *argv[]) {
// initialize petsc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
// Create MoAB database
moab::Core moab_core;
moab::Interface &moab = moab_core;
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab);
MoFEM::Interface &m_field = mofem_core;
// Register DM Manager
DMType dm_name = "DMMOFEM";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "OpTester"));
LogManager::setLog("OpTester");
MOFEM_LOG_TAG("OpTester", "OpTester");
// Simple interface
auto simple = m_field.getInterface<Simple>();
// get options from command line
CHKERR simple->getOptions();
// load mesh file
CHKERR simple->loadFile();
// Scalar fields and vector field is tested. Add more fields, i.e. vector
// field if needed.
CHKERR simple->addDomainField("SCALAR", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addDomainField("VECTOR", H1, AINSWORTH_LEGENDRE_BASE,
// set fields order, i.e. for most first cases order is sufficient.
CHKERR simple->setFieldOrder("SCALAR", 1);
CHKERR simple->setFieldOrder("VECTOR", 1);
// setup problem
CHKERR simple->setUp();
// get operators tester
auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
// OperatorsTester
auto pip = m_field.getInterface<PipelineManager>(); // get interface to
// pipeline manager
auto post_proc = [&](auto dm, auto f_res, auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
auto scal_vec = boost::make_shared<VectorDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("SCALAR", scal_vec, f_res));
post_proc_fe->getOpPtrVector().push_back(
f_res));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"SCALAR", scal_vec}},
{{"VECTOR", vec_mat}},
{}, {})
);
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
// Note: I do not push OPs transferring base to physical coordinates. That
// is not needed to test OPs, and not required, since we like to test only
// operators.
// Test grad-grad operator for scalar and vector field specialisation
auto TestOpGradGrad = [&]() {
MOFEM_LOG("OpTester", Sev::verbose) << "TestOpGradGrad";
pip->getOpDomainLhsPipeline().clear();
pip->getOpDomainRhsPipeline().clear();
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<1>(
"SCALAR", "SCALAR", [](double, double, double) { return 1; }));
pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<SPACE_DIM>(
"VECTOR", "VECTOR", [](double, double, double) { return 1; }));
auto scl_mat = boost::make_shared<MatrixDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
pip->getOpDomainRhsPipeline().push_back(
pip->getOpDomainRhsPipeline().push_back(
vec_mat));
pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<1>(
"SCALAR", scl_mat, [](double, double, double) { return 1; }));
pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<SPACE_DIM>(
"VECTOR", vec_mat, [](double, double, double) { return 1; }));
constexpr double eps = 1e-6;
auto x = opt->setRandomFields(simple->getDM(),
{{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto diff_x = opt->setRandomFields(
simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
pip->getDomainLhsFE(), x, SmartPetscObj<Vec>(), SmartPetscObj<Vec>(),
diff_x, 0, 1, eps);
// Calculate norm of difference between directive calculated from finite
// difference, and tangent matrix.
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C("OpTester", Sev::inform, "TestOpGradGrad %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Norm of directional derivative too large");
if (debug) {
// Example how to plot direction in direction diff_x. If instead
// directionalCentralFiniteDifference(...) diff_res is used, then error
// on directive is plotted.
CHKERR post_proc(simple->getDM(),
//
opt->directionalCentralFiniteDifference(
simple->getDM(), simple->getDomainFEName(),
pip->getDomainRhsFE(), x, SmartPetscObj<Vec>(),
SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
//
"out_directional_directive_gradgrad.h5m");
}
};
// test operator for convection (part of material time directives)
auto TestOpConvection = [&]() {
MOFEM_LOG("OpTester", Sev::verbose) << "TestOpConvection";
pip->getOpDomainLhsPipeline().clear();
pip->getOpDomainRhsPipeline().clear();
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
auto evaluate_field = [&](auto &pipe) {
auto scl_mat = boost::make_shared<MatrixDouble>();
auto vec_mat = boost::make_shared<MatrixDouble>();
auto dot_vec_vel = boost::make_shared<MatrixDouble>();
pipe.push_back(
"VECTOR", vec_mat));
"VECTOR", dot_vec_vel));
return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
};
auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
evaluate_field(pip->getOpDomainRhsPipeline());
pip->getOpDomainRhsPipeline().push_back(
new OpConvectiveTermRhs<1>("SCALAR", rhs_dot_vec_vel, rhs_scl_mat));
pip->getOpDomainRhsPipeline().push_back(
new OpConvectiveTermRhs<SPACE_DIM>("VECTOR", rhs_dot_vec_vel,
rhs_vec_mat));
auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
evaluate_field(pip->getOpDomainLhsPipeline());
// I create op, set scaling function to calculate time directive, and add
// operator pinter to pipeline
auto op_convective_term_lhs_du_scalar =
new OpConvectiveTermLhsDu<1>("SCALAR", "VECTOR", lhs_scl_mat);
op_convective_term_lhs_du_scalar->feScalingFun =
[](const FEMethod *fe_ptr) { return fe_ptr->ts_a; };
pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
pip->getOpDomainLhsPipeline().push_back(
new OpConvectiveTermLhsDy<1>("SCALAR", "SCALAR", lhs_dot_vec_vel));
auto op_convective_term_lhs_du_vec =
new OpConvectiveTermLhsDu<SPACE_DIM>("VECTOR", "VECTOR", lhs_vec_mat);
op_convective_term_lhs_du_vec->feScalingFun = [](const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
pip->getOpDomainLhsPipeline().push_back(
new OpConvectiveTermLhsDy<SPACE_DIM>("VECTOR", "VECTOR",
lhs_dot_vec_vel));
constexpr double eps = 1e-6;
auto x = opt->setRandomFields(simple->getDM(),
{{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto x_t = opt->setRandomFields(
simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto diff_x = opt->setRandomFields(
simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
auto diff_res = opt->checkCentralFiniteDifference(
simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
pip->getDomainLhsFE(), x, x_t, SmartPetscObj<Vec>(), diff_x, 0, 1,
eps);
// Calculate norm of difference between directive calculated from finite
// difference, and tangent matrix.
double fnorm;
CHKERR VecNorm(diff_res, NORM_2, &fnorm);
MOFEM_LOG_C("OpTester", Sev::inform, "TestOpConvection %3.4e", fnorm);
constexpr double err = 1e-9;
if (fnorm > err) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Norm of directional derivative too large");
}
if (debug) {
// Example how to plot direction in direction diff_x. If instead
// directionalCentralFiniteDifference(...) diff_res is used, then error
// on directive is plotted.
CHKERR post_proc(simple->getDM(),
//
opt->directionalCentralFiniteDifference(
simple->getDM(), simple->getDomainFEName(),
pip->getDomainRhsFE(), x, x_t,
SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
//
"out_directional_directive_convection.h5m");
}
};
CHKERR TestOpGradGrad();
CHKERR TestOpConvection();
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}
static Index< 'o', 3 > o
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
static char help[]
int main()
static const double eps
constexpr int SPACE_DIM
constexpr int FIELD_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
static const bool debug
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr IntegrationType I
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.