13 #include <petsc/private/petscimpl.h>
15 using namespace MoFEM;
17 static char help[] =
"...\n\n";
19 #undef PYTHON_INIT_SURFACE
21 #ifdef PYTHON_INIT_SURFACE
22 #include <boost/python.hpp>
23 #include <boost/python/def.hpp>
24 namespace bp = boost::python;
26 struct SurfacePython {
27 SurfacePython() =
default;
28 virtual ~SurfacePython() =
default;
35 auto main_module = bp::import(
"__main__");
36 mainNamespace = main_module.attr(
"__dict__");
37 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
39 surfaceFun = mainNamespace[
"surface"];
41 }
catch (bp::error_already_set
const &) {
51 double x,
double y,
double z,
double eta,
double &s
58 s = bp::extract<double>(surfaceFun(x, y, z,
eta));
60 }
catch (bp::error_already_set
const &) {
69 bp::object mainNamespace;
70 bp::object surfaceFun;
73 static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
125 template <CoordinateTypes COORD_TYPE>
178 double tol = std::numeric_limits<float>::epsilon();
197 constexpr
int sub_stepping = 16;
198 return std::min(1.,
static_cast<double>(ts_step) / sub_stepping);
203 auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
204 auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
207 if (
h >= -1 &&
h < 1)
214 return diff *
h + ave;
223 auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
224 auto get_f_dh = [](
const double h) {
return 4 *
W * (3 *
h *
h - 1); };
230 return md * (1 -
h *
h);
236 const double h2 =
h *
h;
237 const double h3 = h2 *
h;
239 return md * (2 * h3 - 3 * h2 + 1);
241 return md * (-2 * h3 - 3 * h2 + 1);
246 return md * (6 *
h * (
h - 1));
248 return md * (-6 *
h * (
h + 1));
264 constexpr
double R = 0.0125;
265 constexpr
double A =
R * 0.2;
266 const double theta = atan2(
r, y);
267 const double w =
R +
A * cos(
n * theta);
268 const double d = std::sqrt(
r *
r + y * y);
269 return tanh((
w -
d) / (
eta * std::sqrt(2)));
273 constexpr
double y0 = 0.5;
274 constexpr
double R = 0.4;
276 const double d = std::sqrt(
r *
r + y * y);
277 return tanh((
R -
d) / (
eta * std::sqrt(2)));
281 constexpr
double water_height = 0.;
282 return tanh((water_height - y) / (
eta * std::sqrt(2)));
287 return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
297 auto init_h = [](
double r,
double y,
double theta) {
298 #ifdef PYTHON_INIT_SURFACE
300 if (
auto ptr = surfacePythonWeakPtr.lock()) {
302 "error eval python");
317 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
322 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
331 CHKERR moab.add_entities(*out_meshset,
r);
332 CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
333 out_meshset->get_ptr(), 1);
344 std::vector<EntityHandle> ents_vec;
350 auto dofs = prb_ptr->numeredRowDofsPtr;
353 ents_vec.reserve(std::distance(lo_it, hi_it));
355 for (; lo_it != hi_it; ++lo_it) {
356 ents_vec.push_back((*lo_it)->getEnt());
359 std::sort(ents_vec.begin(), ents_vec.end());
360 auto it = std::unique(ents_vec.begin(), ents_vec.end());
362 r.insert_list(ents_vec.begin(), it);
372 std::vector<EntityHandle> ents_vec;
377 auto dofs = prb_ptr->numeredRowDofsPtr;
378 ents_vec.reserve(dofs->size());
380 for (
auto d : *dofs) {
381 ents_vec.push_back(
d->getEnt());
384 std::sort(ents_vec.begin(), ents_vec.end());
385 auto it = std::unique(ents_vec.begin(), ents_vec.end());
387 r.insert_list(ents_vec.begin(), it);
452 PetscReal
a, Mat
A, Mat B,
463 boost::shared_ptr<SnesCtx>
465 boost::shared_ptr<TsCtx>
492 std::vector<Range> findEntitiesCrossedByPhaseInterface(
size_t overlap);
514 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
517 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
528 CHKERR boundaryCondition();
538 MOFEM_LOG(
"FS", Sev::inform) <<
"Read mesh for problem";
541 simple->getParentAdjacencies() =
true;
561 const char *coord_type_names[] = {
"cartesian",
"polar",
"cylindrical",
566 MOFEM_LOG(
"FS", Sev::inform) <<
"Approximation order = " <<
order;
568 <<
"Number of refinement levels nb_levels = " <<
nb_levels;
575 "-rho_m", &
rho_m, PETSC_NULL);
585 "Height of the well in energy functional",
"-W",
599 MOFEM_LOG(
"FS", Sev::inform) <<
"Acceleration a0 = " <<
a0;
600 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase density rho_m = " <<
rho_m;
601 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Minus\" phase viscosity mu_m = " <<
mu_m;
602 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase density rho_p = " <<
rho_p;
603 MOFEM_LOG(
"FS", Sev::inform) <<
"\"Plus\" phase viscosity mu_p = " <<
mu_p;
606 <<
"Height of the well in energy functional W = " <<
W;
607 MOFEM_LOG(
"FS", Sev::inform) <<
"Characteristic mesh size h = " <<
h;
608 MOFEM_LOG(
"FS", Sev::inform) <<
"Mobility md = " <<
md;
612 MOFEM_LOG(
"FS", Sev::inform) <<
"Average viscosity mu_ave = " <<
mu_ave;
613 MOFEM_LOG(
"FS", Sev::inform) <<
"Difference viscosity mu_diff = " <<
mu_diff;
643 auto set_problem_bit = [&]() {
666 #ifdef PYTHON_INIT_SURFACE
667 auto get_py_surface_init = []() {
668 auto py_surf_init = boost::make_shared<SurfacePython>();
669 CHKERR py_surf_init->surfaceInit(
"surface.py");
670 surfacePythonWeakPtr = py_surf_init;
673 auto py_surf_init = get_py_surface_init();
680 auto dm =
simple->getDM();
684 auto reset_bits = [&]() {
688 start_mask[s] =
true;
690 CHKERR bit_mng->lambdaBitRefLevel(
699 auto add_parent_field = [&](
auto fe,
auto op,
auto field) {
700 return setParentDofs(
704 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
712 auto h_ptr = boost::make_shared<VectorDouble>();
713 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
714 auto g_ptr = boost::make_shared<VectorDouble>();
715 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
717 auto set_generic = [&](
auto fe) {
719 auto &pip = fe->getOpPtrVector();
727 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
730 fe_parent->getOpPtrVector(), {H1});
736 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
741 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
749 auto post_proc = [&](
auto exe_test) {
751 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
752 post_proc_fe->exeTestHook = exe_test;
754 CHKERR set_generic(post_proc_fe);
758 post_proc_fe->getOpPtrVector().push_back(
760 new OpPPMap(post_proc_fe->getPostProcMesh(),
761 post_proc_fe->getMapGaussPts(),
763 {{
"H", h_ptr}, {
"G", g_ptr}},
765 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
776 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
781 auto solve_init = [&](
auto exe_test) {
784 pip_mng->getOpDomainRhsPipeline().clear();
785 pip_mng->getOpDomainLhsPipeline().clear();
787 auto set_domain_rhs = [&](
auto fe) {
790 auto &pip = fe->getOpPtrVector();
792 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
793 pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
795 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
796 pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
800 auto set_domain_lhs = [&](
auto fe) {
804 auto &pip = fe->getOpPtrVector();
806 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
807 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
810 CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
813 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
816 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
821 auto create_subdm = [&]() {
822 auto level_ents_ptr = boost::make_shared<Range>();
828 CHKERR DMSetType(subdm,
"DMMOFEM");
834 for (
auto f : {
"H",
"G"}) {
840 if constexpr (
debug) {
844 dm_ents.subset_by_type(MBVERTEX));
845 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
853 auto subdm = create_subdm();
855 pip_mng->getDomainRhsFE().reset();
856 pip_mng->getDomainLhsFE().reset();
859 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
860 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
862 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
863 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
866 auto snes = pip_mng->createSNES(subdm);
869 auto set_section_monitor = [&](
auto solver) {
871 CHKERR SNESMonitorSet(solver,
874 (
void *)(snes_ctx_ptr.get()),
nullptr);
879 auto print_section_field = [&]() {
885 CHKERR PetscSectionGetNumFields(section, &num_fields);
886 for (
int f = 0;
f < num_fields; ++
f) {
890 <<
"Field " <<
f <<
" " << std::string(
field_name);
896 CHKERR set_section_monitor(snes);
897 CHKERR print_section_field();
899 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
903 CHKERR SNESSetFromOptions(snes);
904 CHKERR SNESSolve(snes, PETSC_NULL,
D);
906 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
907 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
918 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
929 pip_mng->getOpDomainRhsPipeline().clear();
930 pip_mng->getOpDomainLhsPipeline().clear();
933 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
935 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
937 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
939 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
941 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
943 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
945 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
967 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
968 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
969 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
970 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
972 pip_mng->getDomainLhsFE().reset();
973 pip_mng->getDomainRhsFE().reset();
974 pip_mng->getBoundaryLhsFE().reset();
975 pip_mng->getBoundaryRhsFE().reset();
980 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
981 return setParentDofs(
985 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
994 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
995 return setParentDofs(
999 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1009 auto create_run_parent_op = [&](
auto parent_fe_ptr,
auto this_fe_ptr,
1011 auto parent_mask = fe_bit;
1019 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr,
auto fe_bit) {
1020 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
1022 parents_elems_ptr_vec.emplace_back(
1023 boost::make_shared<DomianParentEle>(mField));
1025 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1026 create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
1028 return parents_elems_ptr_vec[0];
1032 auto solve_projection = [&](
auto exe_test) {
1035 auto set_domain_rhs = [&](
auto fe) {
1037 auto &pip = fe->getOpPtrVector();
1039 auto u_ptr = boost::make_shared<MatrixDouble>();
1040 auto p_ptr = boost::make_shared<VectorDouble>();
1041 auto h_ptr = boost::make_shared<VectorDouble>();
1042 auto g_ptr = boost::make_shared<VectorDouble>();
1044 auto eval_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1046 auto &pip = eval_fe_ptr->getOpPtrVector();
1052 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1060 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
1063 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
1066 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
1069 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
1073 auto parent_eval_fe_ptr =
1075 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1078 auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(mField);
1080 auto &pip = assemble_fe_ptr->getOpPtrVector();
1086 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1092 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
1095 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
1098 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
1101 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
1105 auto parent_assemble_fe_ptr =
1107 pip.push_back(create_run_parent_op(
1113 auto set_domain_lhs = [&](
auto fe) {
1116 auto &pip = fe->getOpPtrVector();
1124 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1133 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
1135 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
1138 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
1140 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
1143 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
1145 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
1148 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
1150 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
1157 auto set_bdy_rhs = [&](
auto fe) {
1159 auto &pip = fe->getOpPtrVector();
1161 auto l_ptr = boost::make_shared<VectorDouble>();
1163 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1165 auto &pip = eval_fe_ptr->getOpPtrVector();
1170 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1178 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
1182 auto parent_eval_fe_ptr =
1184 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1187 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1189 auto &pip = assemble_fe_ptr->getOpPtrVector();
1194 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1202 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1206 if (lPtr->size() != getGaussPts().size2()) {
1207 lPtr->resize(getGaussPts().size2());
1214 boost::shared_ptr<VectorDouble> lPtr;
1217 pip.push_back(
new OpLSize(l_ptr));
1219 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
1223 auto parent_assemble_fe_ptr =
1225 pip.push_back(create_run_parent_op(
1231 auto set_bdy_lhs = [&](
auto fe) {
1234 auto &pip = fe->getOpPtrVector();
1240 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1249 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
1251 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
1259 auto level_ents_ptr = boost::make_shared<Range>();
1263 auto get_prj_ents = [&]() {
1268 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1269 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1275 auto prj_ents = get_prj_ents();
1279 auto rebuild = [&]() {
1283 std::vector<std::string> fields{
"U",
"P",
"H",
"G",
"L"};
1284 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1286 {
"U", level_ents_ptr},
1287 {
"P", level_ents_ptr},
1288 {
"H", level_ents_ptr},
1289 {
"G", level_ents_ptr},
1290 {
"L", level_ents_ptr}
1294 CHKERR prb_mng->buildSubProblem(
"SUB_SOLVER", fields, fields,
1295 simple->getProblemName(), PETSC_TRUE,
1296 &range_maps, &range_maps);
1299 CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
1302 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
1307 MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
1311 auto dm =
simple->getDM();
1314 CHKERR DMSetType(subdm,
"DMMOFEM");
1319 MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
1324 auto create_dummy_dm = [&]() {
1327 simple->getProblemName().c_str(),
1333 auto subdm = create_subdm();
1336 pip_mng->getDomainRhsFE().reset();
1337 pip_mng->getDomainLhsFE().reset();
1338 pip_mng->getBoundaryRhsFE().reset();
1339 pip_mng->getBoundaryLhsFE().reset();
1344 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1345 pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1348 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1349 pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1353 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1354 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1355 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1356 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1361 auto ksp = pip_mng->createKSP(subdm);
1362 CHKERR KSPSetFromOptions(ksp);
1366 auto get_norm = [&](
auto x) {
1368 CHKERR VecNorm(x, NORM_2, &nrm);
1375 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1377 for (
auto &
v : ent_ptr->getEntFieldData()) {
1383 auto solve = [&](
auto S) {
1385 CHKERR VecZeroEntries(S);
1387 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
1388 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
1390 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1391 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1398 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
1403 auto dummy_dm = create_dummy_dm();
1408 auto apply_restrict = [&]() {
1409 auto get_is = [](
auto v) {
1411 auto create = [&]() {
1415 CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
1416 CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
1423 auto iy = get_is(glob_x);
1427 DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULL, dummy_dm),
1433 DMGetNamedGlobalVector(dummy_dm,
"TSTheta_Xdot", &Xdot),
1436 auto forward_ghost = [](
auto g) {
1438 CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
1439 CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
1446 if constexpr (
debug) {
1448 <<
"Reverse restrict: X0 " << get_norm(X0) <<
" Xdot "
1452 return std::vector<Vec>{X0, Xdot};
1455 auto ts_solver_vecs = apply_restrict();
1457 if (ts_solver_vecs.size()) {
1459 for (
auto v : ts_solver_vecs) {
1460 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1466 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1467 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1468 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
1476 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1479 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1480 &ts_solver_vecs[0]);
1481 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1482 &ts_solver_vecs[1]);
1485 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1486 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1487 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
1496 auto post_proc = [&](
auto exe_test) {
1498 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1499 auto &pip = post_proc_fe->getOpPtrVector();
1500 post_proc_fe->exeTestHook = exe_test;
1508 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1511 fe_parent->getOpPtrVector(), {H1});
1519 auto u_ptr = boost::make_shared<MatrixDouble>();
1520 auto p_ptr = boost::make_shared<VectorDouble>();
1521 auto h_ptr = boost::make_shared<VectorDouble>();
1522 auto g_ptr = boost::make_shared<VectorDouble>();
1524 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1527 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1530 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1533 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1537 post_proc_fe->getOpPtrVector().push_back(
1539 new OpPPMap(post_proc_fe->getPostProcMesh(),
1540 post_proc_fe->getMapGaussPts(),
1542 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1554 auto dm =
simple->getDM();
1556 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1565 if constexpr (
debug) {
1571 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1572 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1573 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1574 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1587 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
1592 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1600 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1618 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1619 auto u_ptr = boost::make_shared<MatrixDouble>();
1620 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1621 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1622 auto h_ptr = boost::make_shared<VectorDouble>();
1623 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1624 auto g_ptr = boost::make_shared<VectorDouble>();
1625 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1626 auto lambda_ptr = boost::make_shared<VectorDouble>();
1627 auto p_ptr = boost::make_shared<VectorDouble>();
1628 auto div_u_ptr = boost::make_shared<VectorDouble>();
1632 auto set_domain_general = [&](
auto fe) {
1634 auto &pip = fe->getOpPtrVector();
1642 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1645 fe_parent->getOpPtrVector(), {H1});
1651 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1671 "Coordinate system not implemented");
1674 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1680 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1685 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1690 auto set_domain_rhs = [&](
auto fe) {
1692 auto &pip = fe->getOpPtrVector();
1694 CHKERR set_domain_general(fe);
1696 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1697 pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1698 grad_h_ptr, g_ptr, p_ptr));
1700 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1701 pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1704 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1705 pip.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
1707 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1709 "P", div_u_ptr, [](
const double r,
const double,
const double) {
1713 "P", p_ptr, [](
const double r,
const double,
const double) {
1719 auto set_domain_lhs = [&](
auto fe) {
1721 auto &pip = fe->getOpPtrVector();
1723 CHKERR set_domain_general(fe);
1725 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1727 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1728 pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
1729 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1731 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1733 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1734 pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
1737 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1739 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1740 pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
1741 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1743 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1747 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1749 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1751 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1755 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1757 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1763 [](
const double r,
const double,
const double) {
1771 [](
const double r,
const double,
const double) {
1778 "Coordinate system not implemented");
1781 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
1782 pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
1790 auto get_block_name = [](
auto name) {
1791 return boost::format(
"%s(.*)") %
"WETTING_ANGLE";
1794 auto get_blocks = [&](
auto &&name) {
1796 std::regex(name.str()));
1799 auto set_boundary_rhs = [&](
auto fe) {
1801 auto &pip = fe->getOpPtrVector();
1807 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1814 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1817 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1822 pip,
mField,
"L", {},
"INFLUX",
1825 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1828 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1829 if (wetting_block.size()) {
1836 op_bdy_side->getOpPtrVector(), {H1});
1839 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1843 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1846 fe_parent->getOpPtrVector(), {H1});
1852 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1854 op_bdy_side->getOpPtrVector().push_back(
1857 pip.push_back(op_bdy_side);
1860 for (
auto &b : wetting_block) {
1862 std::vector<double> attr_vec;
1863 CHKERR b->getMeshsetIdEntitiesByDimension(
1865 b->getAttributes(attr_vec);
1866 if (attr_vec.size() != 1)
1868 "Should be one attribute");
1869 MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
1871 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1873 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1881 auto set_boundary_lhs = [&](
auto fe) {
1883 auto &pip = fe->getOpPtrVector();
1889 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1896 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1897 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
1900 auto wetting_block = get_blocks(get_block_name(
"WETTING_ANGLE"));
1901 if (wetting_block.size()) {
1902 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1903 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1911 op_bdy_side->getOpPtrVector(), {H1});
1914 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1918 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1921 fe_parent->getOpPtrVector(), {H1});
1927 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1929 op_bdy_side->getOpPtrVector().push_back(
1931 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1933 op_bdy_side->getOpPtrVector().push_back(
1937 pip.push_back(op_bdy_side);
1940 for (
auto &b : wetting_block) {
1942 std::vector<double> attr_vec;
1943 CHKERR b->getMeshsetIdEntitiesByDimension(
1945 b->getAttributes(attr_vec);
1946 if (attr_vec.size() != 1)
1948 "Should be one attribute");
1950 <<
"wetting angle edges size " << force_edges.size();
1952 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1954 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1955 boost::make_shared<Range>(force_edges), attr_vec.front()));
1964 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1965 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1966 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1967 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1991 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1992 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1993 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
2000 MOFEM_LOG(
"FS", Sev::verbose) <<
"Monitor";
2003 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
2006 auto post_proc_begin =
2007 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2009 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2017 CHKERR post_proc_end->writeFile(
2018 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
2019 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
2025 MPI_Allreduce(MPI_IN_PLACE, &(*
liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2029 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
2055 auto setup_subdm = [&](
auto dm) {
2059 auto dm =
simple->getDM();
2060 auto level_ents_ptr = boost::make_shared<Range>();
2061 CHKERR bit_mng->getEntitiesByRefLevel(
2064 CHKERR DMSetType(subdm,
"DMMOFEM");
2070 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
2083 auto get_fe_post_proc = [&](
auto post_proc_mesh) {
2084 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2089 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2098 boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
2099 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2104 auto u_ptr = boost::make_shared<MatrixDouble>();
2105 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2106 auto h_ptr = boost::make_shared<VectorDouble>();
2107 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2108 auto p_ptr = boost::make_shared<VectorDouble>();
2109 auto g_ptr = boost::make_shared<VectorDouble>();
2110 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2113 post_proc_fe->getOpPtrVector(), {H1});
2119 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2122 fe_parent->getOpPtrVector(), {H1});
2128 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
2129 post_proc_fe->getOpPtrVector().push_back(
2131 post_proc_fe->getOpPtrVector().push_back(
2135 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
2136 post_proc_fe->getOpPtrVector().push_back(
2138 post_proc_fe->getOpPtrVector().push_back(
2141 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
2142 post_proc_fe->getOpPtrVector().push_back(
2145 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
2146 post_proc_fe->getOpPtrVector().push_back(
2148 post_proc_fe->getOpPtrVector().push_back(
2153 post_proc_fe->getOpPtrVector().push_back(
2156 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2158 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
2160 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
2162 {{
"GRAD_U", grad_u_ptr}},
2170 return post_proc_fe;
2173 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh) {
2174 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2175 return setParentDofs(
2179 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2188 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2189 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2198 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2207 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2208 boost::shared_ptr<MatrixDouble> n_ptr)
2216 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2217 ptrNormal->resize(
SPACE_DIM, getGaussPts().size2(),
false);
2218 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2219 for (
auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2220 t_n(
i) = t_n_fe(
i) * t_l / std::sqrt(t_n_fe(
i) * t_n_fe(
i));
2229 boost::shared_ptr<VectorDouble> ptrL;
2230 boost::shared_ptr<MatrixDouble> ptrNormal;
2233 auto u_ptr = boost::make_shared<MatrixDouble>();
2234 auto p_ptr = boost::make_shared<VectorDouble>();
2235 auto lambda_ptr = boost::make_shared<VectorDouble>();
2236 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2238 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
2239 post_proc_fe->getOpPtrVector().push_back(
2242 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
2243 post_proc_fe->getOpPtrVector().push_back(
2245 post_proc_fe->getOpPtrVector().push_back(
2246 new OpGetNormal(lambda_ptr, normal_l_ptr));
2248 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
2249 post_proc_fe->getOpPtrVector().push_back(
2253 op_ptr->doWorkRhsHook = [&](
DataOperator *base_op_ptr,
int side,
2262 post_proc_fe->getOpPtrVector().push_back(
2264 new OpPPMap(post_proc_fe->getPostProcMesh(),
2265 post_proc_fe->getMapGaussPts(),
2279 return post_proc_fe;
2282 auto get_lift_fe = [&]() {
2283 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2284 return setParentDofs(
2288 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2296 auto fe = boost::make_shared<BoundaryEle>(mField);
2297 fe->exeTestHook = [](
FEMethod *fe_ptr) {
2302 auto lift_ptr = boost::make_shared<VectorDouble>();
2303 auto p_ptr = boost::make_shared<VectorDouble>();
2304 auto ents_ptr = boost::make_shared<Range>();
2310 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2317 std::vector<const CubitMeshSets *> vec_ptr;
2319 std::regex(
"LIFT"), vec_ptr);
2320 for (
auto m_ptr : vec_ptr) {
2325 ents_ptr->merge(ents);
2328 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
2330 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
2331 fe->getOpPtrVector().push_back(
2333 fe->getOpPtrVector().push_back(
2336 return std::make_pair(fe, lift_ptr);
2339 auto set_post_proc_monitor = [&](
auto ts) {
2343 boost::shared_ptr<FEMethod> null_fe;
2344 auto post_proc_mesh = boost::make_shared<moab::Core>();
2345 auto monitor_ptr = boost::make_shared<Monitor>(
2347 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2350 null_fe, monitor_ptr);
2354 auto dm =
simple->getDM();
2358 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2363 ptr->fsRawPtr =
this;
2364 ptr->solverSubDM = create_solver_dm(
simple->getDM());
2368 CHKERR VecAssemblyBegin(ptr->globSol);
2369 CHKERR VecAssemblyEnd(ptr->globSol);
2371 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2373 CHKERR set_post_proc_monitor(sub_ts);
2376 CHKERR TSSetFromOptions(ts);
2380 auto print_fields_in_section = [&]() {
2383 simple->getProblemName());
2384 PetscInt num_fields;
2385 CHKERR PetscSectionGetNumFields(section, &num_fields);
2386 for (
int f = 0;
f < num_fields; ++
f) {
2390 <<
"Field " <<
f <<
" " << std::string(
field_name);
2395 CHKERR print_fields_in_section();
2397 CHKERR TSSolve(ts, ptr->globSol);
2405 #ifdef PYTHON_INIT_SURFACE
2410 const char param_file[] =
"param_file.petsc";
2414 auto core_log = logging::core::get();
2415 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
"FS"));
2416 LogManager::setLog(
"FS");
2422 DMType dm_name =
"DMMOFEM";
2445 #ifdef PYTHON_INIT_SURFACE
2446 if (Py_FinalizeEx() < 0) {
2462 "can not get vertices on bit 0");
2467 Range plus_range, minus_range;
2468 std::vector<EntityHandle> plus, minus;
2471 for (
auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2473 const auto f = p->first;
2474 const auto s = p->second;
2477 const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number,
f);
2478 const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2479 auto it = dofs_mi.lower_bound(lo_uid);
2480 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2484 plus.reserve(std::distance(it, hi_it));
2485 minus.reserve(std::distance(it, hi_it));
2487 for (; it != hi_it; ++it) {
2488 const auto v = (*it)->getFieldData();
2490 plus.push_back((*it)->getEnt());
2492 minus.push_back((*it)->getEnt());
2495 plus_range.insert_list(plus.begin(), plus.end());
2496 minus_range.insert_list(minus.begin(), minus.end());
2501 <<
"Plus range " << plus_range << endl;
2503 <<
"Minus range " << minus_range << endl;
2506 auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
2509 moab::Interface::UNION),
2510 "can not get adjacencies");
2512 "can not filter elements with bit 0");
2516 CHKERR comm_mng->synchroniseEntities(plus_range);
2517 CHKERR comm_mng->synchroniseEntities(minus_range);
2520 ele_plus[0] = get_elems(plus_range,
bit(0),
BitRefLevel().set());
2521 ele_minus[0] = get_elems(minus_range,
bit(0),
BitRefLevel().set());
2522 auto common = intersect(ele_plus[0], ele_minus[0]);
2523 ele_plus[0] = subtract(ele_plus[0], common);
2524 ele_minus[0] = subtract(ele_minus[0], common);
2526 auto get_children = [&](
auto &p,
auto &
c) {
2528 CHKERR bit_mng->updateRangeByChildren(p,
c);
2540 auto get_level = [&](
auto &p,
auto &
m,
auto z,
auto bit,
auto mask) {
2543 bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
SPACE_DIM,
l),
2544 "can not get vertices on bit");
2547 for (
auto f = 0;
f != z; ++
f) {
2551 l = get_elems(conn,
bit, mask);
2556 std::vector<Range> vec_levels(
nb_levels);
2558 vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
2562 if constexpr (
debug) {
2565 std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
2582 start_mask[s] =
true;
2588 Range prev_level_skin;
2592 CHKERR bit_mng->lambdaBitRefLevel(
2599 auto set_levels = [&](
auto &&
2609 auto get_adj = [&](
auto ents) {
2615 ents.subset_by_dimension(
SPACE_DIM),
d,
false, conn,
2616 moab::Interface::UNION),
2630 CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
2633 level_prev = subtract(level_prev, parents);
2635 level_prev.merge(vec_levels[
l]);
2643 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2645 level = get_adj(level);
2654 auto set_skins = [&]() {
2658 ParallelComm *pcomm =
2667 "can't get bit level");
2674 auto get_level_skin = [&]() {
2680 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2681 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2682 PSTATUS_NOT, -1,
nullptr);
2685 skin_level = subtract(skin_level,
2689 Range skin_level_verts;
2692 skin_level.merge(skin_level_verts);
2695 bit_prev.set(
l - 1);
2699 skin.merge(subtract(skin_level, level_prev));
2705 auto resolve_shared = [&](
auto &&skin) {
2706 Range tmp_skin = skin;
2708 map<int, Range> map_procs;
2710 tmp_skin, &map_procs);
2712 Range from_other_procs;
2713 for (
auto &
m : map_procs) {
2715 from_other_procs.merge(
m.second);
2719 auto common = intersect(
2720 skin, from_other_procs);
2721 skin.merge(from_other_procs);
2725 if (!common.empty()) {
2727 skin = subtract(skin, common);
2734 auto get_parent_level_skin = [&](
auto skin) {
2736 CHKERR bit_mng->updateRangeByParent(
2737 skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
2738 Range skin_parent_verts;
2741 skin_parents.merge(skin_parent_verts);
2744 return skin_parents;
2747 auto child_skin = resolve_shared(get_level_skin());
2748 auto parent_skin = get_parent_level_skin(child_skin);
2750 child_skin = subtract(child_skin, parent_skin);
2758 auto set_current = [&]() {
2767 last_level = subtract(last_level, skin_child);
2779 if constexpr (
debug) {
2782 std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
2785 "PARALLEL=WRITE_PART");
2789 "MOAB",
"PARALLEL=WRITE_PART");
2792 "MOAB",
"PARALLEL=WRITE_PART");
2796 "MOAB",
"PARALLEL=WRITE_PART");
2799 "MOAB",
"PARALLEL=WRITE_PART");
2807 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
2810 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2817 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
2819 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
2824 auto fe_ptr_current = get_elem();
2828 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2833 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2836 parent_fe_pt->getOpPtrVector().push_back(
2840 H1, op, fe_ptr_current,
2851 parent_fe_pt->getOpPtrVector().push_back(
2866 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2881 auto cut_off_dofs = [&]() {
2884 auto &m_field = ptr->fsRawPtr->mField;
2886 Range current_verts;
2888 CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
2891 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
2893 for (
auto &
h : ent_ptr->getEntFieldData()) {
2899 auto field_blas = m_field.getInterface<
FieldBlas>();
2900 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts,
"H",
2909 MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
2911 auto &m_field = ptr->fsRawPtr->mField;
2915 auto get_norm = [&](
auto x) {
2917 CHKERR VecNorm(x, NORM_2, &nrm);
2922 auto refine_problem = [&]() {
2924 MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
2926 CHKERR ptr->fsRawPtr->projectData();
2932 auto set_jacobian_operators = [&]() {
2935 CHKERR KSPReset(ptr->subKSP);
2940 auto set_solution = [&]() {
2942 MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
2944 PetscObjectState state;
2955 INSERT_VALUES, SCATTER_FORWARD);
2958 <<
"Set solution, vector norm " << get_norm(ptr->globSol);
2963 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2967 CHKERR set_jacobian_operators();
2972 "Sorry, only TSTheta handling is implemented");
2977 PetscBarrier((PetscObject)ts);
2989 auto &m_field = ptr->fsRawPtr->mField;
3001 auto sub_u = ptr->getSubVector();
3004 auto scatter = ptr->getScatter(sub_u, u,
R);
3006 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3008 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3009 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3010 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3011 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3015 CHKERR apply_scatter_and_update(u, sub_u);
3016 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3019 CHKERR VecScatterBegin(scatter, sub_f,
f, INSERT_VALUES, SCATTER_FORWARD);
3020 CHKERR VecScatterEnd(scatter, sub_f,
f, INSERT_VALUES, SCATTER_FORWARD);
3026 PetscReal
a, Mat
A, Mat B,
3030 auto sub_u = ptr->getSubVector();
3032 auto scatter = ptr->getScatter(sub_u, u,
R);
3034 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
3036 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3037 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3038 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3039 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3043 CHKERR apply_scatter_and_update(u, sub_u);
3044 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3047 ptr->tsCtxPtr.get());
3056 auto get_norm = [&](
auto x) {
3058 CHKERR VecNorm(x, NORM_2, &nrm);
3062 auto sub_u = ptr->getSubVector();
3063 auto scatter = ptr->getScatter(sub_u, u,
R);
3064 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3065 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3066 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3067 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3070 <<
"u norm " << get_norm(u) <<
" u sub nom " << get_norm(sub_u);
3080 MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
3081 ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
3082 CHKERR KSPSetFromOptions(ptr->subKSP);
3090 auto sub_x = ptr->getSubVector();
3092 auto scatter = ptr->getScatter(sub_x, pc_x,
R);
3093 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3095 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3096 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3097 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
3098 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3099 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
3100 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3102 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3109 auto prb_ptr = ptr->fsRawPtr->mField.get_problem(
"SUB_SOLVER");
3110 if (
auto sub_data = prb_ptr->getSubData()) {
3111 auto is = sub_data->getSmartColIs();
3138 auto dm =
simple->getDM();
3146 CHKERR TSGetSNES(ts, &snes);
3149 auto set_section_monitor = [&](
auto snes) {
3151 CHKERR SNESMonitorSet(snes,
3154 (
void *)(snes_ctx_ptr.get()),
nullptr);
3158 CHKERR set_section_monitor(snes);
3160 auto ksp =
createKSP(m_field.get_comm());
3161 CHKERR KSPSetType(ksp, KSPPREONLY);
3163 CHKERR PCSetType(sub_pc, PCSHELL);
3166 CHKERR KSPSetPC(ksp, sub_pc);
3167 CHKERR SNESSetKSP(snes, ksp);