556 {
560
561 CHKERR KSPSetFromOptions(ksp);
562 PC pc;
563 CHKERR KSPGetPC(ksp, &pc);
564
565 PetscBool is_pcfs = PETSC_FALSE;
566 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
567 if (is_pcfs) {
568
569 MOFEM_LOG(
"AT", Sev::inform) <<
"Setup Schur pc";
570
571 auto create_sub_dm = [&]() {
573
574 auto create_dm = [&](
575
576 std::string problem_name,
577 std::vector<std::string> fe_names,
578 std::vector<std::string> fields,
579
580 auto dm_type
581
582 ) {
584 auto create_dm_imp = [&]() {
588 for (auto fe : fe_names) {
590 }
592 for (auto field : fields) {
595 }
598 };
600 create_dm_imp(),
601 "Error in creating schurDM. It is possible that schurDM is "
602 "already created");
603 return dm;
604 };
605
606 auto schur_dm = create_dm(
607
608 "SCHUR",
609
610 {
simple->getDomainFEName(),
simple->getSkeletonFEName()},
611
612 {"HYBRID"},
613
614 "DMMOFEM_MG");
615
616 auto block_dm = create_dm(
617
618 "BLOCK",
619
620 {
simple->getDomainFEName(),
simple->getSkeletonFEName()},
621
622 {"BROKEN", "U"},
623
624 "DMMOFEM");
625
626 return std::make_tuple(schur_dm, block_dm);
627 };
628
629 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
632
633 {
634
635 {
636
637 simple->getDomainFEName(),
638
639 {
640
641 {"BROKEN", "BROKEN"},
642 {"U", "U"},
643 {"BROKEN", "U"},
644 {"U", "BROKEN"}
645
646 }
647
648 },
649
650 {
651
652 simple->getSkeletonFEName(),
653
654 {
655
656 {"BROKEN", "HYBRID"}, {"HYBRID", "BROKEN"}
657
658 }
659
660 }
661
662 }
663
664 );
665
667
668 {schur_dm, block_dm}, block_mat_data,
669
670 {"BROKEN", "U"}, {nullptr, nullptr}, true
671
672 );
673 };
674
675 auto set_ops = [&](auto schur_dm) {
679
680 boost::shared_ptr<BlockStructure> block_data;
682
684 pip_mng->getOpDomainLhsPipeline().push_front(
686 pip_mng->getOpDomainLhsPipeline().push_back(
687
690
691 );
692 }
693
694 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
695 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
696
697 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
700 MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble Begin";
702 };
703
704 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
705 post_proc_schur_lhs_ptr]() {
707 MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble End";
708
710 CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
711 CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
712 if (1) {
713 auto S_from_block =
matDuplicate(
S, MAT_SHARE_NONZERO_PATTERN);
714
716 S_from_block, {"BROKEN", "U"},
717 {nullptr, nullptr}, ao_up);
718 CHKERR MatAssemblyBegin(S_from_block, MAT_FINAL_ASSEMBLY);
719 CHKERR MatAssemblyEnd(S_from_block, MAT_FINAL_ASSEMBLY);
720 CHKERR MatAYPX(S_from_block, -1,
S, DIFFERENT_NONZERO_PATTERN);
721 double norm;
722 CHKERR MatNorm(S_from_block, NORM_FROBENIUS, &norm);
723 MOFEM_LOG(
"AT", Sev::inform) <<
"Norm of difference: " << norm;
724 if (norm > 1e-6)
725 SETERRQ(
727 "Norm of difference between Schur and block matrix is larger "
728 "than accepted");
729 }
730 } else {
732 {"BROKEN", "U"}, {nullptr, nullptr}, ao_up);
733 CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
734 CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
735 }
736
737 MOFEM_LOG(
"AT", Sev::verbose) <<
"Lhs Assemble Finish";
739 };
740
742 ksp_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_schur_lhs_ptr);
743 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_schur_lhs_ptr);
744
746 };
747
748 auto set_pc = [&](auto pc, auto block_dm) {
750 auto block_is =
getDMSubData(block_dm)->getSmartRowIs();
751 CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
752 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
754 };
755
756
757
758 auto [schur_dm, block_dm] = create_sub_dm();
760 auto nested_mat_data = get_nested_mat_data(schur_dm, block_dm);
762 }
764 CHKERR MatSetDM(
S, PETSC_NULLPTR);
765
769
771 CHKERR set_pc(pc, block_dm);
772 DM solver_dm;
773 CHKERR KSPGetDM(ksp, &solver_dm);
775 CHKERR DMSetMatType(solver_dm, MATSHELL);
776
777 auto get_pc = [](auto ksp) {
778 PC pc_raw;
779 CHKERR KSPGetPC(ksp, &pc_raw);
780 return pc_raw;
781 };
782
783 auto set_diagonal_pc = [&](auto pc) {
785
789 CHKERR PCSetOperators(pc,
A, P);
790 }
791
792 KSP *subksp;
793 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
795
798 };
799
800 auto set_mg_for_schur_complement = [&](
auto pc,
auto schur_dm,
auto S) {
802
803 KSP *subksp;
804 CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
805 auto subpc = get_pc(subksp[1]);
806
807 auto set_pc_p_mg = [](
auto dm,
auto pc,
auto S) {
810 PetscBool same = PETSC_FALSE;
811 PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
812 if (same) {
813 MOFEM_LOG(
"TIMER", Sev::inform) <<
"Set up MG";
816 CHKERR PCSetFromOptions(pc);
817 }
819 };
820
821 PetscBool same = PETSC_FALSE;
822 PetscObjectTypeCompare((PetscObject)subpc, PCKSP, &same);
823 if (same) {
824 CHKERR PCSetFromOptions(subpc);
825 KSP inner_ksp;
826 CHKERR PCKSPGetKSP(subpc, &inner_ksp);
827 CHKERR KSPSetFromOptions(inner_ksp);
828 PC inner_pc;
829 CHKERR KSPGetPC(inner_ksp, &inner_pc);
830 CHKERR PCSetFromOptions(inner_pc);
831 CHKERR set_pc_p_mg(schur_dm, inner_pc,
S);
832 }
833
836 };
837
840 CHKERR set_diagonal_pc(pc);
841 CHKERR set_mg_for_schur_complement(pc, schur_dm,
S);
842 }
843
844 } else {
846 "PC is not set to PCFIELDSPLIT");
847 }
849}
@ MOFEM_ATOM_TEST_INVALID
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
#define NBEDGE_L2(P)
Number of base functions on edge from L2 space.
MoFEMErrorCode DMMoFEMGetBlocMatData(DM dm, boost::shared_ptr< BlockStructure > &)
Get data for block mat.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
auto createDMHybridisedL2Matrix(DM dm)
Get smart hybridised L2 matrix from DM.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDMNestSchurMat(DM dm)
MoFEMErrorCode assembleBlockMatSchur(MoFEM::Interface &m_field, Mat B, Mat S, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao)
Assemble Schur matrix.
auto createDMBlockMat(DM dm)
static constexpr int approx_order
constexpr AssemblyType A
[Define dimension]