v0.13.0
SCL-8: Radiation boundary conditions

Table of Contents

Note
Prerequisites of this tutorial include SCL-1: Poisson's equation (homogeneous BC)


Note
Intended learning outcome:
  • implementation of nonlinear boundary conditions
  • implementation of time stepping
  • implementation nonlineaar problem
  • implementation of axisymmetric problem

Source code

/**
* \file lesson6_radiation.cpp
* \example lesson6_radiation.cpp
*
* Using PipelineManager interface calculate the divergence of base functions,
* and integral of flux on the boundary. Since the h-div space is used, volume
* integral and boundary integral should give the same result.
*/
/* 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 <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
#include <BasicFiniteElements.hpp>
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
using OpDomainGradTimesVec = FormsIntegrators<DomainEleOp>::Assembly<
using OpBase = OpBaseImpl<PETSC, EdgeEleOp>;
// Units
// Temperature: Kelvins
// Length: 1 km
// Time: 1 s
constexpr double heat_conductivity = ((0.4+0.7)/2) * 1e3;
constexpr double emissivity = 1;
constexpr double boltzmann_constant = 5.670374419e-2;
constexpr double Beta = emissivity * boltzmann_constant;
constexpr double T_ambient = 2.7;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
static int integrationRule(int, int, int p_data) { return 2 * p_data; };
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode kspSolve();
MoFEMErrorCode postProcess();
MoFEMErrorCode checkResults();
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxGradVals;
struct OpRadiationLhs : public OpBase {
private:
boost::shared_ptr<VectorDouble> approxVals;
public:
OpRadiationLhs(boost::shared_ptr<VectorDouble> &approx_vals)
: OpBase("T", "T", OpBase::OPROWCOL), approxVals(approx_vals) {
this->sYmm = false;
}
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
struct OpRadiationRhs : public OpBase {
private:
boost::shared_ptr<VectorDouble> approxVals;
public:
OpRadiationRhs(boost::shared_ptr<VectorDouble> &approx_vals)
: OpBase("T", "T", OpBase::OPROW), approxVals(approx_vals) {}
MoFEMErrorCode iNtegrate(EntData &row_data);
};
struct OpFluxRhs : public OpBase {
private:
FTensor::Index<'i', 2> i; ///< summit Index
public:
OpFluxRhs() : OpBase("T", "T", OpBase::OPROW) {}
MoFEMErrorCode iNtegrate(EntData &row_data);
};
struct OpCalcSurfaceAverageTemperature : public EdgeEleOp {
private:
boost::shared_ptr<VectorDouble> approxVals;
double &sumTemperature;
double &surfaceArea;
public:
OpCalcSurfaceAverageTemperature(
boost::shared_ptr<VectorDouble> &approx_vals, double &sum_temp,
double &surf)
: EdgeEleOp("T", "T", OpBase::OPROW), approxVals(approx_vals),
sumTemperature(sum_temp), surfaceArea(surf) {}
MoFEMErrorCode doWork(int side, EntityType type,
};
};
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR bC();
CHKERR OPs();
CHKERR kspSolve();
CHKERR postProcess();
CHKERR checkResults();
}
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// Add field
CHKERR simple->addDomainField("T", H1, AINSWORTH_LEGENDRE_BASE, 1);
CHKERR simple->addBoundaryField("T", H1, AINSWORTH_LEGENDRE_BASE, 1);
constexpr int order = 3;
CHKERR simple->setFieldOrder("T", order);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Create common data]
approxVals = boost::make_shared<VectorDouble>();
approxGradVals = boost::make_shared<MatrixDouble>();
}
//! [Create common data]
//! [Boundary condition]
// Set initial values
auto set_initial_temperature = [&](VectorAdaptor &&field_data, double *xcoord,
double *ycoord, double *zcoord) {
field_data[0] = T_ambient;
};
FieldBlas *field_blas;
CHKERR mField.getInterface(field_blas);
CHKERR field_blas->setVertexDofs(set_initial_temperature, "T");
}
//! [Boundary condition]
//! [Push operators to pipeline]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto beta = [](const double r, const double, const double) {
return heat_conductivity * (2 * M_PI * r);
};
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("T", "T", beta));
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integrationRule);
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<2>("T", approxGradVals));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainGradTimesVec("T", approxGradVals, beta));
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpRadiationRhs(approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFluxRhs());
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integrationRule);
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpRadiationLhs(approxVals));
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integrationRule);
}
//! [Push operators to pipeline]
//! [Solve]
Simple *simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
CHKERR TSSetFromOptions(ts);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = smartCreateDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR TSSolve(ts, T);
CHKERR TSGetTime(ts, &ftime);
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetTimeStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
"steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
"%d, linits %d",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
//! [Solve]
//! [Postprocess results]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("T");
pipeline_mng->getDomainRhsFE() = post_proc_fe;
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateScalarFieldValues("T", approxVals));
double sum_temperature;
double surface_area;
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalcSurfaceAverageTemperature(approxVals, sum_temperature,
surface_area));
auto calc_surfcae_area_op = pipeline_mng->getOpBoundaryRhsPipeline().back();
sum_temperature = 0;
surface_area = 0;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_radiation.h5m");
MOFEM_LOG_C("EXAMPLE", Sev::inform, "Surface area %3.4e [km]", surface_area);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
"Average subsurface temperatute %3.4e [K]",
sum_temperature / surface_area);
}
//! [Postprocess results]
//! [Check results]
//! [Check results]
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::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
Simple *simple = m_field.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile("");
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
//! [Radiation Lhs]
EntData &col_data) {
// get element volume
const double vol = this->getMeasure();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_row_base = row_data.getFTensor0N();
// gat temperature at integration points
auto t_val = getFTensor0FromVec(*(approxVals));
// get coordinate at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// Cylinder radius
const double r_cylinder = t_coords(0);
// take into account Jacobean
const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
// loop over rows base functions
for (int rr = 0; rr != nbRows; ++rr) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
// loop over columns
for (int cc = 0; cc != nbCols; cc++) {
if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
locMat(rr, cc) += alpha * t_row_base * t_col_base * 4 * pow(t_val, 3);
}
++t_col_base;
}
++t_row_base;
}
++t_val;
++t_coords;
++t_w; // move to another integration weight
}
}
//! [Radiation Lhs]
//! [Radiation Lhs]
// get element volume
const double vol = getMeasure();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_row_base = row_data.getFTensor0N();
// gat temperature at integration points
auto t_val = getFTensor0FromVec(*(approxVals));
// get coordinate at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// Cylinder radius
const double r_cylinder = t_coords(0);
// take into account Jacobean
const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
// loop over rows base functions
for (int rr = 0; rr != nbRows; ++rr) {
if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
locF[rr] += alpha * t_row_base * (pow(t_val, 4) - pow(T_ambient, 4));
}
++t_row_base;
}
++t_coords;
++t_val;
++t_w; // move to another integration weight
}
}
//! [Radiation Lhs]
//! [Flux Rhs]
// get element volume
const double vol = getMeasure();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_row_base = row_data.getFTensor0N();
// get coordinate at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// // get time
const double time = getFEMethod()->ts_t;
// Look to https://doi.org/10.1016/j.icarus.2014.12.028s
constexpr double flux_p = -0.03e6;
constexpr double flux_c = -0.23e6;
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// Cylinder radius
const double r_cylinder = t_coords(0);
const double r = std::sqrt(t_coords(i) * t_coords(i));
const double s = std::abs(t_coords(1)) / r;
// take into account Jacobean
const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
// loop over rows base functions
for (int rr = 0; rr != nbRows; ++rr) {
locF[rr] += alpha * t_row_base * (s * flux_p + flux_c) * time;
++t_row_base;
}
++t_coords;
++t_w; // move to another integration weight
}
}
//! [Flux Rhs]
//! [Ave Temp]
int side, EntityType type, EntitiesFieldData::EntData &data) {
if (type == MBVERTEX) {
// get element volume
const double vol = getMeasure();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// gat temperature at integration points
auto t_val = getFTensor0FromVec(*(approxVals));
// get coordinate at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// number of integration pts
size_t nb_integration_pts = getGaussPts().size2();
// loop over integration points
for (auto gg = 0; gg != nb_integration_pts; ++gg) {
// Cylinder radius
const double r_cylinder = t_coords(0);
// take into account Jacobean
const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
sumTemperature += alpha * t_val;
surfaceArea += alpha;
++t_coords;
++t_val;
++t_w; // move to another integration weight
}
}
}
//! [Ave Temp]
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
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
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr double alpha
#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 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
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#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< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:126
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double r
rate factor
constexpr double Beta
Definition: radiation.cpp:52
constexpr double boltzmann_constant
Definition: radiation.cpp:51
constexpr double T_ambient
Definition: radiation.cpp:54
constexpr double emissivity
Definition: radiation.cpp:50
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Flux Rhs]
Definition: radiation.cpp:487
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
Definition: radiation.cpp:445
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
Definition: radiation.cpp:363
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
Definition: radiation.cpp:407
[Example]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
Definition: radiation.cpp:243
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:387
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:689
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:319
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)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
double heat_conductivity