Testing broken spaces. Implementations works for 2d and 3d meshes, is aiming to test H-div broken base functions, and L2 base on skeleton.
Also, it test block matrix with Schur complement.
static char help[] =
"...\n\n";
constexpr
bool debug =
false;
static boost::shared_ptr<SetUpSchur>
protected:
};
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
DMType dm_name_mg = "DMMOFEM_MG";
auto core_log = logging::core::get();
core_log->add_sink(
core_log->add_sink(
simple->getAddBoundaryFE() =
true;
auto add_shared_entities_on_skeleton = [&]() {
auto boundary_meshset =
simple->getBoundaryMeshSet();
auto skeleton_meshset =
simple->getSkeletonMeshSet();
bdy_ents, true);
0,
simple->getDim() - 1, skeleton_ents,
true);
skeleton_ents = subtract(skeleton_ents, bdy_ents);
CHKERR m_field.
get_moab().add_entities(skeleton_meshset, skeleton_ents);
};
CHKERR add_shared_entities_on_skeleton();
enum bases {
AINSWORTH,
AINSWORTH_LOBATTO,
DEMKOWICZ,
BERNSTEIN,
LASBASETOP
};
const char *list_bases[] = {"ainsworth", "ainsworth_lobatto", "demkowicz",
"bernstein"};
PetscBool flg;
PetscInt choice_base_value = AINSWORTH;
LASBASETOP, &choice_base_value, &flg);
if (flg != PETSC_TRUE)
if (choice_base_value == AINSWORTH)
if (choice_base_value == AINSWORTH_LOBATTO)
else if (choice_base_value == DEMKOWICZ)
else if (choice_base_value == BERNSTEIN)
enum spaces { hdiv, hcurl, last_space };
const char *list_spaces[] = {"hdiv", "hcurl"};
PetscInt choice_space_value = hdiv;
last_space, &choice_space_value, &flg);
if (flg != PETSC_TRUE)
if (choice_space_value == hdiv)
else if (choice_space_value == hcurl)
PETSC_NULL);
CHKERR simple->addDomainBrokenField(
"BROKEN", space, base, 1);
CHKERR bc_mng->removeSideDOFs(
simple->getProblemName(),
"ZERO_FLUX",
auto assemble_domain_lhs = [&](auto &pip) {
IT>::OpMixDivTimesScalar<SPACE_DIM>;
return 1;
};
pip.push_back(
new OpHdivHdiv(
"BROKEN",
"BROKEN", beta));
auto unity = []() constexpr { return 1; };
pip.push_back(
new OpHdivU(
"BROKEN",
"U", unity,
true));
op_loop_skeleton_side->getOpPtrVector(), {});
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getOpPtrVector(), {HDIV});
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
IT>::OpBrokenSpaceConstrain<1>;
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC("HYBRID", broken_data_ptr, 1., true, false));
constexpr int partition = 1;
op_print->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
if (
auto op_ptr =
dynamic_cast<BdyEleOp *
>(base_op_ptr)) {
auto fe_method = op_ptr->getFEMethod();
auto num_fe = fe_method->numeredEntFiniteElementPtr;
if (num_fe->getPStatus() & PSTATUS_SHARED)
MOFEM_LOG(
"SELF", Sev::inform) <<
"Num FE: " << *num_fe;
}
}
};
op_loop_skeleton_side->getOpPtrVector().push_back(op_print);
};
pip.push_back(op_loop_skeleton_side);
};
auto assemble_domain_rhs = [&](auto &pip) {
AT>::LinearForm<IT>::OpSource<1, 1>;
auto source = [&](
const double x,
const double y,
const double z) constexpr {
return -1;
};
};
CHKERR assemble_domain_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR assemble_domain_rhs(pip_mng->getOpDomainRhsPipeline());
TetPolynomialBase::switchCacheBaseOn<HDIV>(
{pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
TetPolynomialBase::switchCacheBaseOn<L2>(
{pip_mng->getDomainLhsFE().get(), pip_mng->getDomainRhsFE().get()});
auto ksp = pip_mng->createKSP();
CHKERR KSPSetFromOptions(ksp);
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp <= Done";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve <= Done";
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
} else {
auto ksp = pip_mng->createKSP();
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSetUp <= Done";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve";
MOFEM_LOG(
"TIMER", Sev::inform) <<
"KSPSolve <= Done";
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
}
auto check_residual = [&](
auto x,
auto f) {
auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
domain_rhs.clear();
auto div_flux_ptr = boost::make_shared<VectorDouble>();
"BROKEN", div_flux_ptr));
AT>::LinearForm<IT>::OpBaseTimesScalarField<1>;
domain_rhs.push_back(new OpUDivFlux("U", div_flux_ptr, beta));
auto source = [&](
const double x,
const double y,
const double z) constexpr { return 1; };
AT>::LinearForm<IT>::OpSource<1, 1>;
IT>::OpMixDivTimesU<3, 1, SPACE_DIM>;
AT>::LinearForm<IT>::OpBaseTimesVector<3, 3, 1>;
auto flux_ptr = boost::make_shared<MatrixDouble>();
domain_rhs.push_back(
boost::shared_ptr<VectorDouble> u_ptr =
boost::make_shared<VectorDouble>();
domain_rhs.push_back(
new OpHDivH(
"BROKEN", u_ptr, beta));
domain_rhs.push_back(
new OpHdivFlux(
"BROKEN", flux_ptr, beta));
op_loop_skeleton_side->getOpPtrVector(), {});
auto broken_data_ptr =
boost::make_shared<std::vector<BrokenBaseSideData>>();
op_loop_domain_side->getOpPtrVector(), {HDIV});
op_loop_domain_side->getOpPtrVector().push_back(
auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_domain_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
IT>::OpBrokenSpaceConstrainDHybrid<1>;
IT>::OpBrokenSpaceConstrainDFlux<1>;
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dHybrid("HYBRID", broken_data_ptr, 1.));
auto hybrid_ptr = boost::make_shared<MatrixDouble>();
op_loop_skeleton_side->getOpPtrVector().push_back(
op_loop_skeleton_side->getOpPtrVector().push_back(
new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
domain_rhs.push_back(op_loop_skeleton_side);
CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
pip_mng->getDomainRhsFE()->f =
f;
pip_mng->getSkeletonRhsFE()->f =
f;
pip_mng->getDomainRhsFE()->x = x;
pip_mng->getSkeletonRhsFE()->x = x;
pip_mng->getDomainRhsFE());
CHKERR VecGhostUpdateBegin(
f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
f, ADD_VALUES, SCATTER_REVERSE);
double fnrm;
MOFEM_LOG_C(
"AT", Sev::inform,
"Residual %3.4e", fnrm);
constexpr
double eps = 1e-8;
"Residual norm larger than accepted");
};
auto calculate_error = [&]() {
auto &domain_rhs = pip_mng->getOpDomainRhsPipeline();
domain_rhs.clear();
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
auto flux_val_ptr = boost::make_shared<MatrixDouble>();
auto div_val_ptr = boost::make_shared<VectorDouble>();
auto source_ptr = boost::make_shared<VectorDouble>();
domain_rhs.push_back(
domain_rhs.push_back(
"BROKEN", div_val_ptr));
auto source = [&](
const double x,
const double y,
const double z) constexpr { return -1; };
enum { DIV, GRAD,
LAST };
domain_rhs.push_back(
u_grad_ptr, mpi_vec, GRAD, flux_val_ptr));
pip_mng->getDomainRhsFE());
CHKERR VecAssemblyBegin(mpi_vec);
CHKERR VecAssemblyEnd(mpi_vec);
const double *error_ind;
CHKERR VecGetArrayRead(mpi_vec, &error_ind);
<< "Approximation error ||div flux - source||: "
<< std::sqrt(error_ind[DIV]);
MOFEM_LOG(
"AT", Sev::inform) <<
"Approximation error ||grad-flux||: "
<< std::sqrt(error_ind[GRAD]);
CHKERR VecRestoreArrayRead(mpi_vec, &error_ind);
}
};
auto get_post_proc_fe = [&]() {
auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
boost::make_shared<
post_proc_fe->getOpPtrVector().push_back(op_loop_side);
op_loop_side->getOpPtrVector(), {HDIV});
auto u_vec_ptr = boost::make_shared<VectorDouble>();
auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"U", u_vec_ptr}},
{{"BROKEN", flux_mat_ptr}},
{}, {})
);
return post_proc_fe;
};
auto post_proc_fe = get_post_proc_fe();
simple->getBoundaryFEName(), post_proc_fe);
CHKERR post_proc_fe->writeFile(
"out_result.h5m");
}
}
private:
};
CHKERR KSPSetFromOptions(ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
MOFEM_LOG(
"AT", Sev::inform) <<
"Setup Schur pc";
auto create_sub_dm = [&]() {
auto create_dm = [&](
std::string problem_name,
std::vector<std::string> fe_names,
std::vector<std::string> fields,
auto dm_type
) {
auto create_dm_imp = [&]() {
for (auto fe : fe_names) {
}
for (auto field : fields) {
}
};
create_dm_imp(),
"Error in creating schurDM. It is possible that schurDM is "
"already created");
return dm;
};
auto schur_dm = create_dm(
"SCHUR",
{"HYBRID"},
"DMMOFEM_MG");
auto block_dm = create_dm(
"BLOCK",
{"BROKEN", "U"},
"DMMOFEM");
return std::make_tuple(schur_dm, block_dm);
};
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
{
{
simple->getDomainFEName(),
{
{"BROKEN", "BROKEN"},
{"U", "U"},
{"BROKEN", "U"},
{"U", "BROKEN"}
}
},
{
simple->getSkeletonFEName(),
{
{"BROKEN", "HYBRID"}, {"HYBRID", "BROKEN"}
}
}
}
);
{schur_dm, block_dm}, block_mat_data,
{"BROKEN", "U"}, {nullptr, nullptr}, true
);
};
auto set_ops = [&](auto schur_dm) {
boost::shared_ptr<BlockStructure> block_data;
pip_mng->getOpDomainLhsPipeline().push_front(
pip_mng->getOpDomainLhsPipeline().push_back(
S, true, true)
);
auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble Begin";
};
post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
post_proc_schur_lhs_ptr]() {
MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble End";
CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble Finish";
};
ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
};
auto set_pc = [&](auto pc, auto block_dm) {
CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
};
auto set_diagonal_pc = [&](auto pc, auto schur_dm) {
}
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
return pc_raw;
};
auto set_pc_p_mg = [&](auto dm, auto pc) {
PetscBool same = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
if (same) {
}
};
CHKERR set_pc_p_mg(schur_dm, get_pc(subksp[1]));
};
auto [schur_dm, block_dm] = create_sub_dm();
auto nested_mat_data = get_nested_mat_data(schur_dm, block_dm);
}
CHKERR MatSetDM(S, PETSC_NULL);
CHKERR MatSetBlockSize(S, bs);
DM solver_dm;
CHKERR KSPGetDM(ksp, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
CHKERR set_diagonal_pc(pc, schur_dm);
} else {
"PC is not set to PCFIELDSPLIT");
}
}
boost::shared_ptr<SetUpSchur>
}