v0.13.1
Classes | Namespaces | Typedefs | Functions | Variables
reaction_diffusion.cpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  ReactionDiffusionEquation::CommonData
 Common data. More...
 
struct  ReactionDiffusionEquation::OpAssembleMass
 Assemble mass matrix. More...
 
struct  ReactionDiffusionEquation::OpAssembleSlowRhs
 Assemble slow part. More...
 
struct  ReactionDiffusionEquation::OpAssembleStiffRhs< DIM >
 Assemble stiff part. More...
 
struct  ReactionDiffusionEquation::OpAssembleStiffLhs< DIM >
 Assemble stiff part tangent. More...
 
struct  ReactionDiffusionEquation::Monitor
 Monitor solution. More...
 

Namespaces

namespace  ReactionDiffusionEquation
 

Typedefs

using ReactionDiffusionEquation::Ele = FaceElementForcesAndSourcesCore
 
using ReactionDiffusionEquation::OpEle = FaceElementForcesAndSourcesCore::UserDataOperator
 
using ReactionDiffusionEquation::EntData = EntitiesFieldData::EntData
 
using ReactionDiffusionEquation::PostProcEle = PostProcBrokenMeshInMoab< Ele >
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
const double ReactionDiffusionEquation::D = 2e-3
 diffusivity More...
 
const double ReactionDiffusionEquation::r = 1
 rate factor More...
 
const double ReactionDiffusionEquation::k = 1
 caring capacity More...
 
const double ReactionDiffusionEquation::u0 = 0.1
 inital vale on blocksets More...
 
const int ReactionDiffusionEquation::save_every_nth_step = 1
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
reaction_diffusion.cpp.

Definition at line 296 of file reaction_diffusion.cpp.

296 {
297
298 // initialize petsc
299 const char param_file[] = "param_file.petsc";
300 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
301
302 try {
303
304 // Create moab and mofem instances
305 moab::Core mb_instance;
306 moab::Interface &moab = mb_instance;
307 MoFEM::Core core(moab);
308 MoFEM::Interface &m_field = core;
309
310 // Register DM Manager
311 DMType dm_name = "DMMOFEM";
312 CHKERR DMRegister_MoFEM(dm_name);
313
314 // Simple interface
315 Simple *simple_interface;
316 CHKERR m_field.getInterface(simple_interface);
317 CHKERR simple_interface->getOptions();
318 CHKERR simple_interface->loadFile();
319
320 int order = 4; ///< approximation order
321 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
322
323 // add fields
324 CHKERR simple_interface->addDomainField("u", H1, AINSWORTH_LEGENDRE_BASE,
325 1);
326 // set fields order
327 CHKERR simple_interface->setFieldOrder("u", order);
328 // setup problem
329 CHKERR simple_interface->setUp();
330
331 // Create common data structure
332 boost::shared_ptr<CommonData> data(new CommonData());
333 /// Alias pointers to data in common data structure
334 auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
335 auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
336 auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
337
338 // Create finite element instances to integrate the right-hand side of slow
339 // and stiff vector, and the tangent left-hand side for stiff part.
340 boost::shared_ptr<Ele> vol_ele_slow_rhs(new Ele(m_field));
341 boost::shared_ptr<Ele> vol_ele_stiff_rhs(new Ele(m_field));
342 boost::shared_ptr<Ele> vol_ele_stiff_lhs(new Ele(m_field));
343
344 // Push operators to integrate the slow right-hand side vector
345 vol_ele_slow_rhs->getOpPtrVector().push_back(
346 new OpCalculateScalarFieldValues("u", val_ptr));
347 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
348
349 // PETSc IMAX and Explicit solver demans that g = M^-1 G is provided. So
350 // when the slow right-hand side vector (G) is assembled is solved for g
351 // vector.
352 auto solve_for_g = [&]() {
353 MoFEMFunctionBegin;
354 if (vol_ele_slow_rhs->vecAssembleSwitch) {
355 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
356 SCATTER_REVERSE);
357 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
358 SCATTER_REVERSE);
359 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
360 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
361 *vol_ele_slow_rhs->vecAssembleSwitch = false;
362 }
363 CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
364 vol_ele_slow_rhs->ts_F);
365 MoFEMFunctionReturn(0);
366 };
367 // Add hook to the element to calculate g.
368 vol_ele_slow_rhs->postProcessHook = solve_for_g;
369
370 auto det_ptr = boost::make_shared<VectorDouble>();
371 auto jac_ptr = boost::make_shared<MatrixDouble>();
372 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
373 // Add operators to calculate the stiff right-hand side
374 vol_ele_stiff_rhs->getOpPtrVector().push_back(
375 new OpCalculateHOJac<2>(jac_ptr));
376 vol_ele_stiff_rhs->getOpPtrVector().push_back(
377 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
379 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
380 vol_ele_stiff_rhs->getOpPtrVector().push_back(
381 new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
382 vol_ele_stiff_rhs->getOpPtrVector().push_back(
383 new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
384 vol_ele_stiff_rhs->getOpPtrVector().push_back(
385 new OpAssembleStiffRhs<2>(data));
386
387 // Add operators to calculate the stiff left-hand side
388 vol_ele_stiff_lhs->getOpPtrVector().push_back(
389 new OpCalculateHOJac<2>(jac_ptr));
390 vol_ele_stiff_lhs->getOpPtrVector().push_back(
391 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
392 vol_ele_stiff_lhs->getOpPtrVector().push_back(
393 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
394 vol_ele_stiff_lhs->getOpPtrVector().push_back(
395 new OpAssembleStiffLhs<2>(data));
396
397 // Set integration rules
398 auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
399 vol_ele_slow_rhs->getRuleHook = vol_rule;
400 vol_ele_stiff_rhs->getRuleHook = vol_rule;
401 vol_ele_stiff_lhs->getRuleHook = vol_rule;
402
403 // Crate element for post-processing
404
405 auto post_proc = boost::make_shared<PostProcEle>(m_field);
406 boost::shared_ptr<ForcesAndSourcesCore> null;
407
408 using OpPPMap = OpPostProcMapInMoab<2, 2>;
409
410 auto u_ptr = boost::make_shared<VectorDouble>();
411 post_proc->getOpPtrVector().push_back(
412 new OpCalculateScalarFieldValues("u", u_ptr));
413 post_proc->getOpPtrVector().push_back(
414
415 new OpPPMap(
416
417 post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
418
419 {{"u", u_ptr}},
420
421 {},
422
423 {},
424
425 {}
426
427 )
428
429 );
430
431 // Get PETSc discrete manager
432 auto dm = simple_interface->getDM();
433
434 // Get surface entities form blockset, set initial values in those
435 // blocksets. To keep it simple is assumed that inital values are on
436 // blockset 1
437 if (m_field.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
438 Range inner_surface;
439 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
440 1, BLOCKSET, 2, inner_surface, true);
441 if (!inner_surface.empty()) {
442 Range inner_surface_verts;
443 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
444 CHKERR m_field.getInterface<FieldBlas>()->setField(
445 u0, MBVERTEX, inner_surface_verts, "u");
446 }
447 }
448
449 // Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
450 // conditions on that boundary.
451 Range surface;
452 CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
453 Skinner skin(&m_field.get_moab());
454 Range edges;
455 CHKERR skin.find_skin(0, surface, false, edges);
456 Range edges_part;
457 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
458 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
459 PSTATUS_NOT, -1, &edges_part);
460 Range edges_verts;
461 CHKERR moab.get_connectivity(edges_part, edges_verts, false);
462 // Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
463 // the problem.
464 CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
465 simple_interface->getProblemName(), "u",
466 unite(edges_verts, edges_part));
467
468 // Create mass matrix, calculate and assemble
469 CHKERR DMCreateMatrix_MoFEM(dm, data->M);
470 CHKERR MatZeroEntries(data->M);
471 boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
472 vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
473 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
474 vol_mass_ele);
475 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
476 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
477
478 // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
479 // M^-1G(t,u)
480 data->ksp = createKSP(m_field.get_comm());
481 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
482 CHKERR KSPSetFromOptions(data->ksp);
483 CHKERR KSPSetUp(data->ksp);
484
485 // Create and setup TS solver
486 auto ts = createTS(m_field.get_comm());
487 // Use IMEX solver, i.e. implicit/explicit solver
488 CHKERR TSSetType(ts, TSARKIMEX);
489 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
490
491 // Add element to calculate lhs of stiff part
492 CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
493 vol_ele_stiff_lhs, null, null);
494 // Add element to calculate rhs of stiff part
495 CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
496 vol_ele_stiff_rhs, null, null);
497 // Add element to calculate rhs of slow (nonlinear) part
498 CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
499 vol_ele_slow_rhs, null, null);
500
501 // Add monitor to time solver
502 boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
503 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
504 monitor_ptr, null, null);
505
506 // Create solution vector
507 SmartPetscObj<Vec> X;
508 CHKERR DMCreateGlobalVector_MoFEM(dm, X);
509 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
510
511 // Solve problem
512 double ftime = 1;
513 CHKERR TSSetDM(ts, dm);
514 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
515 CHKERR TSSetSolution(ts, X);
516 CHKERR TSSetFromOptions(ts);
517 CHKERR TSSolve(ts, X);
518 }
519 CATCH_ERRORS;
520
521 // finish work cleaning memory, getting statistics, etc.
522 MoFEM::Core::Finalize();
523
524 return 0;
525}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
param_file
std::string param_file
Definition: DefaultParams.hpp:113
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
H1
@ H1
continuous field
Definition: definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
MoFEM::DMMoFEMTSSetIFunction
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:747
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1144
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
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: DMMMoFEM.cpp:533
MoFEM::DMMoFEMTSSetIJacobian
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:800
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1114
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:272
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:247
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:1003
MoFEM::DMMoFEMTSSetRHSFunction
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS the right hand side function
Definition: DMMMoFEM.cpp:829
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:235
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1086
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:262
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
OpRotatingFrameTools::CommonData
OpPlasticTools::CommonData CommonData
Definition: RotatingFrameOperators.hpp:17
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Ele
FaceElementForcesAndSourcesCore Ele
Definition: quad_polynomial_approximation.cpp:13
help
static char help[]
Definition: reaction_diffusion.cpp:10
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
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::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1224
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:83
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:583
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEM::Simple::addDomainField
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:374
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93

Variable Documentation

◆ help

char help[] = "...\n\n"
static
Examples
reaction_diffusion.cpp.

Definition at line 10 of file reaction_diffusion.cpp.

Generated on Thu Sep 22 2022 19:26:31 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow