v0.13.0
SCL-7: Wave equation
Note
Prerequisite of this tutorial is SCL-6: Heat equation


Note
Intended learning outcome:
  • practice of time-dependent problem in MoFEM
  • how to formulate and implement code in MoFEM using PETSc time-stepping solver TS
  • TS solver with implicit scheme

Introduction

Introduction

Plain program

The plain program for both the implementation of the UDOs (*.spp) and the main program (*.cpp) are as follows

Implementation of the main program (*.cpp)

/**
* \file wave_equation.cpp
* \example wave_equation.cpp
*
* \brief Solve the time-dependent Wave 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;
static char help[] = "...\n\n";
template <int DIM> struct ElementsAndOps {};
//! [Define dimension]
constexpr int SPACE_DIM = 2; //< Space dimension of problem, mesh
//! [Define dimension]
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainTimesScalarField = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainGradTimesVec = FormsIntegrators<DomainEleOp>::Assembly<
using OpBoundaryMass = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundaryTimeScalarField = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpBoundarySource = FormsIntegrators<BoundaryEleOp>::Assembly<
constexpr double wave_speed2 = 1;
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
MoFEMErrorCode preProcess() { return 0; }
MoFEMErrorCode operator()() { return 0; }
static constexpr int saveEveryNthStep = 1;
MoFEMErrorCode postProcess() {
if (ts_step % saveEveryNthStep == 0) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
};
struct WaveEquation {
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();
// Main interfaces
// Object to mark boundary entities for the assembling of domain elements
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
WaveEquation::WaveEquation(MoFEM::Interface &m_field) : mField(m_field) {}
auto *simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
}
auto *simple = mField.getInterface<Simple>();
CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
int order = 3;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setUp();
}
auto integration_rule = [](int o_row, int o_col, int approx_order) {
return 2 * approx_order;
};
auto *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
}
}
auto bc_mng = mField.getInterface<BcManager>();
auto *simple = mField.getInterface<Simple>();
CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ESSENTIAL",
"U", 0, 0);
auto &bc_map = bc_mng->getBcMapByBlockName();
boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
boundaryMarker->resize(b.second->bcMarkers.size(), 0);
for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
(*boundaryMarker)[i] |= b.second->bcMarkers[i];
}
}
}
}
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline.push_back(new OpSetHOWeightsOnFace());
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(new OpDomainGradGrad(
"U", "U", [](double, double, double) -> double { return 1; }));
auto get_c = [this](const double, const double, const double) {
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
return wave_speed2 * fe_domain_lhs->ts_aa;
};
pipeline.push_back(new OpDomainMass("U", "U", get_c));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
auto dot2_u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
"U", grad_u_at_gauss_pts));
pipeline.push_back(
new OpCalculateScalarFieldValuesDotDot("U", dot2_u_at_gauss_pts));
pipeline.push_back(new OpDomainGradTimesVec(
"U", grad_u_at_gauss_pts,
[](double, double, double) -> double { return 1; }));
pipeline.push_back(new OpDomainTimesScalarField(
"U", dot2_u_at_gauss_pts,
[](const double, const double, const double) { return wave_speed2; }));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_base_ops = [&](auto &pipeline) {};
auto add_lhs_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
pipeline.push_back(new OpBoundaryMass(
"U", "U", [](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_rhs_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto boundary_function = [&](const double x, const double y,
const double z) {
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_rhs = pipeline_mng->getBoundaryRhsFE();
const double t = fe_rhs->ts_t;
if ((t <= 0.5) && (x < 0.) && (y > -1. / 3) && (y < 1. / 3))
return sin(4 * M_PI * t);
else
return 0.;
};
pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
pipeline.push_back(new OpBoundaryTimeScalarField(
"U", u_at_gauss_pts,
[](const double, const double, const double) { return 1; }));
pipeline.push_back(new OpBoundarySource("U", boundary_function));
pipeline.push_back(new OpUnSetBc("U"));
};
auto pipeline_mng = mField.getInterface<PipelineManager>();
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
post_froc_fe->generateReferenceElementMesh();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_froc_fe->getOpPtrVector().push_back(
new OpCalculateHOJacForFace(jac_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
post_froc_fe->addFieldValuesPostProc("U");
return post_froc_fe;
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, create_post_process_element()));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
CHKERR KSPGetPC(ksp, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb = simple->getProblemName();
auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc); // boundary block
}
};
auto dm = simple->getDM();
auto D = smartCreateDMVector(dm);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
auto ts = pipeline_mng->createTSIM2();
CHKERR TSSetSolution(ts, D);
auto DD = smartVectorDuplicate(D);
CHKERR TS2SetSolution(ts, D, DD);
CHKERR set_time_monitor(dm, ts);
CHKERR TSSetFromOptions(ts);
CHKERR set_fieldsplit_preconditioner(ts);
CHKERR TSSetUp(ts);
CHKERR TSSolve(ts, NULL);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
// Processes to set output results are integrated in solveSystem()
}
}
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example")
// 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
WaveEquation heat_problem(m_field);
CHKERR heat_problem.runProgram();
}
// Finish work: cleaning memory, getting statistics, etc.
return 0;
}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main(int argc, char *argv[])
static char help[]
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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:544
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:43
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
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:994
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_TT > OpCalculateScalarFieldValuesDotDot
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
constexpr double t
plate stiffness
Definition: plate.cpp:76
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Postprocess on face.
WaveEquation(MoFEM::Interface &m_field)
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode setupProblem()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode initialCondition()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
MoFEMErrorCode runProgram()
MoFEMErrorCode setIntegrationRules()
constexpr double wave_speed2