v0.14.0
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

 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[] 
)

< approximation order

Alias pointers to data in common data structure

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 = [&]() {
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);
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 
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
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.
523 
524  return 0;
525 }

Variable Documentation

◆ help

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

Definition at line 10 of file reaction_diffusion.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
surface
Definition: surface.py:1
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:245
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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: DMMoFEM.cpp:527
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:257
ReactionDiffusionEquation::OpAssembleSlowRhs
Assemble slow part.
Definition: reaction_diffusion.cpp:102
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1171
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1201
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
Ele
FaceElementForcesAndSourcesCore Ele
Definition: quad_polynomial_approximation.cpp:13
ReactionDiffusionEquation::OpAssembleStiffLhs
Assemble stiff part tangent.
Definition: reaction_diffusion.cpp:203
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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: DMMoFEM.cpp:857
ReactionDiffusionEquation::OpAssembleStiffRhs
Assemble stiff part.
Definition: reaction_diffusion.cpp:151
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
help
static char help[]
Definition: reaction_diffusion.cpp:10
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
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:264
surface.surface
def surface(x, y, z, eta)
Definition: surface.py:3
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
ReactionDiffusionEquation::Monitor
Monitor solution.
Definition: reaction_diffusion.cpp:269
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
ReactionDiffusionEquation::CommonData
Common data.
Definition: reaction_diffusion.cpp:34
Range
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: DMMoFEM.cpp:886
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
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
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: DMMoFEM.cpp:1060
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
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: DMMoFEM.cpp:804
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
ReactionDiffusionEquation::OpAssembleMass
Assemble mass matrix.
Definition: reaction_diffusion.cpp:47
MoFEM::SmartPetscObj< Vec >
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
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:357
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
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:590
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698