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";
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
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 }
520
521 // finish work cleaning memory, getting statistics, etc.
523
524 return 0;
525}
static Index< 'p', 3 > p
std::string param_file
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
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
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
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1144
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
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
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
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1114
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
auto createKSP(MPI_Comm comm)
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
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS(MPI_Comm comm)
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
OpPlasticTools::CommonData CommonData
const double u0
inital vale on blocksets
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FaceElementForcesAndSourcesCore Ele
static char help[]
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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:112
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
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.
Post post-proc data at points from hash maps.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Variable Documentation

◆ help

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

Definition at line 10 of file reaction_diffusion.cpp.