v0.9.1
Classes | Namespaces | Typedefs | Functions | Variables
reaction_diffusion_equation.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 = DataForcesAndSourcesCore::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_equation.cpp.

Definition at line 310 of file reaction_diffusion_equation.cpp.

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

Variable Documentation

◆ help

char help[] = "...\n\n"
static