#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
struct GenericMaterial {
const std::string &prefix) const = 0;
"Not implemented how to calculate hydraulic conductivity");
}
"Not implemented how to calculate derivative of hydraulic "
"conductivity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
};
struct UnsaturatedFlowElement : public MixTransportElement {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
}
int &block_id) const {
for (MaterialsDoubleMap::const_iterator mit =
dMatMap.begin();
if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
block_id = mit->first;
}
}
"Element not found, no material data");
}
struct BcData {
boost::function<
double(
const double x,
const double y,
const double z)>
};
typedef map<int, boost::shared_ptr<BcData>>
BcMap;
const double x, const double y, const double z,
double &value) {
int block_id = -1;
} else {
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
value =
bcValueMap.at(block_id)->hookFun(x, y, z);
} else {
}
} else {
value = 0;
}
}
const double y, const double z, double &flux) {
int block_id = -1;
} else {
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
flux =
bcFluxMap.at(block_id)->hookFun(x, y, z);
} else {
}
} else {
flux = 0;
}
}
struct OpRhsBcOnValues
boost::shared_ptr<MethodForForceScaling>
valueScale;
boost::shared_ptr<MethodForForceScaling> &value_scale)
if (data.getFieldData().size() == 0)
nF.resize(data.getIndices().size());
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double x, y, z;
double value;
const double beta =
w * (value - z);
noalias(
nF) += beta * prod(data.getVectorN<3>(gg),
getNormal());
}
}
}
};
struct OpResidualFlux
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv = data.getFTensor1N<3>();
auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
divVec.resize(nb_dofs,
false);
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
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);
const double K = block_data->K;
const double scaleZ = block_data->scaleZ;
const double z = t_coords(2);
noalias(
nF) -= alpha * (t_h - z * scaleZ) *
divVec;
for (int rr = 0; rr != nb_dofs; rr++) {
t_nf += alpha * (1 /
K) * (t_n_hdiv(
i) * t_flux(
i));
++t_n_hdiv;
++t_nf;
}
++t_h;
++t_flux;
++t_coords;
++t_w;
}
&*data.getIndices().begin(), &*
nF.begin(),
ADD_VALUES);
}
};
struct OpResidualMass
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
const double scale = block_data->sCale;
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);
const double C = block_data->C;
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;
}
ADD_VALUES);
}
};
struct OpTauDotSigma_HdivHdiv
const std::string flux_name)
}
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);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
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);
const double K = block_data->K;
const double alpha = t_w * vol;
const double beta = alpha * (1 /
K);
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;
}
&*col_data.getIndices().begin(), &*
nN.data().begin(),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
&*row_data.getIndices().begin(),
&*
transNN.data().begin(), ADD_VALUES);
}
}
};
struct OpVU_L2L2
}
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);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
auto t_n_row = row_data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double alpha = t_w * vol *
scale;
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->calDiffC();
const double C = block_data->C;
const double diffC = block_data->diffC;
for (int kk = 0; kk != nb_row; kk++) {
auto t_n_col = col_data.getFTensor0N(gg, 0);
for (int ll = 0; ll != nb_col; ll++) {
t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
++t_n_col;
++t_a;
}
++t_n_row;
}
++t_w;
++t_coords;
++t_h;
++t_h_t;
}
&col_data.getIndices()[0], &*
nN.data().begin(),
ADD_VALUES);
}
};
struct OpVDivSigma_L2Hdiv
const std::string &val_name_row,
const std::string &flux_name_col)
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)
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
nN.resize(nb_row, nb_col,
false);
auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
v = t_base_diff_hdiv(
i,
i);
++t_base_diff_hdiv;
}
noalias(
nN) += alpha * outer_prod(row_data.getN(gg),
divVec);
}
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &
nN(0, 0), ADD_VALUES);
}
};
struct OpDivTauU_HdivL2
const std::string &flux_name_col,
const std::string &val_name_row)
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);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
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->calDiffK();
const double K = block_data->K;
const double diffK = block_data->diffK;
double alpha = t_w * vol;
v = t_base_diff_hdiv(
i,
i);
++t_base_diff_hdiv;
}
noalias(
nN) -= alpha * outer_prod(
divVec, col_data.getN(gg));
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;
}
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &
nN(0, 0), ADD_VALUES);
}
};
struct OpEvaluateInitiallHead
const std::string &val_name)
if (data.getFieldData().size() == 0)
auto nb_dofs = data.getFieldData().size();
if (nb_dofs != data.getN().size2()) {
"wrong number of dofs");
}
nN.resize(nb_dofs, nb_dofs,
false);
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
for (auto gg = 0; gg != nb_gauss_pts; gg++) {
nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
nF += alpha * block_data->initialPcEval() * data.getN(gg);
}
&*
nF.begin(), INSERT_VALUES);
}
};
struct OpIntegrateFluxes
const std::string fluxes_name)
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
auto t_n_hdiv = data.getFTensor1N<3>();
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;
}
}
};
struct OpPostProcMaterial
std::vector<EntityHandle> &map_gauss_pts,
field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
Tag th_id;
int def_block_id = -1;
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_block_id);
double zero = 0;
Tag th_theta;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_se;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_k;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_c;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
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);
CHKERR block_data->calTheta();
double theta = block_data->tHeta;
double Se = block_data->Se;
double K = block_data->K;
double C = block_data->C;
++t_h;
++t_coords;
}
}
};
boost::shared_ptr<PostProcEle>
postProc;
boost::shared_ptr<PostProcEle> &post_proc,
boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
const int frequency)
}
}
int step;
PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
}
double *ghost_flux;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Flux at time %6.4g %6.4g\n",
ts_t,
ghost_flux[0]);
}
};
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
MBTET, fluxes);
MBTET, values);
MBTET, values + "_t");
}
}
const std::string &values_name) {
values_name + "_t");
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
dMatMap[it->getMeshsetId()]->tEts, MBTET,
"MIX");
}
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
bcValueMap[it->getMeshsetId()]->eNts, MBTRI,
"MIX_BCVALUE");
}
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
bcFluxMap[it->getMeshsetId()]->eNts, MBTRI,
"MIX_BCFLUX");
}
}
"MIX", "FLUXES", zero_flux_ents);
PetscSection section;
CHKERR PetscSectionDestroy(§ion);
}
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
feVolLhs;
boost::shared_ptr<MethodForForceScaling>
boost::shared_ptr<MethodForForceScaling>
boost::shared_ptr<VectorDouble>
int operator()(
int,
int,
int p_data)
const {
return 3 * p_data; }
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data;
}
};
struct preProcessVol {
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr
)
CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
fePtr->problemPtr,
"VALUES",
string(
"VALUES") +
"_t",
ROW,
fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
case TSMethod::CTX_TSSETIFUNCTION:
CHKERR VecSetOption(
fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
}
const NumeredDofEntity *dof_ptr;
for (std::vector<int>::iterator it =
cTx.
bcVecIds.begin();
if (auto dof_ptr =
fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
.lock())
dof_ptr->getFieldData() = *vit;
}
} else {
}
break;
default:
break;
}
}
};
struct postProcessVol {
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr
)
case TSMethod::CTX_TSSETIJACOBIAN: {
CHKERR MatAssemblyBegin(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
PETSC_NULL);
CHKERR MatAssemblyBegin(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
} break;
case TSMethod::CTX_TSSETIFUNCTION: {
INSERT_VALUES);
}
} break;
default:
break;
}
}
};
ForcesAndSourcesCore::RuleHookFun face_rule =
FaceRule()) {
feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
new OpEvaluateBcOnFluxes(*this, "FLUXES"));
new OpEvaluateInitiallHead(*this, "VALUES"));
CHKERR AddHOOps<3, 3, 3>::add(
feVolRhs->getOpPtrVector(), {HDIV, L2});
feVolRhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(
new OpValuesAtGaussPts(*this, "VALUES"));
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;
CHKERR AddHOOps<3, 3, 3>::add(
feVolLhs->getOpPtrVector(), {HDIV, L2});
feVolLhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(
new OpValuesAtGaussPts(*this, "VALUES"));
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpVU_L2L2(*
this,
"VALUES"));
new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
boost::shared_ptr<FEMethod> null;
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"));
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");
"-how_often_output",
"frequency how often results are dumped on hard disk", "", frequency,
&frequency, NULL);
ierr = PetscOptionsEnd();
}
CHKERR DMCreateMatrix(dM, &Aij);
CHKERR DMCreateGlobalVector(dM, &D0);
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);
}
CHKERR VecDestroy(&ghostFlux);
}
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
bcIndices.clear();
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
double norm2D0;
CHKERR VecNorm(D0, NORM_2, &norm2D0);
PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
}
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
double norm2D1;
CHKERR VecNorm(D1, NORM_2, &norm2D1);
PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
}
if (set_initial_pc) {
}
TS ts;
CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
CHKERR TSSetType(ts, TSBEULER);
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
SNES snes;
#if PETSC_VERSION_GE(3, 7, 0)
{
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
void *))SNESMonitorFields,
}
#else
{
void *))SNESMonitorFields,
0, 0);
}
#endif
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);
}
};
}
#endif // __UNSATURATD_FLOW_HPP__