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) {
985 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
994 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
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()) {
1460 for (
auto v : ts_solver_vecs) {
1461 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1467 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1468 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1469 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
1477 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1480 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1481 &ts_solver_vecs[0]);
1482 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1483 &ts_solver_vecs[1]);
1486 for (
auto f : {
"U",
"P",
"H",
"G",
"L"}) {
1487 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1488 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs,
f);
1497 auto post_proc = [&](
auto exe_test) {
1499 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1500 auto &pip = post_proc_fe->getOpPtrVector();
1501 post_proc_fe->exeTestHook = exe_test;
1509 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1512 fe_parent->getOpPtrVector(), {H1});
1520 auto u_ptr = boost::make_shared<MatrixDouble>();
1521 auto p_ptr = boost::make_shared<VectorDouble>();
1522 auto h_ptr = boost::make_shared<VectorDouble>();
1523 auto g_ptr = boost::make_shared<VectorDouble>();
1525 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1528 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1531 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1534 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1538 post_proc_fe->getOpPtrVector().push_back(
1540 new OpPPMap(post_proc_fe->getPostProcMesh(),
1541 post_proc_fe->getMapGaussPts(),
1543 {{
"P", p_ptr}, {
"H", h_ptr}, {
"G", g_ptr}},
1555 auto dm =
simple->getDM();
1557 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1566 if constexpr (
debug) {
1572 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
1573 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
1574 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
1575 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());