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

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 = 4
 

Function Documentation

◆ main()

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

Definition at line 309 of file reaction_diffusion.cpp.

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

Variable Documentation

◆ help

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

Definition at line 24 of file reaction_diffusion.cpp.