800 auto dm =
simple->getDM();
801 auto solver = pipeline_mng->createTSIM();
804 auto set_section_monitor = [&](
auto solver) {
807 CHKERR TSGetSNES(solver, &snes);
808 CHKERR SNESMonitorSet(snes,
811 (
void *)(snes_ctx_ptr.get()),
nullptr);
815 auto create_post_process_element = [&]() {
816 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
818 auto block_params = boost::make_shared<BlockedParameters>();
819 auto mDPtr = block_params->getDPtr();
821 "MAT_FLUID", block_params, Sev::verbose);
823 post_proc_fe->getOpPtrVector(), {H1, HDIV});
825 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
826 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
827 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
829 auto h_ptr = boost::make_shared<VectorDouble>();
830 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
832 post_proc_fe->getOpPtrVector().push_back(
834 post_proc_fe->getOpPtrVector().push_back(
837 auto u_ptr = boost::make_shared<MatrixDouble>();
838 post_proc_fe->getOpPtrVector().push_back(
840 post_proc_fe->getOpPtrVector().push_back(
843 post_proc_fe->getOpPtrVector().push_back(
845 post_proc_fe->getOpPtrVector().push_back(
847 "U", mat_strain_ptr, mat_stress_ptr, mDPtr));
851 post_proc_fe->getOpPtrVector().push_back(
855 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
859 {{
"U", u_ptr}, {
"FLUX", mat_flux_ptr}},
863 {{
"STRAIN", mat_strain_ptr}, {
"STRESS", mat_stress_ptr}}
872 auto create_creaction_fe = [&]() {
873 auto fe_ptr = boost::make_shared<DomainEle>(
mField);
874 fe_ptr->getRuleHook = [](
int,
int,
int o) {
return 2 * o; };
876 auto &pip = fe_ptr->getOpPtrVector();
878 auto block_params = boost::make_shared<BlockedParameters>();
879 auto mDPtr = block_params->getDPtr();
884 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
885 auto strain_ptr = boost::make_shared<MatrixDouble>();
886 auto stress_ptr = boost::make_shared<MatrixDouble>();
894 "U", strain_ptr, stress_ptr, mDPtr));
897 fe_ptr->postProcessHook =
903 auto monitor_ptr = boost::make_shared<FEMethod>();
904 auto post_proc_fe = create_post_process_element();
906 auto rections_fe = create_creaction_fe();
908 auto set_time_monitor = [&](
auto dm,
auto solver) {
910 monitor_ptr->preProcessHook = [&]() {
915 monitor_ptr->getCacheWeakPtr());
916 CHKERR post_proc_fe->writeFile(
917 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
920 rections_fe->f = res;
923 monitor_ptr->getCacheWeakPtr());
926 CHKERR VecNorm(res, NORM_2, &nrm);
928 <<
"Residual norm " << nrm <<
" at time step "
929 << monitor_ptr->ts_step;
933 auto scalar_field_ptr = boost::make_shared<VectorDouble>();
934 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
935 auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
938 ->getData<DomainEle>();
942 field_eval_data,
simple->getDomainFEName());
945 field_eval_data,
simple->getDomainFEName());
949 auto no_rule = [](
int,
int,
int) {
return -1; };
951 auto field_eval_ptr = field_eval_data->feMethodPtr.lock();
952 field_eval_ptr->getRuleHook = no_rule;
954 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
955 auto strain_ptr = boost::make_shared<MatrixDouble>();
956 auto stress_ptr = boost::make_shared<MatrixDouble>();
957 auto h_ptr = boost::make_shared<VectorDouble>();
959 auto block_params = boost::make_shared<BlockedParameters>();
960 auto mDPtr = block_params->getDPtr();
962 "MAT_FLUID", block_params, Sev::noisy);
964 field_eval_ptr->getOpPtrVector(), {H1, HDIV});
965 field_eval_ptr->getOpPtrVector().push_back(
968 field_eval_ptr->getOpPtrVector().push_back(
970 field_eval_ptr->getOpPtrVector().push_back(
972 field_eval_ptr->getOpPtrVector().push_back(
974 "U", strain_ptr, stress_ptr, mDPtr));
978 ->evalFEAtThePoint3D(
980 simple->getDomainFEName(), field_eval_data,
985 ->evalFEAtThePoint2D(
987 simple->getDomainFEName(), field_eval_data,
993 <<
"Eval point hydrostatic hight: " << *h_ptr;
995 <<
"Eval point skeleton stress pressure: " << *stress_ptr;
1001 auto null = boost::shared_ptr<FEMethod>();
1007 auto set_fieldsplit_preconditioner = [&](
auto solver) {
1011 CHKERR TSGetSNES(solver, &snes);
1013 CHKERR SNESGetKSP(snes, &ksp);
1015 CHKERR KSPGetPC(ksp, &pc);
1016 PetscBool is_pcfs = PETSC_FALSE;
1017 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1020 if (is_pcfs == PETSC_TRUE) {
1023 auto name_prb =
simple->getProblemName();
1026 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
1029 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
1032 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"H", 0, 0,
1035 CHKERR ISExpand(is_H, is_flux, &is_tmp);
1038 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"FLUIDFLUX",
"FLUX", 0, 0);
1040 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1042 <<
"Field split block size " << is_all_bc_size;
1043 if (is_all_bc_size) {
1045 CHKERR ISDifference(is_Flux, is_all_bc, &is_tmp2);
1047 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1053 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_Flux);
1054 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
1060 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1061 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1062 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1063 auto time_scale = boost::make_shared<TimeScale>();
1065 auto get_bc_hook_rhs = [
this, pre_proc_ptr, time_scale]() {
1067 {time_scale},
false);
1071 auto get_post_proc_hook_rhs = [
this, post_proc_rhs_ptr]() {
1074 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1076 mField, post_proc_rhs_ptr, 1.)();
1079 auto get_post_proc_hook_lhs = [
this, post_proc_lhs_ptr]() {
1081 post_proc_lhs_ptr, 1.);
1084 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1085 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1086 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1089 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1090 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1091 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1092 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1095 CHKERR TSSetSolution(solver,
D);
1096 CHKERR TSSetFromOptions(solver);
1098 CHKERR set_section_monitor(solver);
1099 CHKERR set_fieldsplit_preconditioner(solver);
1100 CHKERR set_time_monitor(dm, solver);
1103 CHKERR TSSolve(solver, NULL);