569 {
571
575
576 auto set_section_monitor = [&](auto solver) {
578 SNES snes;
579 CHKERR TSGetSNES(solver, &snes);
580 PetscViewerAndFormat *vf;
581 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
582 PETSC_VIEWER_DEFAULT, &vf);
584 snes,
586 void *))SNESMonitorFields,
589 };
590
591 auto scatter_create = [&](
auto D,
auto coeff) {
593 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
594 ROW,
"U", coeff, coeff, is);
595 int loc_size;
596 CHKERR ISGetLocalSize(is, &loc_size);
599 VecScatter scatter;
600 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULLPTR, &scatter);
603 };
604
605 auto create_post_process_elements = [&]() {
606 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
607
608 auto push_vol_ops = [this](auto &pip) {
611 "GEOMETRY");
612
614 };
615
616 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
618
619 auto &pip = pp_fe->getOpPtrVector();
620
622
623 auto x_ptr = boost::make_shared<MatrixDouble>();
624 pip.push_back(
626 auto u_ptr = boost::make_shared<MatrixDouble>();
628
629 auto pressure_ptr = boost::make_shared<VectorDouble>();
631
632 auto div_u_ptr = boost::make_shared<VectorDouble>();
634 "U", div_u_ptr));
635
636 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
638 "U", grad_u_ptr));
639
640 auto strain_ptr = boost::make_shared<MatrixDouble>();
642
643 auto stress_ptr = boost::make_shared<MatrixDouble>();
645 mu, stress_ptr, strain_ptr, pressure_ptr));
646
647 pip.push_back(
648
650
651 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
652
653 {{"P", pressure_ptr}},
654
655 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
656
657 {},
658
659 {{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
660
661 )
662
663 );
664
666 };
667
668 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
669 PetscBool post_proc_vol = PETSC_FALSE;
671 &post_proc_vol, PETSC_NULLPTR);
672 if (post_proc_vol == PETSC_FALSE)
673 return boost::shared_ptr<PostProcEle>();
674 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
676 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
677 "push_vol_post_proc_ops");
678 return pp_fe;
679 };
680
681 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
682 PetscBool post_proc_skin = PETSC_TRUE;
684 &post_proc_skin, PETSC_NULLPTR);
685 if (post_proc_skin == PETSC_FALSE)
686 return boost::shared_ptr<SkinPostProcEle>();
687
689 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
692 pp_fe->getOpPtrVector().push_back(op_side);
694 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
695 "push_vol_post_proc_ops");
696 return pp_fe;
697 };
698
699 return std::make_pair(vol_post_proc(), skin_post_proc());
700 };
701
702 boost::shared_ptr<SetPtsData> field_eval_data;
703 boost::shared_ptr<MatrixDouble> vector_field_ptr;
704
706 vector_field_ptr = boost::make_shared<MatrixDouble>();
707 field_eval_data =
710 field_eval_data,
simple->getDomainFEName());
711
713 auto no_rule = [](
int,
int,
int) {
return -1; };
714 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
715 field_eval_fe_ptr->getRuleHook = no_rule;
716 field_eval_fe_ptr->getOpPtrVector().push_back(
718 }
719
720 auto set_time_monitor = [&](auto dm, auto solver) {
722 boost::shared_ptr<MonitorIncompressible> monitor_ptr(
724 create_post_process_elements(),
uXScatter,
726 field_eval_data, vector_field_ptr));
727 boost::shared_ptr<ForcesAndSourcesCore> null;
729 monitor_ptr, null, null);
731 };
732
733 auto set_essential_bc = [&]() {
735
737 auto pre_proc_ptr = boost::make_shared<FEMethod>();
738 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
739
740
741 auto time_scale = boost::make_shared<TimeScale>();
742
744 mField, pre_proc_ptr, {time_scale},
false);
745 post_proc_rhs_ptr->postProcessHook =
747 1.);
748 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
749 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
750 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
752 };
753
754 auto set_schur_pc = [&](auto solver) {
755 SNES snes;
756 CHKERR TSGetSNES(solver, &snes);
757 KSP ksp;
758 CHKERR SNESGetKSP(snes, &ksp);
759 PC pc;
760 CHKERR KSPGetPC(ksp, &pc);
761 PetscBool is_pcfs = PETSC_FALSE;
762 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
763 boost::shared_ptr<SetUpSchur> schur_ptr;
765
768
771 "set operators");
772 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
773 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
774 pre_proc_schur_lhs_ptr->preProcessHook = [pre_proc_schur_lhs_ptr]() {
776 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Zero matrices";
777 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->A);
778 CHKERR MatZeroEntries(pre_proc_schur_lhs_ptr->B);
780 };
781 post_proc_schur_lhs_ptr->postProcessHook = [this,
782 post_proc_schur_lhs_ptr]() {
784 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Assemble Begin";
785 *(post_proc_schur_lhs_ptr->matAssembleSwitch) = false;
786 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
787 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->A, MAT_FINAL_ASSEMBLY);
789 mField, post_proc_schur_lhs_ptr, 1.,
791
792 CHKERR MatAssemblyBegin(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
793 CHKERR MatAssemblyEnd(post_proc_schur_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
794 CHKERR MatAXPY(post_proc_schur_lhs_ptr->B, 1, post_proc_schur_lhs_ptr->A,
795 SAME_NONZERO_PATTERN);
796 MOFEM_LOG(
"INCOMP_ELASTICITY", Sev::verbose) <<
"Lhs Assemble End";
798 };
799 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
800 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
801
802 if (is_pcfs == PETSC_TRUE) {
803 if (
AT == AssemblyType::SCHUR) {
805 CHK_MOAB_THROW(schur_ptr->setUp(solver),
"setup schur preconditioner");
806 } else {
807 auto set_pcfieldsplit_preconditioned_ts = [&](auto solver) {
809 auto name_prb =
simple->getProblemName();
812 name_prb,
ROW,
"P", 0, 1, is_p);
813 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_p);
815 };
817 "set pcfieldsplit preconditioned");
818 }
819 return boost::make_tuple(schur_ptr, A,
B);
820 }
821
822 return boost::make_tuple(schur_ptr, A,
B);
823 };
824
825 auto dm =
simple->getDM();
827
832
833
834
835 CHKERR set_essential_bc();
836
837 auto solver = pip_mng->createTSIM();
838 CHKERR TSSetFromOptions(solver);
839
840 CHKERR set_section_monitor(solver);
841 CHKERR set_time_monitor(dm, solver);
842 auto [schur_pc_ptr,
A,
B] = set_schur_pc(solver);
843
844 CHKERR TSSetSolution(solver,
D);
846 CHKERR TSSolve(solver, NULL);
847
849}
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
const double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
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 enforce essential constrains.
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
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
