13#include <petsc/private/petscimpl.h>
17static char help[] =
"...\n\n";
19#ifdef PYTHON_INIT_SURFACE
20#include <boost/python.hpp>
21#include <boost/python/def.hpp>
22namespace bp = boost::python;
25 SurfacePython() =
default;
26 virtual ~SurfacePython() =
default;
33 auto main_module = bp::import(
"__main__");
34 mainNamespace = main_module.attr(
"__dict__");
35 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
37 surfaceFun = mainNamespace[
"surface"];
39 }
catch (bp::error_already_set
const &) {
49 double x,
double y,
double z,
double eta,
double &s
56 s = bp::extract<double>(surfaceFun(x, y, z,
eta));
58 }
catch (bp::error_already_set
const &) {
67 bp::object mainNamespace;
68 bp::object surfaceFun;
71static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
79 EXECUTABLE_COORD_TYPE;
84 IntegrationType::GAUSS;
96using SideOp = SideEle::UserDataOperator;
160constexpr double a0 = 0;
162constexpr double mu_m = 0.010101;
164constexpr double mu_p = 0.000182;
166constexpr double W = 0.25;
169constexpr double h = 0.0015 / 20;
174constexpr double md = 1e-2;
175constexpr double eps = 1e-12;
176constexpr double tol = std::numeric_limits<float>::epsilon();
196 constexpr int sub_stepping = 16;
197 return std::min(1.,
static_cast<double>(ts_step) / sub_stepping);
202auto my_max = [](
const double x) {
return (x - 1 + std::abs(x + 1)) / 2; };
203auto my_min = [](
const double x) {
return (x + 1 - std::abs(x - 1)) / 2; };
206 if (
h >= -1 &&
h < 1)
222auto get_f = [](
const double h) {
return 4 *
W *
h * (
h *
h - 1); };
230 return md * (1 -
h *
h);
240 const double h2 =
h *
h;
241 const double h3 = h2 *
h;
243 return md * (2 * h3 - 3 * h2 + 1);
245 return md * (-2 * h3 - 3 * h2 + 1);
269 constexpr double R = 0.0125;
270 constexpr double A =
R * 0.2;
271 const double theta = atan2(r, y);
272 const double w =
R +
A * cos(
n * theta);
273 const double d = std::sqrt(r * r + y * y);
274 return tanh((w - d) / (
eta * std::sqrt(2)));
278 constexpr double y0 = 0.5;
279 constexpr double R = 0.4;
281 const double d = std::sqrt(r * r + y * y);
282 return tanh((
R - d) / (
eta * std::sqrt(2)));
286 constexpr double water_height = 0.;
287 return tanh((water_height - y) / (
eta * std::sqrt(2)));
292 return -tanh((-0.039 - x) / (
eta * std::sqrt(2)));
302auto init_h = [](
double r,
double y,
double theta) {
303#ifdef PYTHON_INIT_SURFACE
305 if (
auto ptr = surfacePythonWeakPtr.lock()) {
307 "error eval python");
322 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
327 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
331auto save_range = [](moab::Interface &moab,
const std::string name,
336 CHKERR moab.add_entities(*out_meshset, r);
337 CHKERR moab.write_file(name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART",
338 out_meshset->get_ptr(), 1);
349 std::vector<EntityHandle> ents_vec;
355 auto dofs = prb_ptr->numeredRowDofsPtr;
358 ents_vec.reserve(std::distance(lo_it, hi_it));
360 for (; lo_it != hi_it; ++lo_it) {
361 ents_vec.push_back((*lo_it)->getEnt());
364 std::sort(ents_vec.begin(), ents_vec.end());
365 auto it = std::unique(ents_vec.begin(), ents_vec.end());
367 r.insert_list(ents_vec.begin(), it);
377 std::vector<EntityHandle> ents_vec;
382 auto dofs = prb_ptr->numeredRowDofsPtr;
383 ents_vec.reserve(dofs->size());
385 for (
auto d : *dofs) {
386 ents_vec.push_back(d->getEnt());
389 std::sort(ents_vec.begin(), ents_vec.end());
390 auto it = std::unique(ents_vec.begin(), ents_vec.end());
392 r.insert_list(ents_vec.begin(), it);
457 PetscReal
a, Mat
A, Mat
B,
468 boost::shared_ptr<SnesCtx>
470 boost::shared_ptr<TsCtx>
519 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
520 ForcesAndSourcesCore::UserDataOperator::OpType op,
522 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
544 <<
"Read mesh for problem in " << EXECUTABLE_COORD_TYPE;
547 simple->getParentAdjacencies() =
true;
600 auto set_problem_bit = [&]() {
623#ifdef PYTHON_INIT_SURFACE
624 auto get_py_surface_init = []() {
625 auto py_surf_init = boost::make_shared<SurfacePython>();
626 CHKERR py_surf_init->surfaceInit(
"surface.py");
627 surfacePythonWeakPtr = py_surf_init;
630 auto py_surf_init = get_py_surface_init();
637 auto dm =
simple->getDM();
641 auto reset_bits = [&]() {
645 start_mask[s] =
true;
647 CHKERR bit_mng->lambdaBitRefLevel(
656 auto add_parent_field = [&](
auto fe,
auto op,
auto field) {
661 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
669 auto h_ptr = boost::make_shared<VectorDouble>();
670 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
671 auto g_ptr = boost::make_shared<VectorDouble>();
672 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
674 auto set_generic = [&](
auto fe) {
676 auto &pip = fe->getOpPtrVector();
684 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
687 fe_parent->getOpPtrVector(), {H1});
693 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
698 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
706 auto post_proc = [&](
auto exe_test) {
708 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
709 post_proc_fe->exeTestHook = exe_test;
711 CHKERR set_generic(post_proc_fe);
715 post_proc_fe->getOpPtrVector().push_back(
717 new OpPPMap(post_proc_fe->getPostProcMesh(),
718 post_proc_fe->getMapGaussPts(),
720 {{
"H", h_ptr}, {
"G", g_ptr}},
722 {{
"GRAD_H", grad_h_ptr}, {
"GRAD_G", grad_g_ptr}},
733 CHKERR post_proc_fe->writeFile(
"out_init.h5m");
738 auto solve_init = [&](
auto exe_test) {
741 pip_mng->getOpDomainRhsPipeline().clear();
742 pip_mng->getOpDomainLhsPipeline().clear();
744 auto set_domain_rhs = [&](
auto fe) {
747 auto &pip = fe->getOpPtrVector();
749 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
750 pip.push_back(
new OpRhsH<true>(
"H",
nullptr,
nullptr, h_ptr, grad_h_ptr,
752 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
753 pip.push_back(
new OpRhsG<true>(
"G", h_ptr, grad_h_ptr, g_ptr));
757 auto set_domain_lhs = [&](
auto fe) {
761 auto &pip = fe->getOpPtrVector();
763 CHKERR add_parent_field(fe, UDO::OPROW,
"H");
764 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
767 CHKERR add_parent_field(fe, UDO::OPCOL,
"G");
770 CHKERR add_parent_field(fe, UDO::OPROW,
"G");
773 CHKERR add_parent_field(fe, UDO::OPCOL,
"H");
778 auto create_subdm = [&]() {
779 auto level_ents_ptr = boost::make_shared<Range>();
785 CHKERR DMSetType(subdm,
"DMMOFEM");
791 for (
auto f : {
"H",
"G"}) {
797 if constexpr (
debug) {
801 dm_ents.subset_by_type(MBVERTEX));
802 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
810 auto subdm = create_subdm();
812 pip_mng->getDomainRhsFE().reset();
813 pip_mng->getDomainLhsFE().reset();
816 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
817 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
819 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
820 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
823 auto snes = pip_mng->createSNES(subdm);
826 auto set_section_monitor = [&](
auto solver) {
828 CHKERR SNESMonitorSet(solver,
831 (
void *)(snes_ctx_ptr.get()),
nullptr);
836 auto print_section_field = [&]() {
842 CHKERR PetscSectionGetNumFields(section, &num_fields);
843 for (
int f = 0;
f < num_fields; ++
f) {
847 <<
"Field " <<
f <<
" " << std::string(
field_name);
853 CHKERR set_section_monitor(snes);
854 CHKERR print_section_field();
856 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
860 CHKERR SNESSetFromOptions(snes);
861 CHKERR SNESSolve(snes, PETSC_NULL,
D);
863 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
864 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
874 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
884 pip_mng->getOpDomainRhsPipeline().clear();
885 pip_mng->getOpDomainLhsPipeline().clear();
888 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
890 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_X",
892 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
894 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM_Y",
896 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"U",
898 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX",
"L",
900 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"ZERO",
922 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
923 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
924 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
925 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
927 pip_mng->getDomainLhsFE().reset();
928 pip_mng->getDomainRhsFE().reset();
929 pip_mng->getBoundaryLhsFE().reset();
930 pip_mng->getBoundaryRhsFE().reset();
935 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
940 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
949 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
954 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
964 auto create_run_parent_op = [&](
auto parent_fe_ptr,
auto this_fe_ptr,
966 auto parent_mask = fe_bit;
974 auto get_parents_vel_fe_ptr = [&](
auto this_fe_ptr,
auto fe_bit) {
975 std::vector<boost::shared_ptr<DomianParentEle>> parents_elems_ptr_vec;
977 parents_elems_ptr_vec.emplace_back(
978 boost::make_shared<DomianParentEle>(
mField));
980 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
981 create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
983 return parents_elems_ptr_vec[0];
987 auto solve_projection = [&](
auto exe_test) {
990 auto set_domain_rhs = [&](
auto fe) {
992 auto &pip = fe->getOpPtrVector();
994 auto u_ptr = boost::make_shared<MatrixDouble>();
995 auto p_ptr = boost::make_shared<VectorDouble>();
996 auto h_ptr = boost::make_shared<VectorDouble>();
997 auto g_ptr = boost::make_shared<VectorDouble>();
999 auto eval_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1001 auto &pip = eval_fe_ptr->getOpPtrVector();
1007 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1015 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
1018 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
1021 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
1024 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
1028 auto parent_eval_fe_ptr =
1030 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1033 auto assemble_fe_ptr = boost::make_shared<DomianParentEle>(
mField);
1035 auto &pip = assemble_fe_ptr->getOpPtrVector();
1041 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1047 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
1050 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
1053 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
1056 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
1060 auto parent_assemble_fe_ptr =
1062 pip.push_back(create_run_parent_op(
1068 auto set_domain_lhs = [&](
auto fe) {
1071 auto &pip = fe->getOpPtrVector();
1079 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1088 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
1090 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
1093 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
1095 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
1098 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
1100 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
1103 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
1105 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
1112 auto set_bdy_rhs = [&](
auto fe) {
1114 auto &pip = fe->getOpPtrVector();
1116 auto l_ptr = boost::make_shared<VectorDouble>();
1118 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1120 auto &pip = eval_fe_ptr->getOpPtrVector();
1125 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1133 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
1137 auto parent_eval_fe_ptr =
1139 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1142 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1144 auto &pip = assemble_fe_ptr->getOpPtrVector();
1149 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1157 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1161 if (lPtr->size() != getGaussPts().size2()) {
1162 lPtr->resize(getGaussPts().size2());
1169 boost::shared_ptr<VectorDouble> lPtr;
1172 pip.push_back(
new OpLSize(l_ptr));
1174 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
1178 auto parent_assemble_fe_ptr =
1180 pip.push_back(create_run_parent_op(
1186 auto set_bdy_lhs = [&](
auto fe) {
1189 auto &pip = fe->getOpPtrVector();
1195 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1204 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
1206 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
1214 auto level_ents_ptr = boost::make_shared<Range>();
1218 auto get_prj_ents = [&]() {
1223 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1224 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1230 auto prj_ents = get_prj_ents();
1234 auto rebuild = [&]() {
1238 std::vector<std::string> fields{
"U",
"P",
"H",
"G",
"L"};
1239 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1241 {
"U", level_ents_ptr},
1242 {
"P", level_ents_ptr},
1243 {
"H", level_ents_ptr},
1244 {
"G", level_ents_ptr},
1245 {
"L", level_ents_ptr}
1249 CHKERR prb_mng->buildSubProblem(
"SUB_SOLVER", fields, fields,
1250 simple->getProblemName(), PETSC_TRUE,
1251 &range_maps, &range_maps);
1254 CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
1257 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
1261 MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
1265 auto dm =
simple->getDM();
1268 CHKERR DMSetType(subdm,
"DMMOFEM");
1273 MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
1278 auto create_dummy_dm = [&]() {
1281 simple->getProblemName().c_str(),
1287 auto subdm = create_subdm();
1290 pip_mng->getDomainRhsFE().reset();
1291 pip_mng->getDomainLhsFE().reset();
1292 pip_mng->getBoundaryRhsFE().reset();
1293 pip_mng->getBoundaryLhsFE().reset();
1298 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1299 pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1302 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1303 pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1307 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1308 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1309 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1310 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1315 auto ksp = pip_mng->createKSP(subdm);
1316 CHKERR KSPSetFromOptions(ksp);
1320 auto get_norm = [&](
auto x) {
1322 CHKERR VecNorm(x, NORM_2, &nrm);
1329 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1331 for (
auto &
v : ent_ptr->getEntFieldData()) {
1341 auto cut_off_dofs = [&]() {
1344 Range current_verts;
1347 MBVERTEX, current_verts);
1349 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
1351 for (
auto &
h : ent_ptr->getEntFieldData()) {
1358 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts,
"H",
1363 auto solve = [&](
auto S) {
1365 CHKERR VecZeroEntries(S);
1367 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
1368 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
1370 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1371 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1375 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
1380 auto dummy_dm = create_dummy_dm();
1385 auto apply_restrict = [&]() {
1386 auto get_is = [](
auto v) {
1388 auto create = [&]() {
1392 CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
1393 CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
1400 auto iy = get_is(glob_x);
1404 DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULL, dummy_dm),
1410 DMGetNamedGlobalVector(dummy_dm,
"TSTheta_Xdot", &Xdot),
1413 auto forward_ghost = [](
auto g) {
1415 CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
1416 CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
1423 if constexpr (
debug) {
1425 <<
"Reverse restrict: X0 " << get_norm(X0) <<
" Xdot "
1429 return std::vector<Vec>{X0, Xdot};
1432 auto ts_solver_vecs = apply_restrict();
1434 if (ts_solver_vecs.size()) {
1436 for (
auto v : ts_solver_vecs) {
1437 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1443 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1444 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1445 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1452 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1455 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1456 &ts_solver_vecs[0]);
1457 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1458 &ts_solver_vecs[1]);
1461 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1462 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " << f;
1463 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1473 auto post_proc = [&](
auto exe_test) {
1475 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1476 auto &pip = post_proc_fe->getOpPtrVector();
1477 post_proc_fe->exeTestHook = exe_test;
1485 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1488 fe_parent->getOpPtrVector(), {H1});
1496 auto u_ptr = boost::make_shared<MatrixDouble>();
1497 auto p_ptr = boost::make_shared<VectorDouble>();
1498 auto h_ptr = boost::make_shared<VectorDouble>();
1499 auto g_ptr = boost::make_shared<VectorDouble>();
1501 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1504 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1507 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1510 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1514 post_proc_fe->getOpPtrVector().push_back(
1516 new OpPPMap(post_proc_fe->getPostProcMesh(),
1517 post_proc_fe->getMapGaussPts(),
1519 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1531 auto dm =
simple->getDM();
1533 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1542 if constexpr (
debug) {
1548 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1549 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1550 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1551 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
1564 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
1569 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1577 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
1582 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1595 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
1596 auto u_ptr = boost::make_shared<MatrixDouble>();
1597 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
1598 auto dot_h_ptr = boost::make_shared<VectorDouble>();
1599 auto h_ptr = boost::make_shared<VectorDouble>();
1600 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1601 auto g_ptr = boost::make_shared<VectorDouble>();
1602 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1603 auto lambda_ptr = boost::make_shared<VectorDouble>();
1604 auto p_ptr = boost::make_shared<VectorDouble>();
1605 auto div_u_ptr = boost::make_shared<VectorDouble>();
1609 auto set_domain_general = [&](
auto fe) {
1611 auto &pip = fe->getOpPtrVector();
1619 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1622 fe_parent->getOpPtrVector(), {H1});
1628 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1638 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1644 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1649 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1654 auto set_domain_rhs = [&](
auto fe) {
1656 auto &pip = fe->getOpPtrVector();
1658 CHKERR set_domain_general(fe);
1660 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1661 pip.push_back(
new OpRhsU(
"U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
1662 grad_h_ptr, g_ptr, p_ptr));
1664 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1665 pip.push_back(
new OpRhsH<false>(
"H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
1668 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1669 pip.push_back(
new OpRhsG<false>(
"G", h_ptr, grad_h_ptr, g_ptr));
1671 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1673 "P", div_u_ptr, [](
const double r,
const double,
const double) {
1677 "P", p_ptr, [](
const double r,
const double,
const double) {
1683 auto set_domain_lhs = [&](
auto fe) {
1685 auto &pip = fe->getOpPtrVector();
1687 CHKERR set_domain_general(fe);
1689 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U");
1691 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1692 pip.push_back(
new OpLhsU_dU(
"U", u_ptr, grad_u_ptr, h_ptr));
1693 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1695 new OpLhsU_dH(
"U",
"H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
1697 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1698 pip.push_back(
new OpLhsU_dG(
"U",
"G", grad_h_ptr));
1701 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H");
1703 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1704 pip.push_back(
new OpLhsH_dU(
"H",
"U", grad_h_ptr));
1705 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1707 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1711 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G");
1713 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H");
1715 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G");
1719 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P");
1721 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U");
1724 [](
const double r,
const double,
const double) {
1728 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P");
1729 pip.push_back(
new OpDomainMassP(
"P",
"P", [](
double r,
double,
double) {
1737 auto get_block_name = [](
auto name) {
1738 return boost::format(
"%s(.*)") %
"WETTING_ANGLE";
1741 auto get_blocks = [&](
auto &&name) {
1743 std::regex(name.str()));
1746 auto set_boundary_rhs = [&](
auto fe) {
1748 auto &pip = fe->getOpPtrVector();
1754 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1761 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1764 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1769 pip,
mField,
"L", {},
"INFLUX", Sev::inform);
1771 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"U");
1780 op_bdy_side->getOpPtrVector(), {H1});
1783 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1787 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1790 fe_parent->getOpPtrVector(), {H1});
1796 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1798 op_bdy_side->getOpPtrVector().push_back(
1801 pip.push_back(op_bdy_side);
1804 for (
auto &b : get_blocks(get_block_name(
"WETTING_ANGLE"))) {
1806 std::vector<double> attr_vec;
1807 CHKERR b->getMeshsetIdEntitiesByDimension(
1809 b->getAttributes(attr_vec);
1810 if (attr_vec.size() != 1)
1812 MOFEM_LOG(
"FS", Sev::inform) <<
"Wetting angle: " << attr_vec.front();
1814 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1816 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
1823 auto set_boundary_lhs = [&](
auto fe) {
1825 auto &pip = fe->getOpPtrVector();
1831 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1838 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L");
1839 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"U");
1842 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
1843 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
1851 op_bdy_side->getOpPtrVector(), {H1});
1854 op_bdy_side->getSideFEPtr(),
"", UDO::OPSPACE,
1858 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1861 fe_parent->getOpPtrVector(), {H1});
1867 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
1869 op_bdy_side->getOpPtrVector().push_back(
1871 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
1873 op_bdy_side->getOpPtrVector().push_back(
1877 pip.push_back(op_bdy_side);
1880 for (
auto &b : get_blocks(get_block_name(
"WETTING_ANGLE"))) {
1882 std::vector<double> attr_vec;
1883 CHKERR b->getMeshsetIdEntitiesByDimension(
1885 b->getAttributes(attr_vec);
1886 if (attr_vec.size() != 1)
1889 <<
"wetting angle edges size " << force_edges.size();
1891 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"G");
1893 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
1894 boost::make_shared<Range>(force_edges), attr_vec.front()));
1902 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1903 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1904 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
1905 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
1929 boost::shared_ptr<PostProcEleDomainCont> post_proc,
1930 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
1931 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
1938 MOFEM_LOG(
"FS", Sev::verbose) <<
"Monitor";
1941 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc";
1944 auto post_proc_begin =
1945 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
1947 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1955 CHKERR post_proc_end->writeFile(
1956 "out_step_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
1957 MOFEM_LOG(
"FS", Sev::verbose) <<
"Mesh pre proc done";
1963 MPI_Allreduce(MPI_IN_PLACE, &(*
liftVec)[0],
SPACE_DIM, MPI_DOUBLE, MPI_SUM,
1967 <<
" lift vec x: " << (*liftVec)[0] <<
" y: " << (*liftVec)[1];
1993 auto setup_subdm = [&](
auto dm) {
1997 auto dm =
simple->getDM();
1998 auto level_ents_ptr = boost::make_shared<Range>();
1999 CHKERR bit_mng->getEntitiesByRefLevel(
2002 CHKERR DMSetType(subdm,
"DMMOFEM");
2008 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
2021 auto get_fe_post_proc = [&](
auto post_proc_mesh) {
2022 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field) {
2027 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2036 boost::make_shared<PostProcEleDomainCont>(
mField, post_proc_mesh);
2037 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2042 auto u_ptr = boost::make_shared<MatrixDouble>();
2043 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2044 auto h_ptr = boost::make_shared<VectorDouble>();
2045 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2046 auto p_ptr = boost::make_shared<VectorDouble>();
2047 auto g_ptr = boost::make_shared<VectorDouble>();
2048 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2051 post_proc_fe->getOpPtrVector(), {H1});
2057 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2060 fe_parent->getOpPtrVector(), {H1});
2066 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U");
2067 post_proc_fe->getOpPtrVector().push_back(
2069 post_proc_fe->getOpPtrVector().push_back(
2073 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H");
2074 post_proc_fe->getOpPtrVector().push_back(
2076 post_proc_fe->getOpPtrVector().push_back(
2079 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P");
2080 post_proc_fe->getOpPtrVector().push_back(
2083 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G");
2084 post_proc_fe->getOpPtrVector().push_back(
2086 post_proc_fe->getOpPtrVector().push_back(
2091 post_proc_fe->getOpPtrVector().push_back(
2094 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2096 {{
"H", h_ptr}, {
"P", p_ptr}, {
"G", g_ptr}},
2098 {{
"U", u_ptr}, {
"H_GRAD", grad_h_ptr}, {
"G_GRAD", grad_g_ptr}},
2100 {{
"GRAD_U", grad_u_ptr}},
2108 return post_proc_fe;
2111 auto get_bdy_post_proc_fe = [&](
auto post_proc_mesh) {
2112 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2113 return setParentDofs(
2117 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2126 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2127 post_proc_fe->exeTestHook = [](
FEMethod *fe_ptr) {
2136 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2143 auto u_ptr = boost::make_shared<MatrixDouble>();
2144 auto p_ptr = boost::make_shared<VectorDouble>();
2145 auto lambda_ptr = boost::make_shared<VectorDouble>();
2147 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"U");
2148 post_proc_fe->getOpPtrVector().push_back(
2151 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"L");
2152 post_proc_fe->getOpPtrVector().push_back(
2155 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW,
"P");
2156 post_proc_fe->getOpPtrVector().push_back(
2161 post_proc_fe->getOpPtrVector().push_back(
2163 new OpPPMap(post_proc_fe->getPostProcMesh(),
2164 post_proc_fe->getMapGaussPts(),
2178 return post_proc_fe;
2181 auto get_lift_fe = [&]() {
2182 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field) {
2183 return setParentDofs(
2187 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2195 auto fe = boost::make_shared<BoundaryEle>(mField);
2196 fe->exeTestHook = [](
FEMethod *fe_ptr) {
2201 auto lift_ptr = boost::make_shared<VectorDouble>();
2202 auto p_ptr = boost::make_shared<VectorDouble>();
2203 auto ents_ptr = boost::make_shared<Range>();
2209 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2216 std::vector<const CubitMeshSets *> vec_ptr;
2218 std::regex(
"LIFT"), vec_ptr);
2219 for (
auto m_ptr : vec_ptr) {
2224 ents_ptr->merge(ents);
2227 MOFEM_LOG(
"FS", Sev::noisy) <<
"Lift ents " << (*ents_ptr);
2229 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"P");
2230 fe->getOpPtrVector().push_back(
2232 fe->getOpPtrVector().push_back(
2235 return std::make_pair(fe, lift_ptr);
2238 auto set_post_proc_monitor = [&](
auto ts) {
2242 boost::shared_ptr<FEMethod> null_fe;
2243 auto post_proc_mesh = boost::make_shared<moab::Core>();
2244 auto monitor_ptr = boost::make_shared<Monitor>(
2246 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2249 null_fe, monitor_ptr);
2253 auto dm =
simple->getDM();
2257 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2262 ptr->fsRawPtr =
this;
2263 ptr->solverSubDM = create_solver_dm(
simple->getDM());
2267 CHKERR VecAssemblyBegin(ptr->globSol);
2268 CHKERR VecAssemblyEnd(ptr->globSol);
2270 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2272 CHKERR set_post_proc_monitor(sub_ts);
2278 CHKERR TSSetFromOptions(ts);
2282 auto print_fields_in_section = [&]() {
2285 simple->getProblemName());
2286 PetscInt num_fields;
2287 CHKERR PetscSectionGetNumFields(section, &num_fields);
2288 for (
int f = 0;
f < num_fields; ++
f) {
2292 <<
"Field " <<
f <<
" " << std::string(
field_name);
2297 CHKERR print_fields_in_section();
2299 CHKERR TSSolve(ts, ptr->globSol);
2307#ifdef PYTHON_INIT_SURFACE
2312 const char param_file[] =
"param_file.petsc";
2316 auto core_log = logging::core::get();
2324 DMType dm_name =
"DMMOFEM";
2329 moab::Core mb_instance;
2330 moab::Interface &moab = mb_instance;
2347#ifdef PYTHON_INIT_SURFACE
2348 if (Py_FinalizeEx() < 0) {
2364 "can not get vertices on bit 0");
2369 Range plus_range, minus_range;
2370 std::vector<EntityHandle> plus, minus;
2373 for (
auto p = vertices.pair_begin();
p != vertices.pair_end(); ++
p) {
2375 const auto f =
p->first;
2376 const auto s =
p->second;
2381 auto it = dofs_mi.lower_bound(lo_uid);
2382 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2386 plus.reserve(std::distance(it, hi_it));
2387 minus.reserve(std::distance(it, hi_it));
2389 for (; it != hi_it; ++it) {
2390 const auto v = (*it)->getFieldData();
2392 plus.push_back((*it)->getEnt());
2394 minus.push_back((*it)->getEnt());
2397 plus_range.insert_list(plus.begin(), plus.end());
2398 minus_range.insert_list(minus.begin(), minus.end());
2403 <<
"Plus range " << plus_range << endl;
2405 <<
"Minus range " << minus_range << endl;
2408 auto get_elems = [&](
auto &ents,
auto bit,
auto mask) {
2411 moab::Interface::UNION),
2412 "can not get adjacencies");
2414 "can not filter elements with bit 0");
2418 CHKERR comm_mng->synchroniseEntities(plus_range);
2419 CHKERR comm_mng->synchroniseEntities(minus_range);
2422 ele_plus[0] = get_elems(plus_range,
bit(0),
BitRefLevel().set());
2423 ele_minus[0] = get_elems(minus_range,
bit(0),
BitRefLevel().set());
2424 auto common = intersect(ele_plus[0], ele_minus[0]);
2425 ele_plus[0] = subtract(ele_plus[0], common);
2426 ele_minus[0] = subtract(ele_minus[0], common);
2428 auto get_children = [&](
auto &
p,
auto &
c) {
2430 CHKERR bit_mng->updateRangeByChildren(
p,
c);
2442 auto get_level = [&](
auto &
p,
auto &
m,
auto z,
auto bit,
auto mask) {
2445 bit_mng->getEntitiesByDimAndRefLevel(
bit, mask,
SPACE_DIM,
l),
2446 "can not get vertices on bit");
2449 for (
auto f = 0; f != z; ++f) {
2453 l = get_elems(conn,
bit, mask);
2458 std::vector<Range> vec_levels(
nb_levels);
2460 vec_levels[
l] = get_level(ele_plus[
l], ele_minus[
l], 2 * overlap,
bit(
l),
2464 if constexpr (
debug) {
2467 std::string name = (boost::format(
"out_r%d.h5m") %
l).str();
2486 start_mask[s] =
true;
2492 Range prev_level_skin;
2496 CHKERR bit_mng->lambdaBitRefLevel(
2503 auto set_levels = [&](
auto &&
2513 auto get_adj = [&](
auto ents) {
2519 ents.subset_by_dimension(
SPACE_DIM), d,
false, conn,
2520 moab::Interface::UNION),
2534 CHKERR bit_mng->updateRangeByParent(vec_levels[
l], parents);
2537 level_prev = subtract(level_prev, parents);
2539 level_prev.merge(vec_levels[
l]);
2547 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
2549 level = get_adj(level);
2558 auto set_skins = [&]() {
2562 ParallelComm *pcomm =
2571 "can't get bit level");
2578 auto get_level_skin = [&]() {
2584 CHKERR pcomm->filter_pstatus(skin_level_mesh,
2585 PSTATUS_SHARED | PSTATUS_MULTISHARED,
2586 PSTATUS_NOT, -1,
nullptr);
2589 skin_level = subtract(skin_level,
2593 Range skin_level_verts;
2596 skin_level.merge(skin_level_verts);
2599 bit_prev.set(
l - 1);
2603 skin.merge(subtract(skin_level, level_prev));
2609 auto resolve_shared = [&](
auto &&skin) {
2610 Range tmp_skin = skin;
2612 map<int, Range> map_procs;
2614 tmp_skin, &map_procs);
2616 Range from_other_procs;
2617 for (
auto &
m : map_procs) {
2619 from_other_procs.merge(
m.second);
2623 auto common = intersect(
2624 skin, from_other_procs);
2625 skin.merge(from_other_procs);
2629 if (!common.empty()) {
2631 skin = subtract(skin, common);
2638 auto get_parent_level_skin = [&](
auto skin) {
2640 CHKERR bit_mng->updateRangeByParent(
2641 skin.subset_by_dimension(
SPACE_DIM - 1), skin_parents);
2642 Range skin_parent_verts;
2645 skin_parents.merge(skin_parent_verts);
2648 return skin_parents;
2651 auto child_skin = resolve_shared(get_level_skin());
2652 auto parent_skin = get_parent_level_skin(child_skin);
2654 child_skin = subtract(child_skin, parent_skin);
2662 auto set_current = [&]() {
2671 last_level = subtract(last_level, skin_child);
2683 if constexpr (
debug) {
2686 std::string name = (boost::format(
"out_level%d.h5m") %
l).str();
2689 "PARALLEL=WRITE_PART");
2693 "MOAB",
"PARALLEL=WRITE_PART");
2696 "MOAB",
"PARALLEL=WRITE_PART");
2700 "MOAB",
"PARALLEL=WRITE_PART");
2703 "MOAB",
"PARALLEL=WRITE_PART");
2711 boost::shared_ptr<FEMethod> fe_top, std::string
field_name,
2712 ForcesAndSourcesCore::UserDataOperator::OpType op,
2714 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
2721 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>,
int)>
2723 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt,
int level) {
2728 auto fe_ptr_current = get_elem();
2732 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
2737 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
2740 parent_fe_pt->getOpPtrVector().push_back(
2744 H1, op, fe_ptr_current,
2755 parent_fe_pt->getOpPtrVector().push_back(
2770 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
2780 MOFEM_LOG(
"FS", Sev::inform) <<
"Run step pre proc";
2782 auto &m_field = ptr->fsRawPtr->mField;
2792 auto get_norm = [&](
auto x) {
2794 CHKERR VecNorm(x, NORM_2, &nrm);
2799 auto refine_problem = [&]() {
2801 MOFEM_LOG(
"FS", Sev::inform) <<
"Refine problem";
2803 CHKERR ptr->fsRawPtr->projectData();
2809 auto set_jacobian_operators = [&]() {
2812 CHKERR KSPReset(ptr->subKSP);
2817 auto set_solution = [&]() {
2819 MOFEM_LOG(
"FS", Sev::inform) <<
"Set solution";
2821 PetscObjectState state;
2832 INSERT_VALUES, SCATTER_FORWARD);
2835 <<
"Set solution, vector norm " << get_norm(ptr->globSol);
2840 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
2844 CHKERR set_jacobian_operators();
2849 "Sorry, only TSTheta handling is implemented");
2854 PetscBarrier((PetscObject)ts);
2866 auto &m_field = ptr->fsRawPtr->mField;
2878 auto sub_u = ptr->getSubVector();
2881 auto scatter = ptr->getScatter(sub_u, u,
R);
2883 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
2885 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2886 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2887 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2888 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2892 CHKERR apply_scatter_and_update(u, sub_u);
2893 CHKERR apply_scatter_and_update(u_t, sub_u_t);
2896 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
2897 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
2903 PetscReal
a, Mat
A, Mat
B,
2907 auto sub_u = ptr->getSubVector();
2909 auto scatter = ptr->getScatter(sub_u, u,
R);
2911 auto apply_scatter_and_update = [&](
auto x,
auto sub_x) {
2913 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2914 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
2915 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2916 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
2920 CHKERR apply_scatter_and_update(u, sub_u);
2921 CHKERR apply_scatter_and_update(u_t, sub_u_t);
2924 ptr->tsCtxPtr.get());
2933 auto get_norm = [&](
auto x) {
2935 CHKERR VecNorm(x, NORM_2, &nrm);
2939 auto sub_u = ptr->getSubVector();
2940 auto scatter = ptr->getScatter(sub_u, u,
R);
2941 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
2942 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
2943 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
2944 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
2947 <<
"u norm " << get_norm(u) <<
" u sub nom " << get_norm(sub_u);
2957 MOFEM_LOG(
"FS", Sev::verbose) <<
"SetUP sub PC";
2958 ptr->subKSP =
createKSP(ptr->fsRawPtr->mField.get_comm());
2959 CHKERR KSPSetFromOptions(ptr->subKSP);
2967 auto sub_x = ptr->getSubVector();
2969 auto scatter = ptr->getScatter(sub_x, pc_x,
R);
2970 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
2972 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
2973 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
2974 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve";
2975 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
2976 MOFEM_LOG(
"FS", Sev::verbose) <<
"PCShell solve <- done";
2977 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
2979 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
2986 auto prb_ptr = ptr->fsRawPtr->mField.get_problem(
"SUB_SOLVER");
2987 if (
auto sub_data = prb_ptr->getSubData()) {
2988 auto is = sub_data->getSmartColIs();
3015 auto dm =
simple->getDM();
3023 CHKERR TSGetSNES(ts, &snes);
3026 auto set_section_monitor = [&](
auto snes) {
3028 CHKERR SNESMonitorSet(snes,
3031 (
void *)(snes_ctx_ptr.get()),
nullptr);
3035 CHKERR set_section_monitor(snes);
3037 auto ksp =
createKSP(m_field.get_comm());
3038 CHKERR KSPSetType(ksp, KSPPREONLY);
3040 CHKERR PCSetType(sub_pc, PCSHELL);
3043 CHKERR KSPSetPC(ksp, sub_pc);
3044 CHKERR SNESSetKSP(snes, ksp);
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_DATA_INCONSISTENCY
CoordinateTypes
Coodinate system.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
constexpr CoordinateTypes coord_type
FTensor::Index< 'k', SPACE_DIM > k
constexpr double lambda
surface tension
OpDomainMassH OpDomainMassG
FTensor::Index< 'i', SPACE_DIM > i
auto init_h
Initialisation function.
FTensor::Index< 'l', SPACE_DIM > l
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, coord_type > OpMixScalarTimesDiv
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
auto get_dofs_ents_all
get entities of dofs in the problem - used for debugging
constexpr IntegrationType I
auto get_skin_projection_bit
auto get_dofs_ents_by_field_name
get entities of dofs in the problem - used for debugging
ElementsAndOps< SPACE_DIM >::DomianParentEle DomianParentEle
constexpr double rho_diff
auto wetting_angle_sub_stepping
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
int order
approximation order
SideEle::UserDataOperator SideOp
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
auto get_current_bit
dofs bit used to do calculations
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomianParentEle DomianParentEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
auto createTS(MPI_Comm comm)
auto createPC(MPI_Comm comm)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscObject getPetscObject(T obj)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode runProblem()
[Run programme]
FreeSurface(MoFEM::Interface &m_field)
MoFEMErrorCode projectData()
[Boundary condition]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parents levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
[Data projection]
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
MoFEMErrorCode makeRefProblem()
MoFEMErrorCode solveSystem()
[Solve]
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak ptr object.
Simple interface for fast problem set-up.
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
default operator for TRI element
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
@ OPSPACE
operator do Work is execute on space data
Section manager is used to create indexes and sections.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Interface for managing meshsets containing materials and boundary conditions.
Operator to project base functions from parent entity to child.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< moab::Core > postProcMesh
MoFEMErrorCode postProcess()
function is run at the end of loop
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< moab::Core > post_proc_mesh, boost::shared_ptr< PostProcEleDomainCont > post_proc, boost::shared_ptr< PostProcEleBdyCont > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
boost::shared_ptr< PostProcEleDomainCont > postProc
Set of functions called by PETSc solver used to refine and update mesh.
boost::shared_ptr< TsCtx > tsCtxPtr
virtual ~TSPrePostProc()=default
SmartPetscObj< Vec > globRes
SmartPetscObj< Mat > subB
SmartPetscObj< Vec > globSol
SmartPetscObj< Vec > getSubVector()
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Wrapper for TS monitor.
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
SmartPetscObj< DM > solverSubDM
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Wrapper for SNES Lhs.
boost::shared_ptr< SnesCtx > snesCtxPtr
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
static MoFEMErrorCode pcSetup(PC pc)
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
SmartPetscObj< KSP > subKSP