v0.12.1
Classes | Public Types | Public Member Functions | Public Attributes | List of all members
MixTransport::UnsaturatedFlowElement Struct Reference

Implementation of operators, problem and finite elements for unsaturated flow. More...

#include <users_modules/tutorials/cor-0to1/src/UnsaturatedFlow.hpp>

Inheritance diagram for MixTransport::UnsaturatedFlowElement:
[legend]
Collaboration diagram for MixTransport::UnsaturatedFlowElement:
[legend]

Classes

struct  BcData
 Class storing information about boundary condition. More...
 
struct  FaceRule
 Set integration rule to boundary elements. More...
 
struct  MonitorPostProc
 
struct  OpDivTauU_HdivL2
 
struct  OpEvaluateInitiallHead
 
struct  OpIntegrateFluxes
 
struct  OpPostProcMaterial
 
struct  OpResidualFlux
 Assemble flux residual. More...
 
struct  OpResidualMass
 
struct  OpRhsBcOnValues
 Evaluate boundary condition at the boundary. More...
 
struct  OpTauDotSigma_HdivHdiv
 
struct  OpVDivSigma_L2Hdiv
 
struct  OpVU_L2L2
 
struct  postProcessVol
 Post proces method for volume element Assemble vectors and matrices and apply essential boundary conditions. More...
 
struct  preProcessVol
 Pre-peprocessing Set head pressute rate and get inital essential boundary conditions. More...
 
struct  VolRule
 Set integration rule to volume elements. More...
 

Public Types

typedef std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
 
typedef map< int, boost::shared_ptr< BcData > > BcMap
 

Public Member Functions

 UnsaturatedFlowElement (MoFEM::Interface &m_field)
 
 ~UnsaturatedFlowElement ()
 
virtual MoFEMErrorCode getMaterial (const EntityHandle ent, int &block_id) const
 For given element handle get material block Id. More...
 
MoFEMErrorCode getBcOnValues (const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
 Get value on boundary. More...
 
MoFEMErrorCode getBcOnFluxes (const EntityHandle ent, const double x, const double y, const double z, double &flux)
 essential (Neumann) boundary condition (set fluxes) More...
 
MoFEMErrorCode addFields (const std::string &values, const std::string &fluxes, const int order)
 add fields More...
 
MoFEMErrorCode addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
 add finite elements More...
 
MoFEMErrorCode buildProblem (Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
 Build problem. More...
 
MoFEMErrorCode setFiniteElements (ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
 Create finite element instances. More...
 
MoFEMErrorCode createMatrices ()
 Create vectors and matrices. More...
 
MoFEMErrorCode destroyMatrices ()
 Delete matrices and vector when no longer needed. More...
 
MoFEMErrorCode calculateEssentialBc ()
 Calculate boundary conditions for fluxes. More...
 
MoFEMErrorCode calculateInitialPc ()
 Calculate inital pressure head distribution. More...
 
MoFEMErrorCode solveProblem (bool set_initial_pc=true)
 solve problem More...
 
- Public Member Functions inherited from MixTransport::MixTransportElement
 MixTransportElement (MoFEM::Interface &m_field)
 construction of this data structure More...
 
virtual ~MixTransportElement ()
 destructor More...
 
MoFEMErrorCode getDirichletBCIndices (IS *is)
 get dof indices where essential boundary conditions are applied More...
 
virtual MoFEMErrorCode getSource (const EntityHandle ent, const double x, const double y, const double z, double &flux)
 set source term More...
 
virtual MoFEMErrorCode getResistivity (const EntityHandle ent, const double x, const double y, const double z, MatrixDouble3by3 &inv_k)
 natural (Dirichlet) boundary conditions (set values) More...
 
MoFEMErrorCode addFields (const std::string &values, const std::string &fluxes, const int order)
 Add fields to database. More...
 
MoFEMErrorCode addFiniteElements (const std::string &fluxes_name, const std::string &values_name)
 add finite elements More...
 
MoFEMErrorCode buildProblem (BitRefLevel &ref_level)
 Build problem. More...
 
MoFEMErrorCode postProc (const string out_file)
 Post process results. More...
 
MoFEMErrorCode createMatrices ()
 create matrices More...
 
MoFEMErrorCode solveLinearProblem ()
 solve problem More...
 
MoFEMErrorCode calculateResidual ()
 calculate residual More...
 
MoFEMErrorCode evaluateError ()
 Calculate error on elements. More...
 
MoFEMErrorCode destroyMatrices ()
 destroy matrices More...
 

Public Attributes

DM dM
 Discrete manager for unsaturated flow problem. More...
 
MaterialsDoubleMap dMatMap
 materials database More...
 
BcMap bcValueMap
 Store boundary condition for head capillary pressure. More...
 
EntityHandle lastEvalBcValEnt
 
int lastEvalBcBlockValId
 
BcMap bcFluxMap
 
EntityHandle lastEvalBcFluxEnt
 
int lastEvalBcBlockFluxId
 
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
 Elemnet to calculate essential bc. More...
 
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
 Face element apply natural bc. More...
 
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
 Calculate inital boundary conditions. More...
 
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
 Assemble residual vector. More...
 
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
 Assemble tangent matrix. More...
 
boost::shared_ptr< MethodForForceScalingscaleMethodFlux
 Method scaling fluxes. More...
 
boost::shared_ptr< MethodForForceScalingscaleMethodValue
 Method scaling values. More...
 
boost::shared_ptr< FEMethodtsMonitor
 
boost::shared_ptr< VectorDouble > headRateAtGaussPts
 Vector keeps head rate. More...
 
std::vector< int > bcVecIds
 
VectorDouble bcVecVals
 
VectorDouble vecValsOnBc
 
Vec D1
 Vector with inital head capillary pressure. More...
 
Vec ghostFlux
 Ghost Vector of integrated fluxes. More...
 
- Public Attributes inherited from MixTransport::MixTransportElement
MoFEM::InterfacemField
 
MyVolumeFE feVol
 Instance of volume element. More...
 
MyTriFE feTri
 Instance of surface element. More...
 
VectorDouble valuesAtGaussPts
 values at integration points on element More...
 
MatrixDouble valuesGradientAtGaussPts
 gradients at integration points on element More...
 
VectorDouble divergenceAtGaussPts
 divergence at integration points on element More...
 
MatrixDouble fluxesAtGaussPts
 fluxes at integration points on element More...
 
set< int > bcIndices
 
std::map< int, BlockDatasetOfBlocks
 maps block set id with appropriate BlockData More...
 
Vec D
 
Vec D0
 
Vec F
 
Mat Aij
 
map< double, EntityHandleerrorMap
 
double sumErrorFlux
 
double sumErrorDiv
 
double sumErrorJump
 

Detailed Description

Implementation of operators, problem and finite elements for unsaturated flow.

Definition at line 115 of file UnsaturatedFlow.hpp.

Member Typedef Documentation

◆ BcMap

typedef map<int, boost::shared_ptr<BcData> > MixTransport::UnsaturatedFlowElement::BcMap
Examples
UnsaturatedFlow.hpp.

Definition at line 165 of file UnsaturatedFlow.hpp.

◆ MaterialsDoubleMap

Examples
UnsaturatedFlow.hpp.

Definition at line 131 of file UnsaturatedFlow.hpp.

Constructor & Destructor Documentation

◆ UnsaturatedFlowElement()

MixTransport::UnsaturatedFlowElement::UnsaturatedFlowElement ( MoFEM::Interface m_field)
Examples
UnsaturatedFlow.hpp.

Definition at line 119 of file UnsaturatedFlow.hpp.

120  : MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
122  lastEvalBcBlockFluxId(-1) {}
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
DM dM
Discrete manager for unsaturated flow problem.

◆ ~UnsaturatedFlowElement()

MixTransport::UnsaturatedFlowElement::~UnsaturatedFlowElement ( )
Examples
UnsaturatedFlow.hpp.

Definition at line 124 of file UnsaturatedFlow.hpp.

124  {
125  if (dM != PETSC_NULL) {
126  CHKERR DMDestroy(&dM);
127  CHKERRABORT(PETSC_COMM_WORLD, ierr);
128  }
129  }
#define CHKERR
Inline error check.
Definition: definitions.h:548
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87

Member Function Documentation

◆ addFields()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::addFields ( const std::string &  values,
const std::string &  fluxes,
const int  order 
)

add fields

Examples
UnsaturatedFlow.hpp.

Definition at line 1147 of file UnsaturatedFlow.hpp.

1148  {
1150  // Fields
1153  CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1154  // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1155 
1156  // meshset consisting all entities in mesh
1157  EntityHandle root_set = mField.get_moab().get_root_set();
1158  // add entities to field
1159 
1161  if (it->getName().compare(0, 4, "SOIL") != 0)
1162  continue;
1163  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1164  MBTET, fluxes);
1165  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1166  MBTET, values);
1167  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1168  MBTET, values + "_t");
1169  // CHKERR mField.add_ents_to_field_by_type(
1170  // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1171  // );
1172  }
1173 
1174  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1175  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1176  CHKERR mField.set_field_order(root_set, MBTET, values, order);
1177  CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1178  // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1180  }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
#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.
MaterialsDoubleMap dMatMap
materials database
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0

◆ addFiniteElements()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::addFiniteElements ( const std::string &  fluxes_name,
const std::string &  values_name 
)

add finite elements

Examples
UnsaturatedFlow.hpp.

Definition at line 1183 of file UnsaturatedFlow.hpp.

1184  {
1186 
1187  // Define element "MIX". Note that this element will work with fluxes_name
1188  // and values_name. This reflect bilinear form for the problem
1197  values_name + "_t");
1198  // CHKERR
1199  // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1200 
1202  if (it->getName().compare(0, 4, "SOIL") != 0)
1203  continue;
1205  dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1206  }
1207 
1208  // Define element to integrate natural boundary conditions, i.e. set values.
1209  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1211  fluxes_name);
1213  fluxes_name);
1215  fluxes_name);
1217  values_name);
1218 
1220  if (it->getName().compare(0, 4, "HEAD") != 0)
1221  continue;
1223  bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1224  }
1225 
1226  // Define element to apply essential boundary conditions.
1227  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1229  fluxes_name);
1231  fluxes_name);
1233  fluxes_name);
1235  values_name);
1236 
1238  if (it->getName().compare(0, 4, "FLUX") != 0)
1239  continue;
1241  bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1242  }
1243 
1245  }
@ MF_ZERO
Definition: definitions.h:111
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
BcMap bcValueMap
Store boundary condition for head capillary pressure.

◆ buildProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::buildProblem ( Range  zero_flux_ents,
BitRefLevel  ref_level = BitRefLevel().set(0) 
)

Build problem.

Parameters
ref_levelmesh refinement on which mesh problem you like to built.
Returns
error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1252 of file UnsaturatedFlow.hpp.

1253  {
1255 
1256  // Build fields
1258  // Build finite elements
1260  CHKERR mField.build_finite_elements("MIX_BCFLUX");
1261  CHKERR mField.build_finite_elements("MIX_BCVALUE");
1262  // Build adjacencies of degrees of freedom and elements
1263  CHKERR mField.build_adjacencies(ref_level);
1264 
1265  // create DM instance
1266  CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1267  // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1268  CHKERR DMSetType(dM, "DMMOFEM");
1269  // mesh is portioned, each process keeps only part of problem
1270  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1271  // creates problem in DM
1272  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1273  // discretised problem creates square matrix (that makes some optimizations)
1274  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1275  // set DM options from command line
1276  CHKERR DMSetFromOptions(dM);
1277  // add finite elements
1278  CHKERR DMMoFEMAddElement(dM, "MIX");
1279  CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1280  CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1281  // constructor data structures
1282  CHKERR DMSetUp(dM);
1283 
1284  // remove zero flux dofs
1285  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1286  "MIX", "FLUXES", zero_flux_ents);
1287 
1288  PetscSection section;
1289  CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1290  CHKERR DMSetDefaultSection(dM, section);
1291  CHKERR DMSetDefaultGlobalSection(dM, section);
1292  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1293  CHKERR PetscSectionDestroy(&section);
1294 
1296  }
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1058
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:458
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ calculateEssentialBc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateEssentialBc ( )

Calculate boundary conditions for fluxes.

Returns
Error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1628 of file UnsaturatedFlow.hpp.

1628  {
1630  // clear vectors
1631  CHKERR VecZeroEntries(D0);
1632  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1633  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1634  // clear essential bc indices, it could have dofs from other mesh refinement
1635  bcIndices.clear();
1636  // set operator to calculate essential boundary conditions
1637  CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1638  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
1639  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
1640  CHKERR VecAssemblyBegin(D0);
1641  CHKERR VecAssemblyEnd(D0);
1642  double norm2D0;
1643  CHKERR VecNorm(D0, NORM_2, &norm2D0);
1644  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1645  PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1647  }
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:541
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.

◆ calculateInitialPc()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::calculateInitialPc ( )

Calculate inital pressure head distribution.

Returns
Error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1653 of file UnsaturatedFlow.hpp.

1653  {
1655  // clear vectors
1656  CHKERR VecZeroEntries(D1);
1657  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1658  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1659  // Calculate initial pressure head on each element
1661  // Assemble vector
1662  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
1663  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
1664  CHKERR VecAssemblyBegin(D1);
1665  CHKERR VecAssemblyEnd(D1);
1666  // Calculate norm
1667  double norm2D1;
1668  CHKERR VecNorm(D1, NORM_2, &norm2D1);
1669  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1670  PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1672  }
Vec D1
Vector with inital head capillary pressure.
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.

◆ createMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::createMatrices ( )

Create vectors and matrices.

Returns
Error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1594 of file UnsaturatedFlow.hpp.

1594  {
1596  CHKERR DMCreateMatrix(dM, &Aij);
1597  CHKERR DMCreateGlobalVector(dM, &D0);
1598  CHKERR VecDuplicate(D0, &D1);
1599  CHKERR VecDuplicate(D0, &D);
1600  CHKERR VecDuplicate(D0, &F);
1601  int ghosts[] = {0};
1602  int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1603  int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1604  CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1605  &ghostFlux);
1607  }
Vec ghostFlux
Ghost Vector of integrated fluxes.
virtual int get_comm_rank() const =0

◆ destroyMatrices()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::destroyMatrices ( )

Delete matrices and vector when no longer needed.

Returns
error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1613 of file UnsaturatedFlow.hpp.

1613  {
1615  CHKERR MatDestroy(&Aij);
1616  CHKERR VecDestroy(&D);
1617  CHKERR VecDestroy(&D0);
1618  CHKERR VecDestroy(&D1);
1619  CHKERR VecDestroy(&F);
1620  CHKERR VecDestroy(&ghostFlux);
1622  }

◆ getBcOnFluxes()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnFluxes ( const EntityHandle  ent,
const double  x,
const double  y,
const double  z,
double &  flux 
)
virtual

essential (Neumann) boundary condition (set fluxes)

Parameters
enthandle to finite element entity
xcoord
ycoord
zcoord
fluxreference to flux which is set by function
Returns
[description]

Reimplemented from MixTransport::MixTransportElement.

Examples
UnsaturatedFlow.hpp.

Definition at line 223 of file UnsaturatedFlow.hpp.

224  {
226  int block_id = -1;
227  if (lastEvalBcFluxEnt == ent) {
228  block_id = lastEvalBcBlockFluxId;
229  } else {
230  for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
231  it++) {
232  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
233  block_id = it->first;
234  }
235  }
236  lastEvalBcFluxEnt = ent;
237  lastEvalBcBlockFluxId = block_id;
238  }
239  if (block_id >= 0) {
240  if (bcFluxMap.at(block_id)->hookFun) {
241  flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
242  } else {
243  flux = bcFluxMap.at(block_id)->fixValue;
244  }
245  } else {
246  flux = 0;
247  }
249  }
double flux
impulse magnitude

◆ getBcOnValues()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getBcOnValues ( const EntityHandle  ent,
const int  gg,
const double  x,
const double  y,
const double  z,
double &  value 
)
virtual

Get value on boundary.

Parameters
ententity handle
ggnumber of integration point
xx-coordinate
yy-coordinate
zz-coordinate
valuereturned value
Returns
error code

Reimplemented from MixTransport::MixTransportElement.

Examples
UnsaturatedFlow.hpp.

Definition at line 181 of file UnsaturatedFlow.hpp.

183  {
185  int block_id = -1;
186  if (lastEvalBcValEnt == ent) {
187  block_id = lastEvalBcBlockValId;
188  } else {
189  for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
190  it++) {
191  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
192  block_id = it->first;
193  }
194  }
195  lastEvalBcValEnt = ent;
196  lastEvalBcBlockValId = block_id;
197  }
198  if (block_id >= 0) {
199  if (bcValueMap.at(block_id)->hookFun) {
200  value = bcValueMap.at(block_id)->hookFun(x, y, z);
201  } else {
202  value = bcValueMap.at(block_id)->fixValue;
203  }
204  } else {
205  value = 0;
206  }
208  }

◆ getMaterial()

virtual MoFEMErrorCode MixTransport::UnsaturatedFlowElement::getMaterial ( const EntityHandle  ent,
int &  block_id 
) const
virtual

For given element handle get material block Id.

Parameters
entfinite element entity handle
block_idreference to returned block id
Returns
error code
Examples
UnsaturatedFlow.hpp.

Definition at line 140 of file UnsaturatedFlow.hpp.

141  {
143  for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
144  mit != dMatMap.end(); mit++) {
145  if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
146  block_id = mit->first;
148  }
149  }
151  "Element not found, no material data");
153  }
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
virtual MPI_Comm & get_comm() const =0

◆ setFiniteElements()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::setFiniteElements ( ForcesAndSourcesCore::RuleHookFun  vol_rule = VolRule(),
ForcesAndSourcesCore::RuleHookFun  face_rule = FaceRule() 
)

Create finite element instances.

Parameters
vol_ruleintegration rule for volume element
face_ruleintegration rule for boundary element
Returns
error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1458 of file UnsaturatedFlow.hpp.

1459  {
1461 
1462  // create finite element instances
1463  feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1465  feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1467  feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1469  feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1471  feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1473  // set integration rule to elements
1474  feFaceBc->getRuleHook = face_rule;
1475  feFaceRhs->getRuleHook = face_rule;
1476  feVolInitialPc->getRuleHook = vol_rule;
1477  feVolLhs->getRuleHook = vol_rule;
1478  feVolRhs->getRuleHook = vol_rule;
1479  // set function hook for finite element postprocessing stage
1480  feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1481  feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1482  feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1483  feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1484 
1485  // Set Piola transform
1486  feFaceBc->getOpPtrVector().push_back(
1487  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1488  feFaceRhs->getOpPtrVector().push_back(
1489  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1490 
1491  // create method for setting history for fluxes on boundary
1492  scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
1493  new TimeForceScale("-flux_history", false));
1494 
1495  // create method for setting history for presure heads on boundary
1496  scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
1497  new TimeForceScale("-head_history", false));
1498 
1499 
1500  // Set operator to calculate essential boundary conditions
1501  feFaceBc->getOpPtrVector().push_back(
1502  new OpEvaluateBcOnFluxes(*this, "FLUXES"));
1503 
1504  // Set operator to calculate initial capillary pressure
1505  feVolInitialPc->getOpPtrVector().push_back(
1506  new OpEvaluateInitiallHead(*this, "VALUES"));
1507 
1508  // set residual face from Neumann terms, i.e. applied pressure
1509  feFaceRhs->getOpPtrVector().push_back(
1510  new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1511  // set residual finite element operators
1512  headRateAtGaussPts = boost::make_shared<VectorDouble>();
1513  // resAtGaussPts = boost::make_shared<VectorDouble>();
1514  feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1515  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1516  feVolRhs->getOpPtrVector().push_back(
1517  new OpValuesAtGaussPts(*this, "VALUES"));
1518  feVolRhs->getOpPtrVector().push_back(
1519  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1520  feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1521  feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1522  feVolRhs->getOpPtrVector().back().opType =
1523  ForcesAndSourcesCore::UserDataOperator::OPROW;
1524  // set tangent matrix finite element operators
1525  feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1526  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1527  // feVolLhs->getOpPtrVector().push_back(
1528  // new
1529  // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1530  // );
1531  feVolLhs->getOpPtrVector().push_back(
1532  new OpValuesAtGaussPts(*this, "VALUES"));
1533  feVolLhs->getOpPtrVector().push_back(
1534  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1535  feVolLhs->getOpPtrVector().push_back(
1536  new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1537  feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1538  feVolLhs->getOpPtrVector().push_back(
1539  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1540  feVolLhs->getOpPtrVector().push_back(
1541  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1542 
1543  // Adding finite elements to DM, time solver will ask for it to assemble
1544  // tangent matrices and residuals
1545  boost::shared_ptr<FEMethod> null;
1546  CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1547  CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1548  CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1549 
1550  // setting up post-processing
1551  boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
1553  CHKERR post_process->generateReferenceElementMesh();
1554  CHKERR post_process->addFieldValuesPostProc("VALUES");
1555  CHKERR post_process->addFieldValuesPostProc("VALUES_t");
1556  // CHKERR post_process->addFieldValuesPostProc("FLUXES_residual");
1557  CHKERR post_process->addFieldValuesPostProc("FLUXES");
1558  post_process->getOpPtrVector().push_back(
1559  new OpValuesAtGaussPts(*this, "VALUES"));
1560  post_process->getOpPtrVector().push_back(
1561  new OpPostProcMaterial(*this, post_process->postProcMesh,
1562  post_process->mapGaussPts, "VALUES"));
1563 
1564  // Integrate fluxes on boundary
1565  boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1566  flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1568  flux_integrate->getOpPtrVector().push_back(
1569  new OpIntegrateFluxes(*this, "FLUXES"));
1570  int frequency = 1;
1571  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1572  CHKERR PetscOptionsInt(
1573  "-how_often_output",
1574  "frequency how often results are dumped on hard disk", "", frequency,
1575  &frequency, NULL);
1576  ierr = PetscOptionsEnd();
1577  CHKERRG(ierr);
1578 
1579  tsMonitor = boost::shared_ptr<FEMethod>(
1580  new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1581  TsCtx *ts_ctx;
1582  DMMoFEMGetTsCtx(dM, &ts_ctx);
1583  ts_ctx->getPostProcessMonitor().push_back(tsMonitor);
1585  }
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
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:755
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1077
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:808
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
boost::shared_ptr< FEMethod > tsMonitor
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
Post processing.
Force scale operator for reading two columns.

◆ solveProblem()

MoFEMErrorCode MixTransport::UnsaturatedFlowElement::solveProblem ( bool  set_initial_pc = true)

solve problem

Returns
error code
Examples
UnsaturatedFlow.hpp.

Definition at line 1678 of file UnsaturatedFlow.hpp.

1678  {
1680  if (set_initial_pc) {
1681  // Set initial head
1682  CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1683  }
1684 
1685  // Initiate vector from data on the mesh
1686  CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1687 
1688  // Create time solver
1689  TS ts;
1690  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1691  // Use backward Euler method
1692  CHKERR TSSetType(ts, TSBEULER);
1693  // Set final time
1694  double ftime = 1;
1695  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1696  // Setup solver from commabd line
1697  CHKERR TSSetFromOptions(ts);
1698  // Set DM to TS
1699  CHKERR TSSetDM(ts, dM);
1700 #if PETSC_VERSION_GE(3, 7, 0)
1701  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1702 #endif
1703  // Set-up monitor
1704  TsCtx *ts_ctx;
1705  DMMoFEMGetTsCtx(dM, &ts_ctx);
1706  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1707 
1708  // This add SNES monitor, to show error by fields. It is dirty trick
1709  // to add monitor, so code is hidden from doxygen
1710  CHKERR TSSetSolution(ts, D);
1711  CHKERR TSSetUp(ts);
1712  SNES snes;
1713  CHKERR TSGetSNES(ts, &snes);
1714 
1715 #if PETSC_VERSION_GE(3, 7, 0)
1716  {
1717  PetscViewerAndFormat *vf;
1718  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1719  PETSC_VIEWER_DEFAULT, &vf);
1720  CHKERR SNESMonitorSet(
1721  snes,
1722  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1723  void *))SNESMonitorFields,
1724  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1725  }
1726 #else
1727  {
1728  CHKERR SNESMonitorSet(snes,
1729  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1730  void *))SNESMonitorFields,
1731  0, 0);
1732  }
1733 #endif
1734 
1735  CHKERR TSSolve(ts, D);
1736 
1737  // Get statistic form TS and print it
1738  CHKERR TSGetTime(ts, &ftime);
1739  PetscInt steps, snesfails, rejects, nonlinits, linits;
1740  CHKERR TSGetTimeStepNumber(ts, &steps);
1741  CHKERR TSGetSNESFailures(ts, &snesfails);
1742  CHKERR TSGetStepRejections(ts, &rejects);
1743  CHKERR TSGetSNESIterations(ts, &nonlinits);
1744  CHKERR TSGetKSPIterations(ts, &linits);
1745  PetscPrintf(PETSC_COMM_WORLD,
1746  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1747  "%D, linits %D\n",
1748  steps, rejects, snesfails, ftime, nonlinits, linits);
1749 
1751  }
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:478
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:209

Member Data Documentation

◆ bcFluxMap

BcMap MixTransport::UnsaturatedFlowElement::bcFluxMap
Examples
UnsaturatedFlow.hpp.

Definition at line 210 of file UnsaturatedFlow.hpp.

◆ bcValueMap

BcMap MixTransport::UnsaturatedFlowElement::bcValueMap

Store boundary condition for head capillary pressure.

Examples
UnsaturatedFlow.hpp.

Definition at line 166 of file UnsaturatedFlow.hpp.

◆ bcVecIds

std::vector<int> MixTransport::UnsaturatedFlowElement::bcVecIds
Examples
UnsaturatedFlow.hpp.

Definition at line 1335 of file UnsaturatedFlow.hpp.

◆ bcVecVals

VectorDouble MixTransport::UnsaturatedFlowElement::bcVecVals
Examples
UnsaturatedFlow.hpp.

Definition at line 1336 of file UnsaturatedFlow.hpp.

◆ D1

Vec MixTransport::UnsaturatedFlowElement::D1

Vector with inital head capillary pressure.

Examples
UnsaturatedFlow.hpp.

Definition at line 1587 of file UnsaturatedFlow.hpp.

◆ dM

DM MixTransport::UnsaturatedFlowElement::dM

Discrete manager for unsaturated flow problem.

Examples
UnsaturatedFlow.hpp.

Definition at line 117 of file UnsaturatedFlow.hpp.

◆ dMatMap

MaterialsDoubleMap MixTransport::UnsaturatedFlowElement::dMatMap

materials database

Examples
UnsaturatedFlow.hpp.

Definition at line 132 of file UnsaturatedFlow.hpp.

◆ feFaceBc

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feFaceBc

Elemnet to calculate essential bc.

Examples
UnsaturatedFlow.hpp.

Definition at line 1299 of file UnsaturatedFlow.hpp.

◆ feFaceRhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feFaceRhs

Face element apply natural bc.

Examples
UnsaturatedFlow.hpp.

Definition at line 1301 of file UnsaturatedFlow.hpp.

◆ feVolInitialPc

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolInitialPc

Calculate inital boundary conditions.

Examples
UnsaturatedFlow.hpp.

Definition at line 1303 of file UnsaturatedFlow.hpp.

◆ feVolLhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolLhs

Assemble tangent matrix.

Examples
UnsaturatedFlow.hpp.

Definition at line 1306 of file UnsaturatedFlow.hpp.

◆ feVolRhs

boost::shared_ptr<ForcesAndSourcesCore> MixTransport::UnsaturatedFlowElement::feVolRhs

Assemble residual vector.

Examples
UnsaturatedFlow.hpp.

Definition at line 1305 of file UnsaturatedFlow.hpp.

◆ ghostFlux

Vec MixTransport::UnsaturatedFlowElement::ghostFlux

Ghost Vector of integrated fluxes.

Examples
UnsaturatedFlow.hpp.

Definition at line 1588 of file UnsaturatedFlow.hpp.

◆ headRateAtGaussPts

boost::shared_ptr<VectorDouble> MixTransport::UnsaturatedFlowElement::headRateAtGaussPts

Vector keeps head rate.

Examples
UnsaturatedFlow.hpp.

Definition at line 1315 of file UnsaturatedFlow.hpp.

◆ lastEvalBcBlockFluxId

int MixTransport::UnsaturatedFlowElement::lastEvalBcBlockFluxId
Examples
UnsaturatedFlow.hpp.

Definition at line 212 of file UnsaturatedFlow.hpp.

◆ lastEvalBcBlockValId

int MixTransport::UnsaturatedFlowElement::lastEvalBcBlockValId
Examples
UnsaturatedFlow.hpp.

Definition at line 169 of file UnsaturatedFlow.hpp.

◆ lastEvalBcFluxEnt

EntityHandle MixTransport::UnsaturatedFlowElement::lastEvalBcFluxEnt
Examples
UnsaturatedFlow.hpp.

Definition at line 211 of file UnsaturatedFlow.hpp.

◆ lastEvalBcValEnt

EntityHandle MixTransport::UnsaturatedFlowElement::lastEvalBcValEnt
Examples
UnsaturatedFlow.hpp.

Definition at line 168 of file UnsaturatedFlow.hpp.

◆ scaleMethodFlux

boost::shared_ptr<MethodForForceScaling> MixTransport::UnsaturatedFlowElement::scaleMethodFlux

Method scaling fluxes.

Examples
UnsaturatedFlow.hpp.

Definition at line 1308 of file UnsaturatedFlow.hpp.

◆ scaleMethodValue

boost::shared_ptr<MethodForForceScaling> MixTransport::UnsaturatedFlowElement::scaleMethodValue

Method scaling values.

Examples
UnsaturatedFlow.hpp.

Definition at line 1310 of file UnsaturatedFlow.hpp.

◆ tsMonitor

boost::shared_ptr<FEMethod> MixTransport::UnsaturatedFlowElement::tsMonitor

Element used by TS monitor to postprocess results at time step

Examples
UnsaturatedFlow.hpp.

Definition at line 1311 of file UnsaturatedFlow.hpp.

◆ vecValsOnBc

VectorDouble MixTransport::UnsaturatedFlowElement::vecValsOnBc
Examples
UnsaturatedFlow.hpp.

Definition at line 1336 of file UnsaturatedFlow.hpp.


The documentation for this struct was generated from the following file: