v0.14.0
operators_tests.cpp
Go to the documentation of this file.
1 /**
2  * @file operators_tests.cpp
3  * @example operators_tests.cpp
4  * @brief Test operators in forms integrators
5  * @date 2022-12-11
6  *
7  * @copyright Copyright (c) 2022
8  *
9  * TODO: Add more operators.
10  *
11  */
12 
13 #include <MoFEM.hpp>
14 using namespace MoFEM;
15 
16 static char help[] = "...\n\n";
17 
18 template <int DIM> struct ElementsAndOps {};
19 
20 template <> struct ElementsAndOps<2> {
22 };
23 
24 template <> struct ElementsAndOps<3> {
26 };
27 
28 //! [Define dimension]
29 constexpr int SPACE_DIM = 3;
30 constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
31 constexpr IntegrationType I =
32  IntegrationType::GAUSS; //< selected integration type
33 
38 
39 
40 // Specializations for tested operators
41 
42 template <int FIELD_DIM>
45 
46 template <int FIELD_DIM>
49 
50 template <int FIELD_DIM>
53 
54 template <int FIELD_DIM>
57 
58 template <int FIELD_DIM>
61 
62 constexpr bool debug = false;
63 
64 int main(int argc, char *argv[]) {
65 
66  // initialize petsc
67  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
68 
69  try {
70 
71  // Create MoAB database
72  moab::Core moab_core;
73  moab::Interface &moab = moab_core;
74 
75  // Create MoFEM database and link it to MoAB
76  MoFEM::Core mofem_core(moab);
77  MoFEM::Interface &m_field = mofem_core;
78 
79  // Register DM Manager
80  DMType dm_name = "DMMOFEM";
81  CHKERR DMRegister_MoFEM(dm_name);
82 
83  // Add logging channel for example
84  auto core_log = logging::core::get();
85  core_log->add_sink(
87  LogManager::setLog("OpTester");
88  MOFEM_LOG_TAG("OpTester", "OpTester");
89 
90  // Simple interface
91  auto simple = m_field.getInterface<Simple>();
92 
93  // get options from command line
94  CHKERR simple->getOptions();
95  // load mesh file
96  CHKERR simple->loadFile();
97 
98  // Scalar fields and vector field is tested. Add more fields, i.e. vector
99  // field if needed.
100  CHKERR simple->addDomainField("SCALAR", H1, AINSWORTH_LEGENDRE_BASE, 1);
101  CHKERR simple->addDomainField("VECTOR", H1, AINSWORTH_LEGENDRE_BASE,
102  SPACE_DIM);
103 
104  // set fields order, i.e. for most first cases order is sufficient.
105  CHKERR simple->setFieldOrder("SCALAR", 1);
106  CHKERR simple->setFieldOrder("VECTOR", 1);
107 
108  // setup problem
109  CHKERR simple->setUp();
110 
111  // get operators tester
112  auto opt = m_field.getInterface<OperatorsTester>(); // get interface to
113  // OperatorsTester
114  auto pip = m_field.getInterface<PipelineManager>(); // get interface to
115  // pipeline manager
116 
117  auto post_proc = [&](auto dm, auto f_res, auto out_name) {
119  auto post_proc_fe =
120  boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
121 
123 
124  auto scal_vec = boost::make_shared<VectorDouble>();
125  auto vec_mat = boost::make_shared<MatrixDouble>();
126 
127  post_proc_fe->getOpPtrVector().push_back(
128  new OpCalculateScalarFieldValues("SCALAR", scal_vec, f_res));
129  post_proc_fe->getOpPtrVector().push_back(
130  new OpCalculateVectorFieldValues<SPACE_DIM>("VECTOR", vec_mat,
131  f_res));
132 
133  post_proc_fe->getOpPtrVector().push_back(
134 
135  new OpPPMap(
136 
137  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
138 
139  {{"SCALAR", scal_vec}},
140 
141  {{"VECTOR", vec_mat}},
142 
143  {}, {})
144 
145  );
146 
147  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
148  post_proc_fe);
149  post_proc_fe->writeFile(out_name);
151  };
152 
153  // Note: I do not push OPs transferring base to physical coordinates. That
154  // is not needed to test OPs, and not required, since we like to test only
155  // operators.
156 
157  // Test grad-grad operator for scalar and vector field specialisation
158  auto TestOpGradGrad = [&]() {
160  MOFEM_LOG("OpTester", Sev::verbose) << "TestOpGradGrad";
161 
162  pip->getOpDomainLhsPipeline().clear();
163  pip->getOpDomainRhsPipeline().clear();
164 
165  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
166  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
167 
168  pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<1>(
169  "SCALAR", "SCALAR", [](double, double, double) { return 1; }));
170  pip->getOpDomainLhsPipeline().push_back(new OpGradGrad<SPACE_DIM>(
171  "VECTOR", "VECTOR", [](double, double, double) { return 1; }));
172 
173  auto scl_mat = boost::make_shared<MatrixDouble>();
174  auto vec_mat = boost::make_shared<MatrixDouble>();
175 
176  pip->getOpDomainRhsPipeline().push_back(
177  new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
178  pip->getOpDomainRhsPipeline().push_back(
180  vec_mat));
181 
182  pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<1>(
183  "SCALAR", scl_mat, [](double, double, double) { return 1; }));
184  pip->getOpDomainRhsPipeline().push_back(new OpGradTimesTensor<SPACE_DIM>(
185  "VECTOR", vec_mat, [](double, double, double) { return 1; }));
186 
187  constexpr double eps = 1e-6;
188 
189  auto x = opt->setRandomFields(simple->getDM(),
190  {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
191  auto diff_x = opt->setRandomFields(
192  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
193 
194  auto diff_res = opt->checkCentralFiniteDifference(
195  simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
196  pip->getDomainLhsFE(), x, SmartPetscObj<Vec>(), SmartPetscObj<Vec>(),
197  diff_x, 0, 1, eps);
198 
199  // Calculate norm of difference between directive calculated from finite
200  // difference, and tangent matrix.
201  double fnorm;
202  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
203  MOFEM_LOG_C("OpTester", Sev::inform, "TestOpGradGrad %3.4e", fnorm);
204 
205  constexpr double err = 1e-9;
206  if (fnorm > err)
207  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
208  "Norm of directional derivative too large");
209 
210  if (debug) {
211  // Example how to plot direction in direction diff_x. If instead
212  // directionalCentralFiniteDifference(...) diff_res is used, then error
213  // on directive is plotted.
214  CHKERR post_proc(simple->getDM(),
215  //
216  opt->directionalCentralFiniteDifference(
217  simple->getDM(), simple->getDomainFEName(),
218  pip->getDomainRhsFE(), x, SmartPetscObj<Vec>(),
219  SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
220  //
221  "out_directional_directive_gradgrad.h5m");
222  }
223 
225  };
226 
227  // test operator for convection (part of material time directives)
228  auto TestOpConvection = [&]() {
230  MOFEM_LOG("OpTester", Sev::verbose) << "TestOpConvection";
231 
232  pip->getOpDomainLhsPipeline().clear();
233  pip->getOpDomainRhsPipeline().clear();
234 
235  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
236  pip->setDomainLhsIntegrationRule([](int, int, int o) { return 2 * o; });
237 
238  auto evaluate_field = [&](auto &pipe) {
239  auto scl_mat = boost::make_shared<MatrixDouble>();
240  auto vec_mat = boost::make_shared<MatrixDouble>();
241  auto dot_vec_vel = boost::make_shared<MatrixDouble>();
242  pipe.push_back(
243  new OpCalculateScalarFieldGradient<SPACE_DIM>("SCALAR", scl_mat));
245  "VECTOR", vec_mat));
247  "VECTOR", dot_vec_vel));
248  return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
249  };
250 
251  auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
252  evaluate_field(pip->getOpDomainRhsPipeline());
253  pip->getOpDomainRhsPipeline().push_back(
254  new OpConvectiveTermRhs<1>("SCALAR", rhs_dot_vec_vel, rhs_scl_mat));
255  pip->getOpDomainRhsPipeline().push_back(
256  new OpConvectiveTermRhs<SPACE_DIM>("VECTOR", rhs_dot_vec_vel,
257  rhs_vec_mat));
258 
259  auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
260  evaluate_field(pip->getOpDomainLhsPipeline());
261 
262  // I create op, set scaling function to calculate time directive, and add
263  // operator pinter to pipeline
264  auto op_convective_term_lhs_du_scalar =
265  new OpConvectiveTermLhsDu<1>("SCALAR", "VECTOR", lhs_scl_mat);
266  op_convective_term_lhs_du_scalar->feScalingFun =
267  [](const FEMethod *fe_ptr) { return fe_ptr->ts_a; };
268  pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
269  pip->getOpDomainLhsPipeline().push_back(
270  new OpConvectiveTermLhsDy<1>("SCALAR", "SCALAR", lhs_dot_vec_vel));
271 
272  auto op_convective_term_lhs_du_vec =
273  new OpConvectiveTermLhsDu<SPACE_DIM>("VECTOR", "VECTOR", lhs_vec_mat);
274  op_convective_term_lhs_du_vec->feScalingFun = [](const FEMethod *fe_ptr) {
275  return fe_ptr->ts_a;
276  };
277  pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
278  pip->getOpDomainLhsPipeline().push_back(
279  new OpConvectiveTermLhsDy<SPACE_DIM>("VECTOR", "VECTOR",
280  lhs_dot_vec_vel));
281 
282  constexpr double eps = 1e-6;
283 
284  auto x = opt->setRandomFields(simple->getDM(),
285  {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
286  auto x_t = opt->setRandomFields(
287  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
288  auto diff_x = opt->setRandomFields(
289  simple->getDM(), {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
290 
291  auto diff_res = opt->checkCentralFiniteDifference(
292  simple->getDM(), simple->getDomainFEName(), pip->getDomainRhsFE(),
293  pip->getDomainLhsFE(), x, x_t, SmartPetscObj<Vec>(), diff_x, 0, 1,
294  eps);
295 
296  // Calculate norm of difference between directive calculated from finite
297  // difference, and tangent matrix.
298  double fnorm;
299  CHKERR VecNorm(diff_res, NORM_2, &fnorm);
300  MOFEM_LOG_C("OpTester", Sev::inform, "TestOpConvection %3.4e", fnorm);
301 
302  constexpr double err = 1e-9;
303  if (fnorm > err) {
304  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
305  "Norm of directional derivative too large");
306  }
307 
308  if (debug) {
309  // Example how to plot direction in direction diff_x. If instead
310  // directionalCentralFiniteDifference(...) diff_res is used, then error
311  // on directive is plotted.
312  CHKERR post_proc(simple->getDM(),
313  //
314  opt->directionalCentralFiniteDifference(
315  simple->getDM(), simple->getDomainFEName(),
316  pip->getDomainRhsFE(), x, x_t,
317  SmartPetscObj<Vec>(), diff_res, 0, 1, eps),
318  //
319  "out_directional_directive_convection.h5m");
320  }
321 
323  };
324 
325  CHKERR TestOpGradGrad();
326  CHKERR TestOpConvection();
327  }
328  CATCH_ERRORS;
329 
330  // finish work cleaning memory, getting statistics, etc.
332 
333  return 0;
334 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
Definition: operators_tests.cpp:44
OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
Definition: operators_tests.cpp:52
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
OpConvectiveTermLhsDu
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
Definition: operators_tests.cpp:56
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FIELD_DIM
constexpr int FIELD_DIM
Definition: child_and_parent.cpp:15
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
main
int main(int argc, char *argv[])
Definition: operators_tests.cpp:64
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:310
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
Definition: operators_tests.cpp:60
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
BiLinearForm
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition: OperatorsTester.hpp:21
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
EntData
EntitiesFieldData::EntData EntData
Definition: operators_tests.cpp:34
ElementsAndOps
Definition: child_and_parent.cpp:18
debug
constexpr bool debug
Definition: operators_tests.cpp:62
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
MoFEM::DMoFEMLoopFiniteElements
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:586
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: operators_tests.cpp:29
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
help
static char help[]
Definition: operators_tests.cpp:16