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>
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
// 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);
}
}
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
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);
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>>();
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(
// Push operators for Jacobian of function (RHS) of the stiff part
stiffTangentLhsMatrixPipeline->getOpPtrVector().push_back(
}
{ // 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(
// Add default operator to calculate field values at integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field gradient at integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
// Add default operator to calculate time derivative of field at
// integration points
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
stiffFunctionRhsVectorPipeline->getOpPtrVector().push_back(
}
{ // Push operators to the Pipeline that is responsible for calculating the
// boundary condition
boundaryLhsMatrixPipeline->getOpPtrVector().push_back(
// Add default operator to calculate field values at integration points
boundaryRhsVectorPipeline->getOpPtrVector().push_back(
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)
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;
}