v0.15.0
Loading...
Searching...
No Matches
radiation.cpp
Go to the documentation of this file.
1/**
2 * \file lesson6_radiation.cpp
3 * \example lesson6_radiation.cpp
4 *
5 * Using PipelineManager interface calculate the divergence of base functions,
6 * and integral of flux on the boundary. Since the h-div space is used, volume
7 * integral and boundary integral should give the same result.
8 */
9
10
11
12#include <MoFEM.hpp>
13
14using namespace MoFEM;
15
16static char help[] = "...\n\n";
17
18#include <MoFEM.hpp>
19
21using DomainEleOp = DomainEle::UserDataOperator;
30
31// Units
32// Temperature: Kelvins
33// Length: 1 km
34// Time: 1 s
35
36constexpr double heat_conductivity = ((0.4 + 0.7) / 2) * 1e3;
37
38constexpr double emissivity = 1;
39constexpr double boltzmann_constant = 5.670374419e-2;
40constexpr double Beta = emissivity * boltzmann_constant;
41
42constexpr double T_ambient = 2.7;
43struct Example {
44
45 Example(MoFEM::Interface &m_field) : mField(m_field) {}
46
48
49private:
51
52 static int integrationRule(int, int, int p_data) { return 2 * p_data; };
53
61
62 boost::shared_ptr<VectorDouble> approxVals;
63 boost::shared_ptr<MatrixDouble> approxGradVals;
64
65 struct OpRadiationLhs : public OpBase {
66
67 private:
68 boost::shared_ptr<VectorDouble> approxVals;
69
70 public:
71 OpRadiationLhs(boost::shared_ptr<VectorDouble> &approx_vals)
72 : OpBase("T", "T", OpBase::OPROWCOL), approxVals(approx_vals) {
73 this->sYmm = false;
74 }
75
76 MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
77 };
78
79 struct OpRadiationRhs : public OpBase {
80
81 private:
82 boost::shared_ptr<VectorDouble> approxVals;
83
84 public:
85 OpRadiationRhs(boost::shared_ptr<VectorDouble> &approx_vals)
86 : OpBase("T", "T", OpBase::OPROW), approxVals(approx_vals) {}
87
89 };
90
91 struct OpFluxRhs : public OpBase {
92
93 private:
94 FTensor::Index<'i', 2> i; ///< summit Index
95
96 public:
97 OpFluxRhs() : OpBase("T", "T", OpBase::OPROW) {}
98
100 };
101
103
104 private:
105 boost::shared_ptr<VectorDouble> approxVals;
107 double &surfaceArea;
108
109 public:
111 boost::shared_ptr<VectorDouble> &approx_vals, double &sum_temp,
112 double &surf)
113 : EdgeEleOp("T", "T", OpBase::OPROW), approxVals(approx_vals),
114 sumTemperature(sum_temp), surfaceArea(surf) {}
115
116 MoFEMErrorCode doWork(int side, EntityType type,
118 };
119};
120
125 CHKERR bC();
126 CHKERR OPs();
131}
132
133//! [Set up problem]
137 // Add field
139 CHKERR simple->addBoundaryField("T", H1, AINSWORTH_LEGENDRE_BASE, 1);
140 constexpr int order = 3;
141 CHKERR simple->setFieldOrder("T", order);
142 CHKERR simple->setUp();
144}
145//! [Set up problem]
146
147//! [Create common data]
150 approxVals = boost::make_shared<VectorDouble>();
151 approxGradVals = boost::make_shared<MatrixDouble>();
153}
154//! [Create common data]
155
156//! [Boundary condition]
159 // Set initial values
160 auto set_initial_temperature = [&](VectorAdaptor &&field_data, double *xcoord,
161 double *ycoord, double *zcoord) {
163 field_data[0] = T_ambient;
165 };
166 FieldBlas *field_blas;
167 CHKERR mField.getInterface(field_blas);
168 CHKERR field_blas->setVertexDofs(set_initial_temperature, "T");
170}
171//! [Boundary condition]
172
173//! [Push operators to pipeline]
177 auto beta = [](const double r, const double, const double) {
178 return heat_conductivity * (2 * M_PI * r);
179 };
180
181 auto det_ptr = boost::make_shared<VectorDouble>();
182 auto jac_ptr = boost::make_shared<MatrixDouble>();
183 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
184
185 pipeline_mng->getOpDomainLhsPipeline().push_back(
186 new OpCalculateHOJac<2>(jac_ptr));
187 pipeline_mng->getOpDomainLhsPipeline().push_back(
188 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
189 pipeline_mng->getOpDomainLhsPipeline().push_back(
190 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
191 pipeline_mng->getOpDomainLhsPipeline().push_back(
193
194 pipeline_mng->getOpDomainLhsPipeline().push_back(
195 new OpDomainGradGrad("T", "T", beta));
197
198 pipeline_mng->getOpDomainRhsPipeline().push_back(
199 new OpCalculateHOJac<2>(jac_ptr));
200 pipeline_mng->getOpDomainRhsPipeline().push_back(
201 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
202 pipeline_mng->getOpDomainRhsPipeline().push_back(
203 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
204 pipeline_mng->getOpDomainRhsPipeline().push_back(
206
207 pipeline_mng->getOpDomainRhsPipeline().push_back(
209 pipeline_mng->getOpDomainRhsPipeline().push_back(
210 new OpDomainGradTimesVec("T", approxGradVals, beta));
212
213 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
215 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
216 new OpRadiationRhs(approxVals));
217 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFluxRhs());
219
220 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
222 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
223 new OpRadiationLhs(approxVals));
226}
227//! [Push operators to pipeline]
228
229//! [Solve]
234 auto ts = pipeline_mng->createTSIM();
235
236 double ftime = 1;
237 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
238 CHKERR TSSetFromOptions(ts);
239 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
240
241 auto T = createDMVector(simple->getDM());
242 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
243 SCATTER_FORWARD);
244
245 CHKERR TSSolve(ts, T);
246 CHKERR TSGetTime(ts, &ftime);
247
248 PetscInt steps, snesfails, rejects, nonlinits, linits;
249 CHKERR TSGetTimeStepNumber(ts, &steps);
250 CHKERR TSGetSNESFailures(ts, &snesfails);
251 CHKERR TSGetStepRejections(ts, &rejects);
252 CHKERR TSGetSNESIterations(ts, &nonlinits);
253 CHKERR TSGetKSPIterations(ts, &linits);
254 MOFEM_LOG_C("EXAMPLE", Sev::inform,
255 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
256 "%d, linits %d",
257 steps, rejects, snesfails, ftime, nonlinits, linits);
258
260}
261//! [Solve]
262
263//! [Postprocess results]
267 pipeline_mng->getDomainLhsFE().reset();
268 pipeline_mng->getBoundaryLhsFE().reset();
269 pipeline_mng->getBoundaryRhsFE().reset();
270 auto post_proc_fe =
271 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
272
273 auto t_ptr = boost::make_shared<VectorDouble>();
274 post_proc_fe->getOpPtrVector().push_back(
275 new OpCalculateScalarFieldValues("T", t_ptr));
276
278
279 post_proc_fe->getOpPtrVector().push_back(
280
281 new OpPPMap(post_proc_fe->getPostProcMesh(),
282 post_proc_fe->getMapGaussPts(),
283
284 {{"T", t_ptr}},
285
286 {}, {}, {}
287
288 )
289
290 );
291
292 pipeline_mng->getDomainRhsFE() = post_proc_fe;
293
294 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
295 new OpCalculateScalarFieldValues("T", approxVals));
296
297 double sum_temperature;
298 double surface_area;
299 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
300 new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
301 surface_area));
302 auto calc_surfcae_area_op = pipeline_mng->getOpBoundaryRhsPipeline().back();
303
304 sum_temperature = 0;
305 surface_area = 0;
306 CHKERR pipeline_mng->loopFiniteElements();
307 CHKERR post_proc_fe->writeFile("out_radiation.h5m");
308
309 MOFEM_LOG_C("EXAMPLE", Sev::inform, "Surface area %3.4e [km]", surface_area);
310 MOFEM_LOG_C("EXAMPLE", Sev::inform,
311 "Average subsurface temperatute %3.4e [K]",
312 sum_temperature / surface_area);
313
315}
316//! [Postprocess results]
317
318//! [Check results]
320//! [Check results]
321
322int main(int argc, char *argv[]) {
323
324 // Initialisation of MoFEM/PETSc and MOAB data structures
325 const char param_file[] = "param_file.petsc";
326 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
327
328 // Add logging channel for example
329 auto core_log = logging::core::get();
330 core_log->add_sink(
332 LogManager::setLog("EXAMPLE");
333 MOFEM_LOG_TAG("EXAMPLE", "example");
334
335 try {
336
337 //! [Register MoFEM discrete manager in PETSc]
338 DMType dm_name = "DMMOFEM";
339 CHKERR DMRegister_MoFEM(dm_name);
340 //! [Register MoFEM discrete manager in PETSc
341
342 //! [Create MoAB]
343 moab::Core mb_instance; ///< mesh database
344 moab::Interface &moab = mb_instance; ///< mesh database interface
345 //! [Create MoAB]
346
347 //! [Create MoFEM]
348 MoFEM::Core core(moab); ///< finite element database
349 MoFEM::Interface &m_field = core; ///< finite element database insterface
350 //! [Create MoFEM]
351
352 //! [Load mesh]
353 Simple *simple = m_field.getInterface<Simple>();
355 CHKERR simple->loadFile("");
356 //! [Load mesh]
357
358 //! [Example]
359 Example ex(m_field);
360 CHKERR ex.runProblem();
361 //! [Example]
362 }
364
366}
367
368//! [Radiation Lhs]
370 EntData &col_data) {
372 // get element volume
373 const double vol = this->getMeasure();
374 // get integration weights
375 auto t_w = getFTensor0IntegrationWeight();
376 // get base function gradient on rows
377 auto t_row_base = row_data.getFTensor0N();
378 // gat temperature at integration points
379 auto t_val = getFTensor0FromVec(*(approxVals));
380 // get coordinate at integration points
381 auto t_coords = getFTensor1CoordsAtGaussPts();
382
383 // loop over integration points
384 for (int gg = 0; gg != nbIntegrationPts; gg++) {
385
386 // Cylinder radius
387 const double r_cylinder = t_coords(0);
388
389 // take into account Jacobean
390 const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
391 // loop over rows base functions
392 for (int rr = 0; rr != nbRows; ++rr) {
393 auto t_col_base = col_data.getFTensor0N(gg, 0);
394 // loop over columns
395 for (int cc = 0; cc != nbCols; cc++) {
396 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
397 locMat(rr, cc) += alpha * t_row_base * t_col_base * 4 *
398 std::pow(static_cast<double>(t_val), 3);
399 }
400 ++t_col_base;
401 }
402 ++t_row_base;
403 }
404
405 ++t_val;
406 ++t_coords;
407 ++t_w; // move to another integration weight
408 }
410}
411//! [Radiation Lhs]
412
413//! [Radiation Lhs]
416 // get element volume
417 const double vol = getMeasure();
418 // get integration weights
419 auto t_w = getFTensor0IntegrationWeight();
420 // get base function gradient on rows
421 auto t_row_base = row_data.getFTensor0N();
422 // gat temperature at integration points
423 auto t_val = getFTensor0FromVec(*(approxVals));
424 // get coordinate at integration points
425 auto t_coords = getFTensor1CoordsAtGaussPts();
426
427 // loop over integration points
428 for (int gg = 0; gg != nbIntegrationPts; gg++) {
429
430 // Cylinder radius
431 const double r_cylinder = t_coords(0);
432
433 // take into account Jacobean
434 const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
435 // loop over rows base functions
436 for (int rr = 0; rr != nbRows; ++rr) {
437 if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
438 locF[rr] += alpha * t_row_base *
439 (std::pow(static_cast<double>(t_val), 4) -
440 std::pow(T_ambient, 4));
441 }
442 ++t_row_base;
443 }
444 ++t_coords;
445 ++t_val;
446 ++t_w; // move to another integration weight
447 }
448
450}
451//! [Radiation Lhs]
452
453//! [Flux Rhs]
456 // get element volume
457 const double vol = getMeasure();
458 // get integration weights
459 auto t_w = getFTensor0IntegrationWeight();
460 // get base function gradient on rows
461 auto t_row_base = row_data.getFTensor0N();
462 // get coordinate at integration points
463 auto t_coords = getFTensor1CoordsAtGaussPts();
464 // // get time
465 const double time = getFEMethod()->ts_t;
466
467 // Look to https://doi.org/10.1016/j.icarus.2014.12.028s
468 constexpr double flux_p = -0.03e6;
469 constexpr double flux_c = -0.23e6;
470
471 // loop over integration points
472 for (int gg = 0; gg != nbIntegrationPts; gg++) {
473
474 // Cylinder radius
475 const double r_cylinder = t_coords(0);
476
477 const double r = std::sqrt(t_coords(i) * t_coords(i));
478 const double s = std::abs(t_coords(1)) / r;
479
480 // take into account Jacobean
481 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
482 // loop over rows base functions
483 for (int rr = 0; rr != nbRows; ++rr) {
484 locF[rr] += alpha * t_row_base * (s * flux_p + flux_c) * time;
485 ++t_row_base;
486 }
487 ++t_coords;
488 ++t_w; // move to another integration weight
489 }
490
492}
493//! [Flux Rhs]
494
495//! [Ave Temp]
497 int side, EntityType type, EntitiesFieldData::EntData &data) {
498
500 if (type == MBVERTEX) {
501 // get element volume
502 const double vol = getMeasure();
503 // get integration weights
504 auto t_w = getFTensor0IntegrationWeight();
505 // gat temperature at integration points
506 auto t_val = getFTensor0FromVec(*(approxVals));
507 // get coordinate at integration points
508 auto t_coords = getFTensor1CoordsAtGaussPts();
509 // number of integration pts
510 size_t nb_integration_pts = getGaussPts().size2();
511
512 // loop over integration points
513 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
514
515 // Cylinder radius
516 const double r_cylinder = t_coords(0);
517
518 // take into account Jacobean
519 const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
520
521 sumTemperature += alpha * t_val;
522 surfaceArea += alpha;
523
524 ++t_coords;
525 ++t_val;
526 ++t_w; // move to another integration weight
527 }
528 }
530}
531
532//! [Ave Temp]
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition helmholtz.cpp:24
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
int r
Definition sdf.py:8
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static char help[]
Definition radiation.cpp:16
constexpr double Beta
Definition radiation.cpp:40
constexpr double boltzmann_constant
Definition radiation.cpp:39
constexpr double T_ambient
Definition radiation.cpp:42
constexpr double emissivity
Definition radiation.cpp:38
constexpr double heat_conductivity
Definition radiation.cpp:36
OpCalcSurfaceAverageTemperature(boost::shared_ptr< VectorDouble > &approx_vals, double &sum_temp, double &surf)
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Flux Rhs]
FTensor::Index< 'i', 2 > i
summit Index
Definition radiation.cpp:94
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
boost::shared_ptr< VectorDouble > approxVals
Definition radiation.cpp:68
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
OpRadiationLhs(boost::shared_ptr< VectorDouble > &approx_vals)
Definition radiation.cpp:71
OpRadiationRhs(boost::shared_ptr< VectorDouble > &approx_vals)
Definition radiation.cpp:85
boost::shared_ptr< VectorDouble > approxVals
Definition radiation.cpp:82
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
[Example]
Definition plastic.cpp:216
static int integrationRule(int, int, int p_data)
Definition radiation.cpp:52
boost::shared_ptr< VectorDouble > approxVals
Definition radiation.cpp:62
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition radiation.cpp:45
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
boost::shared_ptr< MatrixDouble > approxGradVals
Definition radiation.cpp:63
MoFEM::Interface & mField
Definition plastic.cpp:226
MoFEMErrorCode postProcess()
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
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:118
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
@ OPROW
operator doWork function is executed on FE rows
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.