v0.14.0
Loading...
Searching...
No Matches
UnsaturatedFlow.hpp
/** \file UnsaturatedFlow.hpp
* \brief Mix implementation of transport element
* \example UnsaturatedFlow.hpp
*
* \ingroup mofem_mix_transport_elem
*
*/
#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
namespace MixTransport {
using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
/**
* \brief Generic material model for unsaturated water transport
*
* \note This is abstact class, no physical material model is implemented
* here.
*/
struct GenericMaterial {
virtual ~GenericMaterial() {}
static double ePsilon0; ///< Regularization parameter
static double ePsilon1; ///< Regularization parameter
static double scaleZ; ///< Scale z direction
double sCale; ///< Scale time dependent eq.
double h; ///< hydraulic head
double h_t; ///< rate of hydraulic head
// double h_flux_residual; ///< residual at point
double K; ///< Hydraulic conductivity [L/s]
double diffK; ///< Derivative of hydraulic conductivity [L/s * L^2/F]
double C; ///< Capacity [S^2/L^2]
double diffC; ///< Derivative of capacity [S^2/L^2 * L^2/F ]
double tHeta; ///< Water content
double Se; ///< Effective saturation
Range tEts; ///< Elements with this material
double x, y, z; ///< in meters (L)
/**
* \brief Initialize head
* @return value of head
*/
virtual double initialPcEval() const = 0;
virtual void printMatParameters(const int id,
const std::string &prefix) const = 0;
virtual MoFEMErrorCode calK() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate hydraulic conductivity");
}
virtual MoFEMErrorCode calDiffK() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate derivative of hydraulic "
"conductivity");
}
virtual MoFEMErrorCode calC() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
virtual MoFEMErrorCode calDiffC() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
virtual MoFEMErrorCode calTheta() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
virtual MoFEMErrorCode calSe() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
};
/**
* \brief Implementation of operators, problem and finite elements for
* unsaturated flow
*/
struct UnsaturatedFlowElement : public MixTransportElement {
DM dM; ///< Discrete manager for unsaturated flow problem
UnsaturatedFlowElement(MoFEM::Interface &m_field)
: MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
if (dM != PETSC_NULL) {
CHKERR DMDestroy(&dM);
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
}
typedef std::map<int, boost::shared_ptr<GenericMaterial>> MaterialsDoubleMap;
MaterialsDoubleMap dMatMap; ///< materials database
/**
* \brief For given element handle get material block Id
* @param ent finite element entity handle
* @param block_id reference to returned block id
* @return error code
*/
virtual MoFEMErrorCode getMaterial(const EntityHandle ent,
int &block_id) const {
for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
mit != dMatMap.end(); mit++) {
if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
block_id = mit->first;
}
}
"Element not found, no material data");
}
/**
* \brief Class storing information about boundary condition
*/
struct BcData {
double fixValue;
boost::function<double(const double x, const double y, const double z)>
BcData() : hookFun(NULL) {}
};
typedef map<int, boost::shared_ptr<BcData>> BcMap;
BcMap bcValueMap; ///< Store boundary condition for head capillary pressure
/**
* \brief Get value on boundary
* @param ent entity handle
* @param gg number of integration point
* @param x x-coordinate
* @param y y-coordinate
* @param z z-coordinate
* @param value returned value
* @return error code
*/
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
const double x, const double y, const double z,
double &value) {
int block_id = -1;
if (lastEvalBcValEnt == ent) {
} else {
for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
if (bcValueMap.at(block_id)->hookFun) {
value = bcValueMap.at(block_id)->hookFun(x, y, z);
} else {
value = bcValueMap.at(block_id)->fixValue;
}
} else {
value = 0;
}
}
/**
* \brief essential (Neumann) boundary condition (set fluxes)
* @param ent handle to finite element entity
* @param x coord
* @param y coord
* @param z coord
* @param flux reference to flux which is set by function
* @return [description]
*/
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
const double y, const double z, double &flux) {
int block_id = -1;
if (lastEvalBcFluxEnt == ent) {
} else {
for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
if (bcFluxMap.at(block_id)->hookFun) {
flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
} else {
flux = bcFluxMap.at(block_id)->fixValue;
}
} else {
flux = 0;
}
}
/**
* \brief Evaluate boundary condition at the boundary
*/
struct OpRhsBcOnValues
UnsaturatedFlowElement &cTx;
boost::shared_ptr<MethodForForceScaling> valueScale;
/**
* \brief Constructor
*/
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
boost::shared_ptr<MethodForForceScaling> &value_scale)
fluxes_name, UserDataOperator::OPROW),
cTx(ctx), valueScale(value_scale) {}
VectorDouble nF; ///< Vector of residuals
/**
* \brief Integrate boundary condition
* @param side local index of entity
* @param type type of entity
* @param data data on entity
* @return error code
*/
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
if (data.getFieldData().size() == 0)
// Get EntityHandle of the finite element
// Resize and clear vector
nF.resize(data.getIndices().size());
nF.clear();
// Get number of integration points
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double x, y, z;
x = getCoordsAtGaussPts()(gg, 0);
y = getCoordsAtGaussPts()(gg, 1);
z = getCoordsAtGaussPts()(gg, 2);
double value;
// get value of boundary condition
CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
const double w = getGaussPts()(2, gg) * 0.5;
const double beta = w * (value - z);
noalias(nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
}
// Scale vector if history evaluating method is given
Vec f = getFEMethod()->ts_F;
if (valueScale) {
}
// Assemble vector
CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
&nF[0], ADD_VALUES);
}
};
/**
* \brief Assemble flux residual
*/
struct OpResidualFlux
UnsaturatedFlowElement &cTx;
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
flux_name, UserDataOperator::OPROW),
cTx(ctx) {}
VectorDouble divVec, nF;
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs, false);
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// block_data->printMatParameters(block_id,"Read material");
// Get base function
auto t_n_hdiv = data.getFTensor1N<3>();
// Get pressure
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// Get flux
auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get integration weight
// Get volume
double vol = getVolume();
auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
divVec.resize(nb_dofs, false);
// Get material parameters
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get divergence
for (auto &v : divVec) {
v = t_base_diff_hdiv(i, i);
++t_base_diff_hdiv;
}
const double alpha = t_w * vol;
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calK();
const double K = block_data->K;
const double scaleZ = block_data->scaleZ;
const double z = t_coords(2); /// z-coordinate at Gauss pt
// Calculate pressure gradient
noalias(nF) -= alpha * (t_h - z * scaleZ) * divVec;
// Calculate presure gradient from flux
FTensor::Tensor0<double *> t_nf(&*nF.begin());
for (int rr = 0; rr != nb_dofs; rr++) {
t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
++t_n_hdiv; // move to next base function
++t_nf; // move to next element in vector
}
++t_h; // move to next integration point
++t_flux;
++t_coords;
++t_w;
}
// Assemble residual
CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
&*data.getIndices().begin(), &*nF.begin(),
ADD_VALUES);
}
};
/**
* Assemble mass residual
*/
struct OpResidualMass
UnsaturatedFlowElement &cTx;
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
val_name, UserDataOperator::OPROW),
cTx(ctx) {}
VectorDouble nF;
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
// Resize local element vector
nF.resize(nb_dofs, false);
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// Get pressure rate
auto t_h_t = getFTensor0FromVec(*cTx.headRateAtGaussPts);
// Flux divergence
auto t_div_flux = getFTensor0FromVec(cTx.divergenceAtGaussPts);
// Get integration weight
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Scale eq.
const double scale = block_data->sCale;
// Get volume
const double vol = getVolume();
// Get number of integration points
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
const double alpha = t_w * vol * scale;
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calC();
const double C = block_data->C;
// Calculate flux conservation
noalias(nF) += (alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
++t_h;
++t_h_t;
++t_div_flux;
++t_coords;
++t_w;
}
// Assemble local vector
Vec f = getFEMethod()->ts_F;
CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
ADD_VALUES);
}
};
struct OpTauDotSigma_HdivHdiv
UnsaturatedFlowElement &cTx;
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx,
const std::string flux_name)
flux_name, flux_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
MatrixDouble nN, transNN;
/**
* \brief Assemble matrix
* @param row_side local index of row entity on element
* @param col_side local index of col entity on element
* @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param row_data data for row
* @param col_data data for col
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const int nb_row = row_data.getIndices().size();
const int nb_col = col_data.getIndices().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col, false);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get base functions
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
// Get integration weight
// Get volume
const double vol = getVolume();
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calK();
const double K = block_data->K;
// get integration weight and multiply by element volume
const double alpha = t_w * vol;
const double beta = alpha * (1 / K);
FTensor::Tensor0<double *> t_a(&*nN.data().begin());
for (int kk = 0; kk != nb_row; kk++) {
auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
for (int ll = 0; ll != nb_col; ll++) {
t_a += beta * (t_n_hdiv_row(j) * t_n_hdiv_col(j));
++t_n_hdiv_col;
++t_a;
}
++t_n_hdiv_row;
}
++t_coords;
++t_h;
++t_w;
}
Mat a = getFEMethod()->ts_B;
CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
&*col_data.getIndices().begin(), &*nN.data().begin(),
ADD_VALUES);
// matrix is symmetric, assemble other part
if (row_side != col_side || row_type != col_type) {
transNN.resize(nb_col, nb_row);
noalias(transNN) = trans(nN);
CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
&*row_data.getIndices().begin(),
&*transNN.data().begin(), ADD_VALUES);
}
}
};
struct OpVU_L2L2
UnsaturatedFlowElement &cTx;
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
value_name, value_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
MatrixDouble nN;
/**
* \brief Assemble matrix
* @param row_side local index of row entity on element
* @param col_side local index of col entity on element
* @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param row_data data for row
* @param col_data data for col
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
int nb_row = row_data.getIndices().size();
int nb_col = col_data.getIndices().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col, false);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// // Get pressure
// auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
// Get pressure rate
auto t_h_t = getFTensor0FromVec(*cTx.headRateAtGaussPts);
// Get integration weight
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Scale eq.
const double scale = block_data->sCale;
// Time step factor
double ts_a = getFEMethod()->ts_a;
// get volume
const double vol = getVolume();
int nb_gauss_pts = row_data.getN().size1();
// get base functions on rows
auto t_n_row = row_data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get integration weight and multiply by element volume
double alpha = t_w * vol * scale;
// evaluate material model at integration points
// to calculate capacity and tangent of capacity term
block_data->h = t_h;
block_data->h_t = t_h_t;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calC();
CHKERR block_data->calDiffC();
const double C = block_data->C;
const double diffC = block_data->diffC;
// assemble local entity tangent matrix
&*nN.data().begin());
// iterate base functions on rows
for (int kk = 0; kk != nb_row; kk++) {
// get first base function on column at integration point gg
auto t_n_col = col_data.getFTensor0N(gg, 0);
// iterate base functions on columns
for (int ll = 0; ll != nb_col; ll++) {
// assemble elements of local matrix
t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
++t_n_col; // move to next base function on column
++t_a; // move to next element in local tangent matrix
}
++t_n_row; // move to next base function on row
}
++t_w; // move to next integration weight
++t_coords; // move to next coordinate at integration point
++t_h; // move to next capillary head at integration point
// ++t_flux_residual;
++t_h_t; // move to next capillary head rate at integration point
}
Mat a = getFEMethod()->ts_B;
CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &*nN.data().begin(),
ADD_VALUES);
}
};
struct OpVDivSigma_L2Hdiv
UnsaturatedFlowElement &cTx;
/**
* \brief Constructor
*/
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx,
const std::string &val_name_row,
const std::string &flux_name_col)
val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
MatrixDouble nN;
VectorDouble divVec;
/**
* \brief Do calculations
* @param row_side local index of entity on row
* @param col_side local index of entity on column
* @param row_type type of row entity
* @param col_type type of col entity
* @param row_data row data structure carrying information about base
* functions, DOFs indices, etc.
* @param col_data column data structure carrying information about base
* functions, DOFs indices, etc.
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
int nb_row = row_data.getFieldData().size();
int nb_col = col_data.getFieldData().size();
if (nb_row == 0)
if (nb_col == 0)
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
nN.resize(nb_row, nb_col, false);
divVec.resize(nb_col, false);
nN.clear();
// take derivatives of base functions
auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
// Scale eq.
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double alpha = getGaussPts()(3, gg) * getVolume() * scale;
for (auto &v : divVec) {
v = t_base_diff_hdiv(i, i);
++t_base_diff_hdiv;
}
noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
}
CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
}
};
struct OpDivTauU_HdivL2
UnsaturatedFlowElement &cTx;
/**
* \brief Constructor
*/
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx,
const std::string &flux_name_col,
const std::string &val_name_row)
flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
MatrixDouble nN;
VectorDouble divVec;
/**
* \brief Do calculations
* @param row_side local index of entity on row
* @param col_side local index of entity on column
* @param row_type type of row entity
* @param col_type type of col entity
* @param row_data row data structure carrying information about base
* functions, DOFs indices, etc.
* @param col_data column data structure carrying information about base
* functions, DOFs indices, etc.
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
int nb_row = row_data.getFieldData().size();
int nb_col = col_data.getFieldData().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col, false);
divVec.resize(nb_row, false);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// Get flux
auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get integration weight
// Get base function
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
// Get derivatives of base functions
auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
// Get volume
double vol = getVolume();
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calK();
CHKERR block_data->calDiffK();
const double K = block_data->K;
// const double z = t_coords(2);
const double KK = K * K;
const double diffK = block_data->diffK;
double alpha = t_w * vol;
for (auto &v : divVec) {
v = t_base_diff_hdiv(i, i);
++t_base_diff_hdiv;
}
noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
FTensor::Tensor0<double *> t_a(&*nN.data().begin());
for (int rr = 0; rr != nb_row; rr++) {
double beta = alpha * (-diffK / KK) * (t_n_hdiv_row(i) * t_flux(i));
auto t_n_col = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col; cc++) {
t_a += beta * t_n_col;
++t_n_col;
++t_a;
}
++t_n_hdiv_row;
}
++t_w;
++t_coords;
++t_h;
++t_flux;
}
CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
}
};
struct OpEvaluateInitiallHead
UnsaturatedFlowElement &cTx;
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx,
const std::string &val_name)
val_name, UserDataOperator::OPROW),
cTx(ctx) {}
MatrixDouble nN;
VectorDouble nF;
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
if (data.getFieldData().size() == 0)
auto nb_dofs = data.getFieldData().size();
if (nb_dofs != data.getN().size2()) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"wrong number of dofs");
}
nN.resize(nb_dofs, nb_dofs, false);
nF.resize(nb_dofs, false);
nN.clear();
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// loop over integration points
auto nb_gauss_pts = getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; gg++) {
// get coordinates at integration point
block_data->x = getCoordsAtGaussPts()(gg, 0);
block_data->y = getCoordsAtGaussPts()(gg, 1);
block_data->z = getCoordsAtGaussPts()(gg, 2);
// get weight for integration rule
double alpha = getGaussPts()(2, gg);
nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
nF += alpha * block_data->initialPcEval() * data.getN(gg);
}
// factor matrix
// solve local problem
cholesky_solve(nN, nF, ublas::lower());
// set solution to vector
CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
&*nF.begin(), INSERT_VALUES);
}
};
struct OpIntegrateFluxes
UnsaturatedFlowElement &cTx;
/**
* \brief Constructor
*/
OpIntegrateFluxes(UnsaturatedFlowElement &ctx,
const std::string fluxes_name)
fluxes_name, UserDataOperator::OPROW),
cTx(ctx) {}
/**
* \brief Integrate boundary condition
* @param side local index of entity
* @param type type of entity
* @param data data on entity
* @return error code
*/
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
// Get base function
auto t_n_hdiv = data.getFTensor1N<3>();
// get normal of face
auto t_normal = getFTensor1Normal();
// Integration weight
double flux_on_entity = 0;
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int rr = 0; rr != nb_dofs; rr++) {
flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
++t_n_hdiv;
++t_data;
}
++t_w;
}
CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
}
};
/**
* Operator used to post-process results for unsaturated infiltration problem.
* Operator should with element for post-processing results, i.e.
* PostProcVolumeOnRefinedMesh
*/
struct OpPostProcMaterial
UnsaturatedFlowElement &cTx;
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
OpPostProcMaterial(UnsaturatedFlowElement &ctx,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
const std::string field_name)
field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Set bloc Id
Tag th_id;
int def_block_id = -1;
CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_block_id);
// Create mesh tag. Tags are created on post-processing mesh and
// visable in post-processor, e.g. Paraview
double zero = 0;
Tag th_theta;
CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_se;
CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
// Tag th_ks;
// CHKERR postProcMesh.tag_get_handle(
// "Ks",1,MB_TYPE_DOUBLE,th_ks,
// MB_TAG_CREAT|MB_TAG_SPARSE,&zero
// );
// CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
Tag th_k;
CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_c;
CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
// Get pressure at integration points
auto t_h = getFTensor0FromVec(cTx.valuesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
// Calculate theta (water content) and save it on mesh tags
CHKERR block_data->calTheta();
double theta = block_data->tHeta;
CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
CHKERR block_data->calSe();
// Calculate Se (effective saturation and save it on the mesh tags)
double Se = block_data->Se;
CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
// Calculate K (hydraulic conductivity) and save it on the mesh tags
CHKERR block_data->calK();
double K = block_data->K;
CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
// Calculate water capacity and save it on the mesh tags
CHKERR block_data->calC();
double C = block_data->C;
CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
// Set block iD
CHKERR postProcMesh.tag_set_data(th_id, &mapGaussPts[gg], 1, &block_id);
++t_h;
++t_coords;
}
}
};
/**
* Finite element implementation called by TS monitor. Element calls other
* finite elements to evaluate material properties and save results on the
* mesh.
*
* \note Element overloaded only FEMethod::postProcess methods where other
* elements are called.
*/
struct MonitorPostProc : public MoFEM::FEMethod {
UnsaturatedFlowElement &cTx;
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
const int fRequency;
MonitorPostProc(UnsaturatedFlowElement &ctx,
boost::shared_ptr<PostProcEle> &post_proc,
boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
const int frequency)
: cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
fRequency(frequency) {}
MoFEMErrorCode preProcess() {
}
MoFEMErrorCode operator()() {
}
MoFEMErrorCode postProcess() {
// Get time step
int step;
CHKERR TSGetTimeStepNumber(ts, &step);
if ((step) % fRequency == 0) {
// Post-process results and save in the file
PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
CHKERR DMoFEMLoopFiniteElements(cTx.dM, "MIX", postProc);
CHKERR postProc->writeFile(
string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
}
// Integrate fluxes on faces where pressure head is applied
CHKERR VecZeroEntries(cTx.ghostFlux);
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
// Run finite element to integrate fluxes
CHKERR DMoFEMLoopFiniteElements(cTx.dM, "MIX_BCVALUE", fluxIntegrate);
CHKERR VecAssemblyBegin(cTx.ghostFlux);
CHKERR VecAssemblyEnd(cTx.ghostFlux);
// accumulate errors from processors
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
// scatter errors to all processors
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
double *ghost_flux;
CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
ghost_flux[0]);
CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
}
};
/// \brief add fields
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
const int order) {
// Fields
// CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
// meshset consisting all entities in mesh
EntityHandle root_set = mField.get_moab().get_root_set();
// add entities to field
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, fluxes);
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, values);
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, values + "_t");
// CHKERR mField.add_ents_to_field_by_type(
// dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
// );
}
CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
CHKERR mField.set_field_order(root_set, MBTET, values, order);
CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
// CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
}
/// \brief add finite elements
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
const std::string &values_name) {
// Define element "MIX". Note that this element will work with fluxes_name
// and values_name. This reflect bilinear form for the problem
values_name + "_t");
// CHKERR
// mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
}
// Define element to integrate natural boundary conditions, i.e. set values.
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
}
// Define element to apply essential boundary conditions.
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
}
}
/**
* \brief Build problem
* @param ref_level mesh refinement on which mesh problem you like to built.
* @return error code
*/
MoFEMErrorCode buildProblem(Range zero_flux_ents,
BitRefLevel ref_level = BitRefLevel().set(0)) {
// Build fields
// Build finite elements
// Build adjacencies of degrees of freedom and elements
// create DM instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
// setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
CHKERR DMSetType(dM, "DMMOFEM");
// mesh is portioned, each process keeps only part of problem
CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
// creates problem in DM
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
// discretised problem creates square matrix (that makes some optimizations)
CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
// set DM options from command line
CHKERR DMSetFromOptions(dM);
// add finite elements
CHKERR DMMoFEMAddElement(dM, "MIX");
CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
// constructor data structures
CHKERR DMSetUp(dM);
// remove zero flux dofs
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
"MIX", "FLUXES", zero_flux_ents);
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
CHKERR DMSetSection(dM, section);
// CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
CHKERR PetscSectionDestroy(&section);
}
boost::shared_ptr<ForcesAndSourcesCore>
feFaceBc; ///< Elemnet to calculate essential bc
boost::shared_ptr<ForcesAndSourcesCore>
feFaceRhs; ///< Face element apply natural bc
boost::shared_ptr<ForcesAndSourcesCore>
feVolInitialPc; ///< Calculate inital boundary conditions
boost::shared_ptr<ForcesAndSourcesCore>
feVolRhs; ///< Assemble residual vector
boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
boost::shared_ptr<MethodForForceScaling>
scaleMethodFlux; ///< Method scaling fluxes
boost::shared_ptr<MethodForForceScaling>
scaleMethodValue; ///< Method scaling values
boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
///< postprocess results at time step
boost::shared_ptr<VectorDouble>
headRateAtGaussPts; ///< Vector keeps head rate
// boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
/**
* \brief Set integration rule to volume elements
*
*/
struct VolRule {
int operator()(int, int, int p_data) const { return 3 * p_data; }
};
/**
* \brief Set integration rule to boundary elements
*
*/
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const {
return 2 * p_data;
}
};
std::vector<int> bcVecIds;
VectorDouble bcVecVals, vecValsOnBc;
/**
* \brief Pre-peprocessing
* Set head pressure rate and get inital essential boundary conditions
*/
struct preProcessVol {
UnsaturatedFlowElement &cTx;
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
preProcessVol(
UnsaturatedFlowElement &ctx,
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
MoFEMErrorCode operator()() {
// Update pressure rates
CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
switch (fePtr->ts_ctx) {
case TSMethod::CTX_TSSETIFUNCTION:
CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
if (!cTx.bcIndices.empty()) {
double scale;
CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
cTx.bcIndices.end());
cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
}
CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(),
&*cTx.vecValsOnBc.begin());
VectorDouble::iterator vit = cTx.bcVecVals.begin();
const NumeredDofEntity *dof_ptr;
for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
it != cTx.bcVecIds.end(); it++, vit++) {
if (auto dof_ptr =
fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
.lock())
dof_ptr->getFieldData() = *vit;
}
} else {
cTx.bcVecIds.resize(0);
cTx.bcVecVals.resize(0);
cTx.vecValsOnBc.resize(0);
}
break;
default:
// don nothing
break;
}
}
};
/**
* \brief Post proces method for volume element
* Assemble vectors and matrices and apply essential boundary conditions
*/
struct postProcessVol {
UnsaturatedFlowElement &cTx;
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
postProcessVol(
UnsaturatedFlowElement &ctx,
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
MoFEMErrorCode operator()() {
switch (fePtr->ts_ctx) {
case TSMethod::CTX_TSSETIJACOBIAN: {
CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
// MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
// std::string wait;
// std::cin >> wait;
CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(), 1, PETSC_NULL,
PETSC_NULL);
CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
// MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
// std::string wait;
// std::cin >> wait;
} break;
case TSMethod::CTX_TSSETIFUNCTION: {
CHKERR VecAssemblyBegin(fePtr->ts_F);
CHKERR VecAssemblyEnd(fePtr->ts_F);
if (!cTx.bcVecIds.empty()) {
CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
INSERT_VALUES);
}
CHKERR VecAssemblyBegin(fePtr->ts_F);
CHKERR VecAssemblyEnd(fePtr->ts_F);
} break;
default:
// don nothing
break;
}
}
};
/**
* \brief Create finite element instances
* @param vol_rule integration rule for volume element
* @param face_rule integration rule for boundary element
* @return error code
*/
MoFEMErrorCode
setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
// create finite element instances
feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
// set integration rule to elements
feFaceBc->getRuleHook = face_rule;
feFaceRhs->getRuleHook = face_rule;
feVolInitialPc->getRuleHook = vol_rule;
feVolLhs->getRuleHook = vol_rule;
feVolRhs->getRuleHook = vol_rule;
// set function hook for finite element postprocessing stage
feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
// Set Piola transform
CHKERR AddHOOps<2, 3, 3>::add(feFaceBc->getOpPtrVector(), {HDIV});
CHKERR AddHOOps<2, 3, 3>::add(feFaceRhs->getOpPtrVector(), {HDIV});
// create method for setting history for fluxes on boundary
scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
new TimeForceScale("-flux_history", false));
// create method for setting history for presure heads on boundary
scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
new TimeForceScale("-head_history", false));
// Set operator to calculate essential boundary conditions
feFaceBc->getOpPtrVector().push_back(
new OpEvaluateBcOnFluxes(*this, "FLUXES"));
// Set operator to calculate initial capillary pressure
CHKERR AddHOOps<3, 3, 3>::add(feVolInitialPc->getOpPtrVector(), {HDIV, L2});
feVolInitialPc->getOpPtrVector().push_back(
new OpEvaluateInitiallHead(*this, "VALUES"));
// set residual face from Neumann terms, i.e. applied pressure
feFaceRhs->getOpPtrVector().push_back(
new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
// set residual finite element operators
headRateAtGaussPts = boost::make_shared<VectorDouble>();
// resAtGaussPts = boost::make_shared<VectorDouble>();
CHKERR AddHOOps<3, 3, 3>::add(feVolRhs->getOpPtrVector(), {HDIV, L2});
feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
string("VALUES") + "_t", headRateAtGaussPts, MBTET));
feVolRhs->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
feVolRhs->getOpPtrVector().push_back(
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
feVolRhs->getOpPtrVector().back().opType =
ForcesAndSourcesCore::UserDataOperator::OPROW;
// set tangent matrix finite element operators
CHKERR AddHOOps<3, 3, 3>::add(feVolLhs->getOpPtrVector(), {HDIV, L2});
feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
string("VALUES") + "_t", headRateAtGaussPts, MBTET));
// feVolLhs->getOpPtrVector().push_back(
// new
// OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
// );
feVolLhs->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
feVolLhs->getOpPtrVector().push_back(
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
feVolLhs->getOpPtrVector().push_back(
new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
// Adding finite elements to DM, time solver will ask for it to assemble
// tangent matrices and residuals
boost::shared_ptr<FEMethod> null;
CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
// setting up post-processing
auto get_post_process_ele = [&]() {
auto post_process = boost::make_shared<PostProcEle>(mField);
CHKERR AddHOOps<3, 3, 3>::add(post_process->getOpPtrVector(), {HDIV, L2});
auto values_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto flux_ptr = boost::make_shared<MatrixDouble>();
post_process->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("VALUES", values_ptr));
post_process->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
post_process->getOpPtrVector().push_back(
new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
using OpPPMap = OpPostProcMapInMoab<3, 3>;
post_process->getOpPtrVector().push_back(
new OpPPMap(post_process->getPostProcMesh(),
post_process->getMapGaussPts(),
{{"VALUES", values_ptr}},
{{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
{}, {}
));
return post_process;
};
auto post_process = get_post_process_ele();
post_process->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
post_process->getOpPtrVector().push_back(
new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
post_process->getMapGaussPts(), "VALUES"));
// Integrate fluxes on boundary
boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
CHKERR AddHOOps<2, 3, 3>::add(flux_integrate->getOpPtrVector(), {HDIV});
flux_integrate->getOpPtrVector().push_back(
new OpIntegrateFluxes(*this, "FLUXES"));
int frequency = 1;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
CHKERR PetscOptionsInt(
"-how_often_output",
"frequency how often results are dumped on hard disk", "", frequency,
&frequency, NULL);
ierr = PetscOptionsEnd();
tsMonitor = boost::shared_ptr<FEMethod>(
new MonitorPostProc(*this, post_process, flux_integrate, frequency));
TsCtx *ts_ctx;
DMMoFEMGetTsCtx(dM, &ts_ctx);
}
Vec D1; ///< Vector with inital head capillary pressure
Vec ghostFlux; ///< Ghost Vector of integrated fluxes
/**
* \brief Create vectors and matrices
* @return Error code
*/
MoFEMErrorCode createMatrices() {
CHKERR DMCreateMatrix(dM, &Aij);
CHKERR DMCreateGlobalVector(dM, &D0);
CHKERR VecDuplicate(D0, &D1);
CHKERR VecDuplicate(D0, &D);
CHKERR VecDuplicate(D0, &F);
int ghosts[] = {0};
int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
&ghostFlux);
}
/**
* \brief Delete matrices and vector when no longer needed
* @return error code
*/
MoFEMErrorCode destroyMatrices() {
CHKERR MatDestroy(&Aij);
CHKERR VecDestroy(&D);
CHKERR VecDestroy(&D0);
CHKERR VecDestroy(&D1);
CHKERR VecDestroy(&F);
CHKERR VecDestroy(&ghostFlux);
}
/**
* \brief Calculate boundary conditions for fluxes
* @return Error code
*/
MoFEMErrorCode calculateEssentialBc() {
// clear vectors
CHKERR VecZeroEntries(D0);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
// clear essential bc indices, it could have dofs from other mesh refinement
bcIndices.clear();
// set operator to calculate essential boundary conditions
CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
CHKERR VecAssemblyBegin(D0);
CHKERR VecAssemblyEnd(D0);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
double norm2D0;
CHKERR VecNorm(D0, NORM_2, &norm2D0);
// CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
}
/**
* \brief Calculate inital pressure head distribution
* @return Error code
*/
MoFEMErrorCode calculateInitialPc() {
// clear vectors
CHKERR VecZeroEntries(D1);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
// Calculate initial pressure head on each element
CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
// Assemble vector
CHKERR VecAssemblyBegin(D1);
CHKERR VecAssemblyEnd(D1);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
// Calculate norm
double norm2D1;
CHKERR VecNorm(D1, NORM_2, &norm2D1);
// CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
}
/**
* \brief solve problem
* @return error code
*/
MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
if (set_initial_pc) {
// Set initial head
CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
}
// Initiate vector from data on the mesh
CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
// Create time solver
TS ts;
CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
// Use backward Euler method
CHKERR TSSetType(ts, TSBEULER);
// Set final time
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
// Setup solver from commabd line
CHKERR TSSetFromOptions(ts);
// Set DM to TS
CHKERR TSSetDM(ts, dM);
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
// Set-up monitor
TsCtx *ts_ctx;
DMMoFEMGetTsCtx(dM, &ts_ctx);
CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
// This add SNES monitor, to show error by fields. It is dirty trick
// to add monitor, so code is hidden from doxygen
CHKERR TSSetSolution(ts, D);
CHKERR TSSetUp(ts);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
#if PETSC_VERSION_GE(3, 7, 0)
{
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
CHKERR SNESMonitorSet(
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
}
#else
{
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
0, 0);
}
#endif
CHKERR TSSolve(ts, D);
// Get statistic form TS and print it
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);
PetscPrintf(PETSC_COMM_WORLD,
"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
"%D, linits %D\n",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
};
} // namespace MixTransport
#endif // __UNSATURATD_FLOW_HPP__
static Index< 'K', 3 > K
constexpr double a
static PetscErrorCode ierr
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
@ ROW
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
@ F
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 build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition: plastic.cpp:163
constexpr auto field_name
Set integration rule to boundary elements.
virtual MoFEMErrorCode calC()
virtual MoFEMErrorCode calSe()
virtual void printMatParameters(const int id, const std::string &prefix) const =0
virtual MoFEMErrorCode calTheta()
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static double ePsilon0
Regularization parameter.
virtual MoFEMErrorCode calK()
virtual MoFEMErrorCode calDiffC()
double sCale
Scale time dependent eq.
double K
Hydraulic conductivity [L/s].
Range tEts
Elements with this material.
double C
Capacity [S^2/L^2].
double Se
Effective saturation.
static double scaleZ
Scale z direction.
virtual double initialPcEval() const =0
Initialize head.
double h_t
rate of hydraulic head
static double ePsilon1
Regularization parameter.
virtual MoFEMErrorCode calDiffK()
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
VectorDouble valuesAtGaussPts
values at integration points on element
VectorDouble divergenceAtGaussPts
divergence at integration points on element
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
boost::function< double(const double x, const double y, const double z)> hookFun
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
MoFEMErrorCode preProcess()
function is run at the beginning of loop
UnsaturatedFlowElement & cTx
MoFEMErrorCode postProcess()
function is run at the end of loop
const int fRequency
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< MethodForForceScaling > valueScale
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
int operator()(int, int, int p_data) const
boost::shared_ptr< ForcesAndSourcesCore > fePtr
boost::shared_ptr< ForcesAndSourcesCore > fePtr
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
Vec D1
Vector with inital head capillary pressure.
map< int, boost::shared_ptr< BcData > > BcMap
MaterialsDoubleMap dMatMap
materials database
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
Vec ghostFlux
Ghost Vector of integrated fluxes.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
boost::shared_ptr< FEMethod > tsMonitor
BcMap bcValueMap
Store boundary condition for head capillary pressure.
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
DM dM
Discrete manager for unsaturated flow problem.
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
TS ts
time solver
PetscReal ts_t
time
Vec & ts_F
residual vector
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Mat & ts_B
Preconditioner for ts_A.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:154
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Force scale operator for reading two columns.
Set integration rule.