v0.9.0
UnsaturatedFlow.hpp
/** \file UnsaturatedFlow.hpp
* \brief Mix implementation of transport element
* \example UnsaturatedFlow.hpp
*
* \ingroup mofem_mix_transport_elem
*
*/
/*
* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
namespace MixTransport {
/**
* \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
double sCale; ///< Scale time dependent eq.
// double Ks; ///< Saturated hydraulic conductivity [m/day]
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 initalPcEval() 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");
}
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");
}
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
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
: 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
*/
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 {
Range eNts;
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
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,
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
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
flux_name, UserDataOperator::OPROW),
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
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
// 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();
// Get material parameters
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get divergence
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 z = t_coords(2); /// z-coordinate at Gauss pt
// Calculate pressure gradient
noalias(nF) -= alpha * (t_h - z) * 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
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
val_name, UserDataOperator::OPROW),
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
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
// Get pressure rate
// Flux divergence
// 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
const std::string flux_name)
flux_name, flux_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
/**
* \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,
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
// 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
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
value_name, value_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
/**
* \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,
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
// // Get pressure
// auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
// Get pressure rate
// 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
/**
* \brief Constructor
*/
const std::string &val_name_row,
const std::string &flux_name_col)
val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
/**
* \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,
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();
// 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;
CHKERR getDivergenceOfHDivBaseFunctions(col_side, col_type, col_data,
gg, divVec);
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
/**
* \brief Constructor
*/
const std::string &flux_name_col,
const std::string &val_name_row)
flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
/**
* \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,
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
// 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 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;
CHKERR getDivergenceOfHDivBaseFunctions(row_side, row_type, row_data,
gg, divVec);
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
const std::string &val_name)
val_name, UserDataOperator::OPROW),
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (data.getFieldData().size() == 0)
int nb_dofs = data.getFieldData().size();
int nb_gauss_pts = data.getN().size1();
if (nb_dofs != static_cast<int>(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
for (int 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) * getVolume();
nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
nF += alpha * block_data->initalPcEval() * 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
/**
* \brief Constructor
*/
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,
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
auto t_w = getFTensor0IntegrationWeight();
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
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
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,
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);
CHKERR postProcMesh.tag_set_data(th_id, &fe_ent, 1, &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
// 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);
++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 methos where other
* elements are called.
*/
struct MonitorPostProc : public MoFEM::FEMethod {
boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
const int fRequency;
boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc,
boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
const int frequency)
: cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
fRequency(frequency) {}
}
}
// 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 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 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
*/
// 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
// creates problem in DM
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
// discretised problem creates square matrix (that makes some optimizations)
// set DM options from command line
CHKERR DMSetFromOptions(dM);
// add finite elements
CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
// constructor data structures
CHKERR DMSetUp(dM);
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
CHKERR DMSetDefaultSection(dM, section);
CHKERR DMSetDefaultGlobalSection(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 2 * p_data + 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;
/**
* \brief Pre-peprocessing
* Set head pressute rate and get inital essential boundary conditions
*/
struct preProcessVol {
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
// 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:
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());
cTx.bcVecVals *= scale;
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++) {
CHKERR fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it,
&dof_ptr);
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 {
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
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()) {
// cerr << mArk << endl;
// cerr << "a " << cTx.vecValsOnBc << endl;
// cerr << "a " << cTx.bcVecVals << endl;
&*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
INSERT_VALUES);
}
CHKERR VecAssemblyBegin(fePtr->ts_F);
CHKERR VecAssemblyEnd(fePtr->ts_F);
// CHKERR VecView(fePtr->ts_F,PETSC_VIEWER_STDOUT_WORLD);
// CHKERR
// fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
// fePtr->problemPtr,"VALUES",string("FLUXES")+"_residual",
// ROW,fePtr->ts_F,INSERT_VALUES,SCATTER_REVERSE
// );
} 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
*/
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>(
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
// 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);
// 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
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>();
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
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
boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
CHKERR post_process->generateReferenceElementMesh();
CHKERR post_process->addFieldValuesPostProc("VALUES");
CHKERR post_process->addFieldValuesPostProc("VALUES_t");
// CHKERR post_process->addFieldValuesPostProc("FLUXES_residual");
CHKERR post_process->addFieldValuesPostProc("FLUXES");
post_process->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
post_process->getOpPtrVector().push_back(
new OpPostProcMaterial(*this, post_process->postProcMesh,
post_process->mapGaussPts, "VALUES"));
// Integrate fluxes on boundary
boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
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);
ts_ctx->get_postProcess_to_do_Monitor().push_back(tsMonitor);
}
Vec D1; ///< Vector with inital head capillary pressure
Vec ghostFlux; ///< Ghost Vector of integrated fluxes
/**
* \brief Create vectors and matrices
* @return Error code
*/
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,
}
/**
* \brief Delete matrices and vector when no longer needed
* @return error code
*/
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
*/
// 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 VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(D0);
CHKERR VecAssemblyEnd(D0);
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
*/
// 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
// Assemble vector
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(D1);
CHKERR VecAssemblyEnd(D1);
// 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__