v0.14.0
Loading...
Searching...
No Matches
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>
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
18template <int DIM> struct ElementsAndOps {};
19
20template <> struct ElementsAndOps<2> {
22};
23
24template <> struct ElementsAndOps<3> {
26};
27
28//! [Define dimension]
29constexpr int SPACE_DIM = 3;
30constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
31constexpr IntegrationType I =
32 IntegrationType::GAUSS; //< selected integration type
33
38
39
40// Specializations for tested operators
41
42template <int FIELD_DIM>
45
46template <int FIELD_DIM>
49
50template <int FIELD_DIM>
53
54template <int FIELD_DIM>
57
58template <int FIELD_DIM>
61
62constexpr bool debug = false;
63
64int 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 }
329
330 // finish work cleaning memory, getting statistics, etc.
332
333 return 0;
334}
static Index< 'o', 3 > o
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
static const double eps
constexpr int FIELD_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ 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 ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
static char help[]
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
constexpr bool debug
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
constexpr IntegrationType I
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
constexpr AssemblyType A
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
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)
structure for User Loop Methods on finite elements
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
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector 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.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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.