Project solution data between mesh levels.
Projects field data from coarse to fine mesh levels during refinement. Handles both time stepping vectors (for theta method) and regular fields. Uses L2 projection with mass matrix assembly.
1388 {
1390
1391
1392
1393
1394
1399
1400
1401
1402 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
1403 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
1404 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
1405 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
1406
1407 pip_mng->getDomainLhsFE().reset();
1408 pip_mng->getDomainRhsFE().reset();
1409 pip_mng->getBoundaryLhsFE().reset();
1410 pip_mng->getBoundaryRhsFE().reset();
1411
1412 using UDO = ForcesAndSourcesCore::UserDataOperator;
1413
1414
1415 auto add_parent_field_domain = [&](
auto fe,
auto op,
auto field,
auto bit) {
1418
1419 [&]() {
1420 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1422 return fe_parent;
1423 },
1424
1426 };
1427
1428
1429 auto add_parent_field_bdy = [&](
auto fe,
auto op,
auto field,
auto bit) {
1432
1433 [&]() {
1434 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1436 return fe_parent;
1437 },
1438
1440 };
1441
1442
1443
1444 auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
1445 auto fe_bit) {
1446 auto parent_mask = fe_bit;
1447 parent_mask.flip();
1450 Sev::inform);
1451 };
1452
1453
1454 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
1455 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
1457 parents_elems_ptr_vec.emplace_back(
1458 boost::make_shared<DomainParentEle>(
mField));
1460 parents_elems_ptr_vec[
l - 1]->getOpPtrVector().push_back(
1461 create_run_parent_op(parents_elems_ptr_vec[
l], this_fe_ptr, fe_bit));
1462 }
1463 return parents_elems_ptr_vec[0];
1464 };
1465
1466
1467 auto solve_projection = [&](auto exe_test) {
1469
1470 auto set_domain_rhs = [&](auto fe) {
1472 auto &pip = fe->getOpPtrVector();
1473
1474 auto u_ptr = boost::make_shared<MatrixDouble>();
1475 auto p_ptr = boost::make_shared<VectorDouble>();
1476 auto h_ptr = boost::make_shared<VectorDouble>();
1477 auto g_ptr = boost::make_shared<VectorDouble>();
1478
1479 auto eval_fe_ptr = boost::make_shared<DomainParentEle>(
mField);
1480 {
1481 auto &pip = eval_fe_ptr->getOpPtrVector();
1485
1486 [&]() {
1487 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1489 return fe_parent;
1490 },
1491
1493
1494
1495 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"U",
1498 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"P",
1501 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"H",
1504 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW,
"G",
1507 }
1508 auto parent_eval_fe_ptr =
1510 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1512
1513 auto assemble_fe_ptr = boost::make_shared<DomainParentEle>(
mField);
1514 {
1515 auto &pip = assemble_fe_ptr->getOpPtrVector();
1519
1520 [&]() {
1521 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1523 return fe_parent;
1524 },
1525
1527 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"U",
1530 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"P",
1533 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"H",
1536 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW,
"G",
1539 }
1540 auto parent_assemble_fe_ptr =
1542 pip.push_back(create_run_parent_op(
1544
1546 };
1547
1548 auto set_domain_lhs = [&](auto fe) {
1550
1551 auto &pip = fe->getOpPtrVector();
1552
1554
1557
1558 [&]() {
1559 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1561 return fe_parent;
1562 },
1563
1565
1566
1567
1568 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"U",
1570 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"U",
1573 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"P",
1575 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"P",
1578 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"H",
1580 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"H",
1583 CHKERR add_parent_field_domain(fe, UDO::OPROW,
"G",
1585 CHKERR add_parent_field_domain(fe, UDO::OPCOL,
"G",
1588
1590 };
1591
1592 auto set_bdy_rhs = [&](auto fe) {
1594 auto &pip = fe->getOpPtrVector();
1595
1596 auto l_ptr = boost::make_shared<VectorDouble>();
1597
1598 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1599 {
1600 auto &pip = eval_fe_ptr->getOpPtrVector();
1603
1604 [&]() {
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1607 return fe_parent;
1608 },
1609
1611
1612
1613 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW,
"L",
1616 }
1617 auto parent_eval_fe_ptr =
1619 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1621
1622 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(
mField);
1623 {
1624 auto &pip = assemble_fe_ptr->getOpPtrVector();
1627
1628 [&]() {
1629 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1631 return fe_parent;
1632 },
1633
1635
1637 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1641 if (lPtr->size() != getGaussPts().size2()) {
1642 lPtr->resize(getGaussPts().size2());
1643 lPtr->clear();
1644 }
1646 }
1647
1648 private:
1649 boost::shared_ptr<VectorDouble> lPtr;
1650 };
1651
1652 pip.push_back(new OpLSize(l_ptr));
1653
1654 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW,
"L",
1657 }
1658 auto parent_assemble_fe_ptr =
1660 pip.push_back(create_run_parent_op(
1662
1664 };
1665
1666 auto set_bdy_lhs = [&](auto fe) {
1668
1669 auto &pip = fe->getOpPtrVector();
1670
1673
1674 [&]() {
1675 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1677 return fe_parent;
1678 },
1679
1681
1682
1683
1684 CHKERR add_parent_field_bdy(fe, UDO::OPROW,
"L",
1686 CHKERR add_parent_field_bdy(fe, UDO::OPCOL,
"L",
1689
1691 };
1692
1694 auto level_ents_ptr = boost::make_shared<Range>();
1697
1698 auto get_prj_ents = [&]() {
1703 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1704 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1706
1707 return prj_mesh;
1708 };
1709
1710 auto prj_ents = get_prj_ents();
1711
1713
1714 auto rebuild = [&]() {
1717
1718 std::vector<std::string> fields{"U", "P", "H", "G", "L"};
1719 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1720
1721 {"U", level_ents_ptr},
1722 {"P", level_ents_ptr},
1723 {"H", level_ents_ptr},
1724 {"G", level_ents_ptr},
1725 {"L", level_ents_ptr}
1726
1727 };
1728
1730 simple->getProblemName(), PETSC_TRUE,
1731 &range_maps, &range_maps);
1732
1733
1734 CHKERR prb_mng->partitionFiniteElements(
"SUB_SOLVER",
true, 0,
1736
1737 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh(
"SUB_SOLVER");
1738
1740 };
1741
1742 MOFEM_LOG(
"FS", Sev::verbose) <<
"Create projection problem";
1743
1745
1746 auto dm =
simple->getDM();
1747 DM subdm;
1749 CHKERR DMSetType(subdm,
"DMMOFEM");
1752 }
1753
1754 MOFEM_LOG(
"FS", Sev::inform) <<
"Nothing to project";
1755
1757 };
1758
1759 auto create_dummy_dm = [&]() {
1762 simple->getProblemName().c_str(),
1764 "create dummy dm");
1765 return dummy_dm;
1766 };
1767
1768 auto subdm = create_subdm();
1769 if (subdm) {
1770
1771 pip_mng->getDomainRhsFE().reset();
1772 pip_mng->getDomainLhsFE().reset();
1773 pip_mng->getBoundaryRhsFE().reset();
1774 pip_mng->getBoundaryLhsFE().reset();
1779 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1780 pip_mng->getDomainRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1782 };
1783 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1784 pip_mng->getBoundaryRhsFE()->exeTestHook = [](
FEMethod *fe_ptr) {
1786 };
1787
1788 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1789 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1790 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1791 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1792
1795
1796 auto ksp = pip_mng->createKSP(subdm);
1797 CHKERR KSPSetFromOptions(ksp);
1799
1800
1801 auto get_norm = [&](auto x) {
1802 double nrm;
1803 CHKERR VecNorm(x, NORM_2, &nrm);
1804 return nrm;
1805 };
1806
1807
1808
1809
1810 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1812 for (
auto &
v : ent_ptr->getEntFieldData()) {
1814 }
1816 };
1817
1818 auto solve = [&](
auto S) {
1820 CHKERR VecZeroEntries(S);
1822 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
1823 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
1825 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1826 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1827
1828
1829
1831 };
1832
1833 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection";
1835
1838 auto dummy_dm = create_dummy_dm();
1839
1840
1841
1842
1843 auto apply_restrict = [&]() {
1844 auto get_is = [](
auto v) {
1845 IS iy;
1846 auto create = [&]() {
1850 CHKERR VecGetOwnershipRange(
v, &ystart, NULL);
1851 CHKERR ISCreateStride(PETSC_COMM_SELF,
n, ystart, 1, &iy);
1853 };
1856 };
1857
1858 auto iy = get_is(glob_x);
1860
1862 DMSubDomainRestrict(
simple->getDM(), s, PETSC_NULLPTR, dummy_dm),
1863 "restrict");
1866 "get X0");
1868 DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1869 "get Xdot");
1870
1871 auto forward_ghost = [](
auto g) {
1873 CHKERR VecGhostUpdateBegin(
g, INSERT_VALUES, SCATTER_FORWARD);
1874 CHKERR VecGhostUpdateEnd(
g, INSERT_VALUES, SCATTER_FORWARD);
1876 };
1877
1880
1881 if constexpr (
debug) {
1883 << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1884 << get_norm(Xdot);
1885 }
1886
1887 return std::vector<Vec>{X0, Xdot};
1888 };
1889
1890 auto ts_solver_vecs = apply_restrict();
1891
1892 if (ts_solver_vecs.size()) {
1893
1894 for (
auto v : ts_solver_vecs) {
1895 MOFEM_LOG(
"FS", Sev::inform) <<
"Solve projection vector";
1896
1898 SCATTER_REVERSE);
1900
1901 for (auto f : {"U", "P", "H", "G", "L"}) {
1902 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1903 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1904 }
1905
1907 SCATTER_REVERSE);
1909 SCATTER_FORWARD);
1910
1911 MOFEM_LOG(
"FS", Sev::inform) <<
"Norm V " << get_norm(
v);
1912 }
1913
1914 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_X0",
1915 &ts_solver_vecs[0]);
1916 CHKERR DMRestoreNamedGlobalVector(dummy_dm,
"TSTheta_Xdot",
1917 &ts_solver_vecs[1]);
1918 }
1919
1920 for (auto f : {"U", "P", "H", "G", "L"}) {
1921 MOFEM_LOG(
"WORLD", Sev::verbose) <<
"Zero field " <<
f;
1922 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1923 }
1925 }
1926
1928 };
1929
1930
1931 auto post_proc = [&](auto exe_test) {
1933 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(
mField);
1934 auto &pip = post_proc_fe->getOpPtrVector();
1935 post_proc_fe->exeTestHook = exe_test;
1936
1938
1941
1942 [&]() {
1943 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1946 fe_parent->getOpPtrVector(), {H1});
1947 return fe_parent;
1948 },
1949
1951
1953
1954 auto u_ptr = boost::make_shared<MatrixDouble>();
1955 auto p_ptr = boost::make_shared<VectorDouble>();
1956 auto h_ptr = boost::make_shared<VectorDouble>();
1957 auto g_ptr = boost::make_shared<VectorDouble>();
1958
1959 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"U",
1962 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"P",
1965 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"H",
1968 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW,
"G",
1971
1972 post_proc_fe->getOpPtrVector().push_back(
1973
1974 new OpPPMap(post_proc_fe->getPostProcMesh(),
1975 post_proc_fe->getMapGaussPts(),
1976
1977 {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1978
1979 {{"U", u_ptr}},
1980
1981 {},
1982
1983 {}
1984
1985 )
1986
1987 );
1988
1989 auto dm =
simple->getDM();
1991 CHKERR post_proc_fe->writeFile(
"out_projection.h5m");
1992
1994 };
1995
1998 });
1999
2000 if constexpr (
debug) {
2003 });
2004 }
2005
2006 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
2007 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
2008 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
2009 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
2010
2012}
OpDomainMassH OpDomainMassG
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
auto get_skin_projection_bit
auto get_global_size
Get global size across all processors.
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
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.
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
const double n
refractive index of diffusive medium
const FTensor::Tensor2< T, Dim, Dim > Vec
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.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Data on single entity (This is passed as argument to DataOperator::doWork)
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Problem manager is used to build and partition problems.