v0.13.1
Loading...
Searching...
No Matches
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.
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
// 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) {}
private:
static int integrationRule(int, int, int p_data) { return 2 * p_data; };
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) {}
};
struct OpFluxRhs : public OpBase {
private:
FTensor::Index<'i', 2> i; ///< summit Index
public:
OpFluxRhs() : OpBase("T", "T", OpBase::OPROW) {}
};
struct OpCalcSurfaceAverageTemperature : public EdgeEleOp {
private:
boost::shared_ptr<VectorDouble> approxVals;
double &sumTemperature;
double &surfaceArea;
public:
boost::shared_ptr<VectorDouble> &approx_vals, double &sum_temp,
double &surf)
: EdgeEleOp("T", "T", OpBase::OPROW), approxVals(approx_vals),
sumTemperature(sum_temp), surfaceArea(surf) {}
};
};
}
//! [Set up problem]
// 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]
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 OpCalculateHOJac<2>(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("T", "T", beta));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHOJac<2>(jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpRadiationRhs(approxVals));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFluxRhs());
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpRadiationLhs(approxVals));
}
//! [Push operators to pipeline]
//! [Solve]
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]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
auto t_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("T", t_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"T", t_ptr}},
{}, {}, {}
)
);
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::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
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]
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:304
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
static char help[]
int main()
Definition: adol-c_atom.cpp:46
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
const double r
rate factor
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double Beta
Definition: radiation.cpp:40
constexpr double boltzmann_constant
Definition: radiation.cpp:39
constexpr double T_ambient
Definition: radiation.cpp:42
constexpr double emissivity
Definition: radiation.cpp:38
FTensor::Index< 'i', 3 > i
OpCalcSurfaceAverageTemperature(boost::shared_ptr< VectorDouble > &approx_vals, double &sum_temp, double &surf)
Definition: radiation.cpp:110
boost::shared_ptr< VectorDouble > approxVals
Definition: radiation.cpp:105
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Flux Rhs]
Definition: radiation.cpp:493
FTensor::Index< 'i', 2 > i
summit Index
Definition: radiation.cpp:94
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
Definition: radiation.cpp:451
boost::shared_ptr< VectorDouble > approxVals
Definition: radiation.cpp:68
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
Definition: radiation.cpp:369
OpRadiationLhs(boost::shared_ptr< VectorDouble > &approx_vals)
Definition: radiation.cpp:71
OpRadiationRhs(boost::shared_ptr< VectorDouble > &approx_vals)
Definition: radiation.cpp:85
boost::shared_ptr< VectorDouble > approxVals
Definition: radiation.cpp:82
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
Definition: radiation.cpp:413
[Example]
Definition: plastic.cpp:137
static int integrationRule(int, int, int p_data)
Definition: radiation.cpp:52
boost::shared_ptr< VectorDouble > approxVals
Definition: radiation.cpp:62
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
Definition: radiation.cpp:230
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:257
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:453
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:172
boost::shared_ptr< MatrixDouble > approxGradVals
Definition: radiation.cpp:63
MoFEM::Interface & mField
Definition: plastic.cpp:144
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:655
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:184
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:391
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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.
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:320
@ OPROW
operator doWork function is executed on FE rows
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double heat_conductivity