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.

                                 {
 
  // initialize petsc
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  try {
 
    // Create moab and mofem instances
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Simple interface
    Simple *simple_interface;
    CHKERR m_field.getInterface(simple_interface);
    CHKERR simple_interface->getOptions();
    CHKERR simple_interface->loadFile();
 
    int order = 4; ///< approximation order
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
 
    // add fields
    CHKERR simple_interface->addDomainField("u", H1, AINSWORTH_LEGENDRE_BASE,
                                            1);
    // set fields order
    CHKERR simple_interface->setFieldOrder("u", order);
    // setup problem
    CHKERR simple_interface->setUp();
 
    // Create common data structure
    boost::shared_ptr<CommonData> data(new CommonData());
    /// Alias pointers to data in common data structure
    auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
    auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
    auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
 
    // Create finite element instances to integrate the right-hand side of slow
    // and stiff vector, and the tangent left-hand side for stiff part.
    boost::shared_ptr<Ele> vol_ele_slow_rhs(new Ele(m_field));
    boost::shared_ptr<Ele> vol_ele_stiff_rhs(new Ele(m_field));
    boost::shared_ptr<Ele> vol_ele_stiff_lhs(new Ele(m_field));
 
    // Push operators to integrate the slow right-hand side vector
    vol_ele_slow_rhs->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("u", val_ptr));
    vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
 
    // PETSc IMAX and Explicit solver demans that g = M^-1 G is provided. So
    // when the slow right-hand side vector (G) is assembled is solved for g
    // vector.
    auto solve_for_g = [&]() {
      MoFEMFunctionBegin;
      if (vol_ele_slow_rhs->vecAssembleSwitch) {
        CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
                                   SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
                                 SCATTER_REVERSE);
        CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
        CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
        *vol_ele_slow_rhs->vecAssembleSwitch = false;
      }
      CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
                      vol_ele_slow_rhs->ts_F);
      MoFEMFunctionReturn(0);
    };
    // Add hook to the element to calculate g.
    vol_ele_slow_rhs->postProcessHook = solve_for_g;
 
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    // Add operators to calculate the stiff right-hand side
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpCalculateHOJacForFace(jac_ptr));
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpSetInvJacH1ForFace(inv_jac_ptr));
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
    vol_ele_stiff_rhs->getOpPtrVector().push_back(
        new OpAssembleStiffRhs<2>(data));
 
    // Add operators to calculate the stiff left-hand side
    vol_ele_stiff_lhs->getOpPtrVector().push_back(
        new OpCalculateHOJacForFace(jac_ptr));
    vol_ele_stiff_lhs->getOpPtrVector().push_back(
        new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    vol_ele_stiff_lhs->getOpPtrVector().push_back(
        new OpSetInvJacH1ForFace(inv_jac_ptr));
    vol_ele_stiff_lhs->getOpPtrVector().push_back(
        new OpAssembleStiffLhs<2>(data));
 
    // Set integration rules
    auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
    vol_ele_slow_rhs->getRuleHook = vol_rule;
    vol_ele_stiff_rhs->getRuleHook = vol_rule;
    vol_ele_stiff_lhs->getRuleHook = vol_rule;
 
    // Crate element for post-processing
    boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc =
        boost::shared_ptr<PostProcFaceOnRefinedMesh>(
            new PostProcFaceOnRefinedMesh(m_field));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    // Genarte post-processing mesh
    post_proc->generateReferenceElementMesh();
    // Postprocess only field values
    post_proc->addFieldValuesPostProc("u");
 
    // Get PETSc discrete manager
    auto dm = simple_interface->getDM();
 
    // Get surface entities form blockset, set initial values in those
    // blocksets. To keep it simple is assumed that inital values are on
    // blockset 1
    if (m_field.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
      Range inner_surface;
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          1, BLOCKSET, 2, inner_surface, true);
      if (!inner_surface.empty()) {
        Range inner_surface_verts;
        CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
        CHKERR m_field.getInterface<FieldBlas>()->setField(
            u0, MBVERTEX, inner_surface_verts, "u");
      }
    }
 
    // Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
    // conditions on that boundary.
    Range surface;
    CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
    Skinner skin(&m_field.get_moab());
    Range edges;
    CHKERR skin.find_skin(0, surface, false, edges);
    Range edges_part;
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &edges_part);
    Range edges_verts;
    CHKERR moab.get_connectivity(edges_part, edges_verts, false);
    // Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
    // the problem.
    CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
        simple_interface->getProblemName(), "u",
        unite(edges_verts, edges_part));
 
    // Create mass matrix, calculate and assemble
    CHKERR DMCreateMatrix_MoFEM(dm, data->M);
    CHKERR MatZeroEntries(data->M);
    boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
    vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
    CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
                                    vol_mass_ele);
    CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
 
    // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
    // M^-1G(t,u)
    data->ksp = createKSP(m_field.get_comm());
    CHKERR KSPSetOperators(data->ksp, data->M, data->M);
    CHKERR KSPSetFromOptions(data->ksp);
    CHKERR KSPSetUp(data->ksp);
 
    // Create and setup TS solver
    auto ts = createTS(m_field.get_comm());
    // Use IMEX solver, i.e. implicit/explicit solver
    CHKERR TSSetType(ts, TSARKIMEX);
    CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
 
    // Add element to calculate lhs of stiff part
    CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
                                 vol_ele_stiff_lhs, null, null);
    // Add element to calculate rhs of stiff part
    CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
                                 vol_ele_stiff_rhs, null, null);
    // Add element to calculate rhs of slow (nonlinear) part
    CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
                                   vol_ele_slow_rhs, null, null);
 
    // Add monitor to time solver
    boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
    CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
                               monitor_ptr, null, null);
 
    // Create solution vector
    SmartPetscObj<Vec> X;
    CHKERR DMCreateGlobalVector_MoFEM(dm, X);
    CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
 
    // Solve problem
    double ftime = 1;
    CHKERR TSSetDM(ts, dm);
    CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
    CHKERR TSSetSolution(ts, X);
    CHKERR TSSetFromOptions(ts);
    CHKERR TSSolve(ts, X);
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}

Variable Documentation

help

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

static

Examples

reaction_diffusion.cpp.

Definition at line 24 of file reaction_diffusion.cpp.