792 {
796
797 auto dm =
simple->getDM();
800
801 auto set_section_monitor = [&](auto solver) {
803 SNES snes;
804 CHKERR TSGetSNES(solver, &snes);
805 CHKERR SNESMonitorSet(snes,
808 (void *)(snes_ctx_ptr.get()), nullptr);
810 };
811
812 auto create_post_process_element = [&]() {
813 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
814
815 auto block_params = boost::make_shared<BlockedParameters>();
816 auto mDPtr = block_params->getDPtr();
818 "MAT_FLUID", block_params, Sev::verbose);
820 post_proc_fe->getOpPtrVector(), {H1, HDIV});
821
822 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
823 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
824 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
825
826 auto h_ptr = boost::make_shared<VectorDouble>();
827 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
828
829 post_proc_fe->getOpPtrVector().push_back(
831 post_proc_fe->getOpPtrVector().push_back(
833
834 auto u_ptr = boost::make_shared<MatrixDouble>();
835 post_proc_fe->getOpPtrVector().push_back(
837 post_proc_fe->getOpPtrVector().push_back(
839 mat_grad_ptr));
840 post_proc_fe->getOpPtrVector().push_back(
842 post_proc_fe->getOpPtrVector().push_back(
844 mat_strain_ptr, mat_stress_ptr, mDPtr));
845
847
848 post_proc_fe->getOpPtrVector().push_back(
849
851
852 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
853
854 {{"H", h_ptr}},
855
856 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
857
858 {},
859
860 {{"STRAIN", mat_strain_ptr}, {"STRESS", mat_stress_ptr}}
861
862 )
863
864 );
865
866 return post_proc_fe;
867 };
868
869 auto create_creaction_fe = [&]() {
870 auto fe_ptr = boost::make_shared<DomainEle>(
mField);
871 fe_ptr->getRuleHook = [](int, int, int o) { return 2 * o; };
872
873 auto &pip = fe_ptr->getOpPtrVector();
874
875 auto block_params = boost::make_shared<BlockedParameters>();
876 auto mDPtr = block_params->getDPtr();
878 Sev::verbose);
880
881 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
882 auto strain_ptr = boost::make_shared<MatrixDouble>();
883 auto stress_ptr = boost::make_shared<MatrixDouble>();
884
886 "U", u_grad_ptr));
888
889
891 strain_ptr, stress_ptr, mDPtr));
893
894 fe_ptr->postProcessHook =
896
897 return fe_ptr;
898 };
899
900 auto monitor_ptr = boost::make_shared<FEMethod>();
901 auto post_proc_fe = create_post_process_element();
903 auto rections_fe = create_creaction_fe();
904
905 auto set_time_monitor = [&](auto dm, auto solver) {
907 monitor_ptr->preProcessHook = [&]() {
909
911 post_proc_fe,
912 monitor_ptr->getCacheWeakPtr());
913 CHKERR post_proc_fe->writeFile(
914 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
915 ".h5m");
916
917 rections_fe->f = res;
919 rections_fe,
920 monitor_ptr->getCacheWeakPtr());
921
922 CHKERR VecAssemblyBegin(res);
923 CHKERR VecAssemblyEnd(res);
924 double nrm;
925 CHKERR VecNorm(res, NORM_2, &nrm);
927 << "Residual norm " << nrm << " at time step "
928 << monitor_ptr->ts_step;
929
931
932 auto scalar_field_ptr = boost::make_shared<VectorDouble>();
933 auto vector_field_ptr = boost::make_shared<MatrixDouble>();
934 auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
935
937 ->getData<DomainEle>();
938
940 ->buildTree<SPACE_DIM>(field_eval_data,
simple->getDomainFEName());
941
943 auto no_rule = [](int, int, int) { return -1; };
944
945 auto field_eval_ptr = field_eval_data->feMethodPtr.lock();
946 field_eval_ptr->getRuleHook = no_rule;
947
948 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
949 auto strain_ptr = boost::make_shared<MatrixDouble>();
950 auto stress_ptr = boost::make_shared<MatrixDouble>();
951 auto h_ptr = boost::make_shared<VectorDouble>();
952
953 auto block_params = boost::make_shared<BlockedParameters>();
954 auto mDPtr = block_params->getDPtr();
956 "MAT_FLUID", block_params, Sev::noisy);
958 field_eval_ptr->getOpPtrVector(), {H1, HDIV});
959 field_eval_ptr->getOpPtrVector().push_back(
961 "U", u_grad_ptr));
962 field_eval_ptr->getOpPtrVector().push_back(
964 field_eval_ptr->getOpPtrVector().push_back(
966 field_eval_ptr->getOpPtrVector().push_back(
968 strain_ptr, stress_ptr, mDPtr));
969
971 ->evalFEAtThePoint<SPACE_DIM>(
973 simple->getDomainFEName(), field_eval_data,
976
978 << "Eval point hydrostatic hight: " << *h_ptr;
980 << "Eval point skeleton stress pressure: " << *stress_ptr;
982 }
983
985 };
986 auto null = boost::shared_ptr<FEMethod>();
988 monitor_ptr, null);
990 };
991
992 auto set_fieldsplit_preconditioner = [&](auto solver) {
994
995 SNES snes;
996 CHKERR TSGetSNES(solver, &snes);
997 KSP ksp;
998 CHKERR SNESGetKSP(snes, &ksp);
999 PC pc;
1000 CHKERR KSPGetPC(ksp, &pc);
1001 PetscBool is_pcfs = PETSC_FALSE;
1002 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1003
1004
1005 if (is_pcfs == PETSC_TRUE) {
1008 auto name_prb =
simple->getProblemName();
1009
1011 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
1014 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
1015 is_flux);
1017 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"H", 0, 0,
1018 is_H);
1019 IS is_tmp;
1020 CHKERR ISExpand(is_H, is_flux, &is_tmp);
1022
1023 auto is_all_bc = bc_mng->getBlockIS(name_prb, "FLUIDFLUX", "FLUX", 0, 0);
1024 int is_all_bc_size;
1025 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1027 << "Field split block size " << is_all_bc_size;
1028 if (is_all_bc_size) {
1029 IS is_tmp2;
1030 CHKERR ISDifference(is_Flux, is_all_bc, &is_tmp2);
1032 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR,
1033 is_all_bc);
1034 }
1035
1038 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_Flux);
1039 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1040 }
1041
1043 };
1044
1045 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1046 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1047 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1048 auto time_scale = boost::make_shared<TimeScale>();
1049
1050 auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
1052 {time_scale}, false);
1053 return hook;
1054 };
1055
1056 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1059 mField, post_proc_rhs_ptr,
nullptr, Sev::verbose)();
1061 mField, post_proc_rhs_ptr, 1.)();
1063 };
1064 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1066 post_proc_lhs_ptr, 1.);
1067 };
1068
1069 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1070 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1071 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1072
1074 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1075 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1076 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1077 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1078
1080 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1082 CHKERR TSSetSolution(solver,
D);
1083 CHKERR TSSetFromOptions(solver);
1084
1085 CHKERR set_section_monitor(solver);
1086 CHKERR set_fieldsplit_preconditioner(solver);
1087 CHKERR set_time_monitor(dm, solver);
1088
1090 CHKERR TSSolve(solver, NULL);
1091
1093}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
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.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to calculate residual side diagonal.
Class (Function) to enforce essential constrains.
Field evaluator interface.
Section manager is used to create indexes and sections.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
intrusive_ptr for managing petsc objects