1290 {
1292
1296
1297 auto dm =
simple->getDM();
1298 auto solver = pipeline_mng->createTSIM();
1300
1301 auto set_section_monitor = [&](auto solver) {
1303 SNES snes;
1304 CHKERR TSGetSNES(solver, &snes);
1305 CHKERR SNESMonitorSet(snes,
1308 (void *)(snes_ctx_ptr.get()), nullptr);
1310 };
1311
1312 auto create_post_process_elements = [&]() {
1313 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
1314
1315 auto block_thermoelastic_params =
1316 boost::make_shared<BlockedThermoElasticParameters>();
1317 auto coeff_expansion_ptr =
1318 block_thermoelastic_params->getCoeffExpansionPtr();
1319 auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
1320
1321 auto u_ptr = boost::make_shared<MatrixDouble>();
1322 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1323 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
1324 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
1325 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1326 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1327
1328 auto push_domain_ops = [&](auto &pp_fe) {
1330 auto &pip = pp_fe->getOpPtrVector();
1331
1333 Sev::verbose);
1335 pip, "MAT_THERMOELASTIC", block_thermoelastic_params, Sev::verbose);
1336
1338
1340 pip.push_back(
1342
1345 "U", mat_grad_ptr));
1346 auto elastic_common_ptr = commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
1347 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
1349 pip.push_back(
1350 new
1351 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
1352 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
1353 ref_temp_ptr));
1356 elastic_common_ptr->getMatFirstPiolaStress(), mat_stress_ptr));
1357 pip.push_back(
1359 } else {
1360 mat_stress_ptr = elastic_common_ptr->getMatFirstPiolaStress();
1361 mat_strain_ptr = elastic_common_ptr->getMatLogC();
1362 }
1363
1365 };
1366
1367 auto push_post_proc_ops = [&](auto &pp_fe) {
1369 auto &pip = pp_fe->getOpPtrVector();
1371
1373 pip.push_back(
1374
1376
1377 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1378
1379 {{"T", vec_temp_ptr}},
1380
1381 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1382
1383 {},
1384
1385 {{"CAUCHY", mat_stress_ptr}, {"STRAIN", mat_strain_ptr}}
1386
1387 )
1388
1389 );
1390 } else {
1391 pip.push_back(
1392
1394
1395 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1396
1397 {{"T", vec_temp_ptr}},
1398
1399 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1400
1401 {{"PIOLA", mat_stress_ptr}},
1402
1403 {{"HENCKY_STRAIN", mat_strain_ptr}}
1404
1405 )
1406
1407 );
1408 }
1409
1411 };
1412
1413 auto domain_post_proc = [&]() {
1415 return boost::shared_ptr<PostProcEle>();
1416 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
1418 "push domain ops to domain element");
1420 "push post proc ops to domain element");
1421 return pp_fe;
1422 };
1423
1424 auto skin_post_proc = [&]() {
1426 return boost::shared_ptr<SkinPostProcEle>();
1427 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
1432 "push domain ops to side element");
1433 pp_fe->getOpPtrVector().push_back(op_side);
1435 "push post proc ops to skin element");
1436 return pp_fe;
1437 };
1438
1439 return std::make_pair(domain_post_proc(), skin_post_proc());
1440 };
1441
1442 auto monitor_ptr = boost::make_shared<FEMethod>();
1443
1444 auto set_time_monitor = [&](auto dm, auto solver, auto domain_post_proc_fe,
1445 auto skin_post_proc_fe) {
1447 monitor_ptr->preProcessHook = [&]() {
1449
1453 domain_post_proc_fe,
1454 monitor_ptr->getCacheWeakPtr());
1455 CHKERR domain_post_proc_fe->writeFile(
1456 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
1457 ".h5m");
1458 }
1461 skin_post_proc_fe,
1462 monitor_ptr->getCacheWeakPtr());
1463 CHKERR skin_post_proc_fe->writeFile(
1464 "out_skin_" +
1465 boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
1466 }
1467 }
1468
1469 struct AtomTestResult {
1470 bool fail = false;
1471 std::string msg = "";
1472 };
1473 AtomTestResult atom_test_result;
1474
1475 auto fail_atom_test = [&atom_test_result](const std::string &msg) {
1476 atom_test_result.fail = true;
1477 atom_test_result.msg = msg;
1478 };
1479
1480 struct AtomTestData{
1481 double expected = 0.0;
1483 };
1484
1486
1488 ->evalFEAtThePoint<SPACE_DIM>(
1493
1495 auto eval_num_vec =
1497 CHKERR VecZeroEntries(eval_num_vec);
1499 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
1500 }
1501 CHKERR VecAssemblyBegin(eval_num_vec);
1502 CHKERR VecAssemblyEnd(eval_num_vec);
1503
1504 double eval_num;
1505 CHKERR VecSum(eval_num_vec, &eval_num);
1506 if (!(int)eval_num) {
1507 fail_atom_test(
1508 "did not find elements to evaluate the field, check the "
1509 "coordinates");
1510 }
1511 }
1512
1515 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1516 << "Eval point T: " << t_temp;
1517 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1518 AtomTestData atom_test_data;
1520 case 1:
1521 case 2:
1522 case 3:
1523 atom_test_data = {554.48, 1e-2};
1524 break;
1525 case 4:
1526 atom_test_data = {250.0, 1e-2};
1527 break;
1528 case 5:
1529 atom_test_data = {1.0, 1e-2};
1530 break;
1531 default:
1532 fail_atom_test("unknown atom test number");
1533 }
1534 if (fabs(t_temp - atom_test_data.expected) > atom_test_data.tol) {
1535 fail_atom_test("wrong temperature value");
1536 }
1537 }
1538 }
1541 auto t_flux = getFTensor1FromMat<SPACE_DIM>(*
fluxFieldPtr);
1542 auto flux_mag = sqrt(t_flux(
i) * t_flux(
i));
1543 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1544 << "Eval point FLUX magnitude: " << flux_mag;
1545 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1546 AtomTestData atom_test_data;
1548 case 1:
1549 case 2:
1550 case 3:
1551 atom_test_data = {27008.0, 2e1};
1552 break;
1553 case 4:
1554 atom_test_data = {150e3, 1e-1};
1555 break;
1556 case 5:
1557 atom_test_data = {0.0, 1e-6};
1558 break;
1559 default:
1560 fail_atom_test("unknown atom test number");
1561 }
1562 if (fabs(flux_mag - atom_test_data.expected) > atom_test_data.tol) {
1563 fail_atom_test("wrong flux value");
1564 }
1565 }
1566 }
1569 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*
dispFieldPtr);
1570 auto disp_mag = sqrt(t_disp(
i) * t_disp(
i));
1571 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1572 << "Eval point U magnitude: " << disp_mag;
1573 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1574 AtomTestData atom_test_data;
1576 case 1:
1577 atom_test_data = {0.00345, 1e-5};
1578 break;
1579 case 2:
1580 case 3:
1581 atom_test_data = {0.00265, 1e-5};
1582 break;
1583 case 4:
1584 atom_test_data = {0.00061, 1e-5};
1585 break;
1586 case 5:
1587 atom_test_data = {std::sqrt(2) * (std::sqrt(std::exp(0.2)) - 1),
1588 1e-5};
1589 break;
1590 default:
1591 fail_atom_test("unknown atom test number");
1592 }
1593 if (fabs(disp_mag - atom_test_data.expected) > atom_test_data.tol) {
1594 fail_atom_test("wrong displacement value");
1595 }
1596 }
1597 }
1600 auto t_strain =
1602 auto t_strain_trace = t_strain(
i,
i);
1603 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1604 AtomTestData atom_test_data;
1606 case 1:
1607 atom_test_data = {0.00679, 1e-5};
1608 break;
1609 case 2:
1610 case 3:
1611 atom_test_data = {0.00522, 1e-5};
1612 break;
1613 case 4:
1614 atom_test_data = {-0.00085, 1e-5};
1615 break;
1616 case 5:
1617 atom_test_data = {0.2, 1e-5};
1618 break;
1619 default:
1620 fail_atom_test("unknown atom test number");
1621 }
1622 if (fabs(t_strain_trace - atom_test_data.expected) >
1623 atom_test_data.tol) {
1624 fail_atom_test("wrong strain value");
1625 }
1626 }
1627 }
1629 double von_mises_stress;
1630 auto getVonMisesStress = [&](auto t_stress) {
1632 von_mises_stress = std::sqrt(
1633 0.5 *
1634 ((t_stress(0, 0) - t_stress(1, 1)) *
1635 (t_stress(0, 0) - t_stress(1, 1)) +
1636 (
SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
1637 (t_stress(1, 1) - t_stress(2, 2))
1638 : 0) +
1639 (
SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
1640 (t_stress(2, 2) - t_stress(0, 0))
1641 : 0) +
1642 6 * (t_stress(0, 1) * t_stress(0, 1) +
1643 (
SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
1644 (
SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
1646 };
1648 auto t_stress =
1650 CHKERR getVonMisesStress(t_stress);
1651 } else {
1652 auto t_stress =
1654 CHKERR getVonMisesStress(t_stress);
1655 }
1656 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1657 << "Eval point von Mises Stress: " << von_mises_stress;
1658 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1659 AtomTestData atom_test_data;
1661 case 1:
1662 atom_test_data = {523.0, 5e-1};
1663 break;
1664 case 2:
1665 atom_test_data = {16.3, 5e-2};
1666 break;
1667 case 3:
1668 atom_test_data = {14.9, 5e-2};
1669 break;
1670 case 4:
1671 atom_test_data = {92.3, 2e-1};
1672 break;
1673 case 5:
1674 atom_test_data = {0.0, 5e-2};
1675 break;
1676 default:
1677 fail_atom_test("unknown atom test number");
1678 }
1679 if (fabs(von_mises_stress - atom_test_data.expected) >
1680 atom_test_data.tol) {
1681 fail_atom_test("wrong von Mises stress value");
1682 }
1683 }
1684 }
1686 }
1687
1688 if (atom_test_result.fail) {
1691 atom_test_result.msg.c_str());
1692 }
1693
1695 };
1696 auto null = boost::shared_ptr<FEMethod>();
1698 monitor_ptr, null);
1700 };
1701
1702 auto set_fieldsplit_preconditioner = [&](auto solver) {
1704
1705 SNES snes;
1706 CHKERR TSGetSNES(solver, &snes);
1707 KSP ksp;
1708 CHKERR SNESGetKSP(snes, &ksp);
1709 PC pc;
1710 CHKERR KSPGetPC(ksp, &pc);
1711 PetscBool is_pcfs = PETSC_FALSE;
1712 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1713
1714
1715 if (is_pcfs == PETSC_TRUE) {
1717 auto name_prb =
simple->getProblemName();
1718
1720 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
1723 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
1724 is_flux);
1726 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"T", 0, 0,
1727 is_T);
1729 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"TBC", 0, 0,
1730 is_TBC);
1731 IS is_tmp, is_tmp2;
1732 CHKERR ISExpand(is_T, is_flux, &is_tmp);
1733 CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
1734 CHKERR ISDestroy(&is_tmp);
1736
1739 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_TFlux);
1740 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1741 }
1742
1744 };
1745
1747 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1750 CHKERR TSSetSolution(solver,
D);
1751 CHKERR TSSetFromOptions(solver);
1752
1753 CHKERR set_section_monitor(solver);
1754 CHKERR set_fieldsplit_preconditioner(solver);
1755
1756 auto [domain_post_proc_fe, skin_post_proc_fe] =
1757 create_post_process_elements();
1758 CHKERR set_time_monitor(dm, solver, domain_post_proc_fe, skin_post_proc_fe);
1759
1761 CHKERR TSSolve(solver, NULL);
1762
1764}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
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.
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Section manager is used to create indexes and sections.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
intrusive_ptr for managing petsc objects
