v0.10.0
heat_equation.cpp

Solve the time-dependent Heat Equation

\[ \begin{aligned} \frac{\partial u(\mathbf{x}, t)}{\partial t}-\Delta u(\mathbf{x}, t) &=f(\mathbf{x}, t) & & \forall \mathbf{x} \in \Omega, t \in(0, T), \\ u(\mathbf{x}, 0) &=u_{0}(\mathbf{x}) & & \forall \mathbf{x} \in \Omega, \\ u(\mathbf{x}, t) &=g(\mathbf{x}, t) & & \forall \mathbf{x} \in \partial \Omega, t \in(0, T). \end{aligned} \]

/**
* \file heat_equation.cpp
* \example heat_equation.cpp
*
* \brief Solve the time-dependent Heat Equation
\f[
\begin{aligned}
\frac{\partial u(\mathbf{x}, t)}{\partial t}-\Delta u(\mathbf{x}, t)
&=f(\mathbf{x}, t) & & \forall \mathbf{x} \in \Omega, t \in(0, T), \\
u(\mathbf{x}, 0) &=u_{0}(\mathbf{x}) & & \forall \mathbf{x} \in \Omega, \\
u(\mathbf{x}, t) &=g(\mathbf{x}, t) & & \forall \mathbf{x} \in \partial \Omega,
t \in(0, T). \end{aligned}
\f]
**/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <stdlib.h>
#include <cmath>
#include <BasicFiniteElements.hpp>
using namespace MoFEM;
using namespace HeatEquationOperators;
static char help[] = "...\n\n";
struct HeatEquation {
public:
// Declaration of the main function to run analysis
MoFEMErrorCode runProgram();
private:
// Declaration of other main functions called in runProgram()
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode setIntegrationRules();
MoFEMErrorCode initialCondition();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
// Function to calculate the Source term
static double sourceTermFunction(const double x, const double y,
const double z, const double t) {
return 0.1 * pow(M_E, -M_PI * M_PI * t) * sin(1. * M_PI * x) *
sin(2. * M_PI * y);
// return 0;
}
// Function to calculate the Boundary term
static double boundaryFunction(const double x, const double y, const double z,
const double t) {
return abs(0.1 * pow(M_E, -M_PI * M_PI * t) * sin(2. * M_PI * x) *
sin(3. * M_PI * y));
// return 0;
}
// initial value of u (initial condition)
const double initU = 0.1;
// Main interfaces
Simple *simpleInterface;
// mpi parallel communicator
MPI_Comm mpiComm;
// Number of processors
const int mpiRank;
// Discrete Manager and time-stepping solver using SmartPetscObj
SmartPetscObj<DM> dM;
SmartPetscObj<TS> tsSolver;
// Field name and approximation order
std::string domainField;
int oRder;
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<bool>> boundaryMarker;
// MoFEM working Pipelines for stiff part function and tangent (Jacobian)
boost::shared_ptr<FaceEle> stiffTangentLhsMatrixPipeline;
boost::shared_ptr<FaceEle> stiffFunctionRhsVectorPipeline;
boost::shared_ptr<EdgeEle> boundaryLhsMatrixPipeline;
boost::shared_ptr<EdgeEle> boundaryRhsVectorPipeline;
// Objects needed for solution updates in Newton's method
boost::shared_ptr<DataAtGaussPoints> previousUpdate;
boost::shared_ptr<VectorDouble> fieldValuePtr;
boost::shared_ptr<MatrixDouble> fieldGradPtr;
boost::shared_ptr<VectorDouble> fieldDotPtr;
// Object needed for postprocessing
boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcFace;
// Boundary entities marked for fieldsplit (block) solver - optional
Range boundaryEntitiesForFieldsplit;
};
: domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
mpiRank(mField.get_comm_rank()) {
stiffTangentLhsMatrixPipeline =
boost::shared_ptr<FaceEle>(new FaceEle(mField));
stiffFunctionRhsVectorPipeline =
boost::shared_ptr<FaceEle>(new FaceEle(mField));
boundaryLhsMatrixPipeline = boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
boundaryRhsVectorPipeline = boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
previousUpdate =
boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldValue);
fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
&previousUpdate->fieldGrad);
fieldDotPtr = boost::shared_ptr<VectorDouble>(previousUpdate,
&previousUpdate->fieldDot);
}
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
CHKERR simpleInterface->addBoundaryField(domainField, H1,
int oRder = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
}
auto stiff_tangent_rule = [](int, int, int p) -> int { return 2 * p; };
auto stiff_function_rule = [](int, int, int p) -> int { return 2 * p; };
stiffTangentLhsMatrixPipeline->getRuleHook = stiff_tangent_rule;
stiffFunctionRhsVectorPipeline->getRuleHook = stiff_function_rule;
auto boundary_lhs_rule = [](int, int, int p) -> int { return 2 * p; };
auto boundary_rhs_rule = [](int, int, int p) -> int { return 2 * p; };
boundaryLhsMatrixPipeline->getRuleHook = boundary_lhs_rule;
boundaryRhsVectorPipeline->getRuleHook = boundary_rhs_rule;
}
// Get surface entities form blockset, set initial values in those
// blocksets. To keep it simple, it is assumed that inital values are on
// blockset 1
if (mField.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
Range inner_surface;
CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
1, BLOCKSET, 2, inner_surface, true);
if (!inner_surface.empty()) {
Range inner_surface_verts;
CHKERR mField.get_moab().get_connectivity(inner_surface,
inner_surface_verts, false);
CHKERR mField.getInterface<FieldBlas>()->setField(
initU, MBVERTEX, inner_surface_verts, domainField);
}
}
}
auto get_ents_on_mesh_skin = [&]() {
Range boundary_entities;
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
boundary_entities, true);
}
}
// Add vertices
Range boundary_vertices;
CHKERR mField.get_moab().get_connectivity(boundary_entities,
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
// cerr << boundary_entities;
boundaryEntitiesForFieldsplit = boundary_entities;
// Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
// the problem.
// CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
// simpleInterface->getProblemName(), domainField, boundary_entities);
return boundary_entities;
};
auto mark_boundary_dofs = [&](Range &&skin_edges) {
auto problem_manager = mField.getInterface<ProblemsManager>();
auto marker_ptr = boost::make_shared<std::vector<bool>>();
problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
skin_edges, *marker_ptr);
return marker_ptr;
};
// Get global local vector of marked DOFs. Is global, since is set for all
// DOFs on processor. Is local since only DOFs on processor are in the
// vector. To access DOFs use local indices.
boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
}
{ // Push operators to the Pipeline that is responsible for calculating the
// Tangent (Jacobian) of function (LHS) of the stiff part
// Add default operators to calculate inverse of Jacobian (needed for
// implementation of 2D problem but not 3D ones)
stiffTangentLhsMatrixPipeline->getOpPtrVector().push_back(
stiffTangentLhsMatrixPipeline->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(invJac));
// Push operators for Jacobian of function (RHS) of the stiff part
stiffTangentLhsMatrixPipeline->getOpPtrVector().push_back(
new OpStiffTangentLHS(domainField, domainField, previousUpdate,
}
{ // Push operators to the Pipeline that is responsible for calculating the
// function (RHS) of the stiff part
// Add default operators to calculate inverse of Jacobian (needed for
// implementation of 2D problem but not 3D ones)
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(invJac));
// Add default operator to calculate field values at integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(domainField, fieldValuePtr));
// Add default operator to calculate field gradient at integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<2>(domainField, fieldGradPtr));
// Add default operator to calculate time derivative of field at
// integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
new OpCalculateScalarFieldValuesDot(domainField, fieldDotPtr));
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
new OpStiffFunctionRhs(domainField, sourceTermFunction, previousUpdate,
}
{ // Push operators to the Pipeline that is responsible for calculating the
// boundary condition
boundaryLhsMatrixPipeline->getOpPtrVector().push_back(
new OpBoundaryLhs(domainField, domainField));
// Add default operator to calculate field values at integration points
boundaryRhsVectorPipeline->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(domainField, fieldValuePtr));
boundaryRhsVectorPipeline->getOpPtrVector().push_back(
}
}
// Create and set Time Stepping solver
// Use fully implicit type (backward Euler) for Heat Equation
CHKERR TSSetType(tsSolver, TSBEULER);
// CHKERR TSSetType(tsSolver, TSARKIMEX);
// CHKERR TSARKIMEXSetType(tsSolver, TSARKIMEXA2);
// get Discrete Manager (SmartPetscObj)
dM = simpleInterface->getDM();
boost::shared_ptr<ForcesAndSourcesCore> null;
// Add element to calculate Jacobian of function dF/du (LHS) of stiff part
// and boundary term (LHS)
// Add element to calculate function F (RHS) of stiff part
// and boundary term (RHS)
// Add element to calculate function G (RHS) of slow (nonlinear) part
// Note: G(t,y) = 0 in the heat equation with fully implicit scheme
// CHKERR DMMoFEMTSSetRHSFunction(dM, simpleInterface->getDomainFEName(),
// vol_ele_slow_rhs, null, null);
{ // Set output of the results
// Create element for post-processing
postProcFace = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
// Generate post-processing mesh
postProcFace->generateReferenceElementMesh();
// Postprocess only field values
postProcFace->addFieldValuesPostProc(domainField);
// Add monitor to time solver
boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dM, postProcFace));
monitor_ptr, null, null);
}
// Create global solution vector
SmartPetscObj<Vec> global_solution;
// Scatter result data on the mesh
// CHKERR DMoFEMMeshToGlobalVector(dM, global_solution, INSERT_VALUES,
// SCATTER_REVERSE);
CHKERR DMoFEMMeshToLocalVector(dM, global_solution, INSERT_VALUES,
SCATTER_FORWARD);
// Solve problem
// double max_time = 1;
// double max_steps = 1000;
CHKERR TSSetDM(tsSolver, dM);
// CHKERR TSSetMaxTime(tsSolver, max_time);
// CHKERR TSSetMaxSteps(tsSolver, max_steps);
CHKERR TSSetSolution(tsSolver, global_solution);
CHKERR TSSetFromOptions(tsSolver);
CHKERR TSSolve(tsSolver, global_solution);
}
// Processes to set output results are integrated in solveSystem()
}
// outputResults();
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Error handling
try {
// Register MoFEM discrete manager in PETSc
DMType dm_name = "DMMOFEM";
// Create MOAB instance
moab::Core mb_instance; // mesh database
moab::Interface &moab = mb_instance; // mesh database interface
// Create MoFEM instance
MoFEM::Core core(moab); // finite element database
MoFEM::Interface &m_field = core; // finite element interface
// Run the main analysis
HeatEquation heat_problem(m_field);
CHKERR heat_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
static Index< 'p', 3 > p
std::string param_file
MatrixDouble invJac
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY > EdgeEle
@ ROW
Definition: definitions.h:192
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ BLOCKSET
Definition: definitions.h:217
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
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:719
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:445
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
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:772
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1026
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
int main(int argc, char *argv[])
static char help[]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:915
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1129
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
boost::shared_ptr< FaceEle > stiffFunctionRhsVectorPipeline
SmartPetscObj< DM > dM
boost::shared_ptr< EdgeEle > boundaryRhsVectorPipeline
SmartPetscObj< TS > tsSolver
boost::shared_ptr< std::vector< bool > > boundaryMarker
MoFEM::Interface & mField
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
boost::shared_ptr< FaceEle > stiffTangentLhsMatrixPipeline
boost::shared_ptr< EdgeEle > boundaryLhsMatrixPipeline
Range boundaryEntitiesForFieldsplit
boost::shared_ptr< VectorDouble > fieldDotPtr
boost::shared_ptr< VectorDouble > fieldValuePtr
HeatEquation(MoFEM::Interface &m_field)
Simple * simpleInterface
boost::shared_ptr< DataAtGaussPoints > previousUpdate
MoFEMErrorCode setupProblem()
static double sourceTermFunction(const double x, const double y, const double z, const double t)
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFace
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
const double initU
MoFEMErrorCode boundaryCondition()
MatrixDouble invJac
std::string domainField
MoFEMErrorCode runProgram()
MoFEMErrorCode initialCondition()
static double boundaryFunction(const double x, const double y, const double z, const double t)
boost::shared_ptr< MatrixDouble > fieldGradPtr
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
[Push operators to pipeline]
Postprocess on face.