v0.14.0
Public Member Functions | Private Member Functions | Private Attributes | List of all members
SetUpSchurImpl Struct Reference
Inheritance diagram for SetUpSchurImpl:
[legend]
Collaboration diagram for SetUpSchurImpl:
[legend]

Public Member Functions

 SetUpSchurImpl (MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (KSP solver)
 
MoFEMErrorCode preProc ()
 
MoFEMErrorCode postProc ()
 
 SetUpSchurImpl (MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (TS solver)
 
MoFEMErrorCode preProc ()
 
MoFEMErrorCode postProc ()
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< TS > solver)
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< TS > solver)
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< KSP > solver)
 

Private Member Functions

MoFEMErrorCode setEntities ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 
SmartPetscObj< DM > createSubDM ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 
SmartPetscObj< DM > createSubDM ()
 
MoFEMErrorCode setEntities ()
 
MoFEMErrorCode setUpSubDM ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 

Private Attributes

SmartPetscObj< Mat > M
 
SmartPetscObj< Mat > P
 
SmartPetscObj< Mat > S
 
SmartPetscObj< AO > aoUp
 
SmartPetscObj< IS > isA00
 
SmartPetscObj< IS > isShift
 
MoFEM::InterfacemField
 
EshelbianCore * epCorePtr
 
SmartPetscObj< DM > subDM
 field split sub dm More...
 
SmartPetscObj< IS > fieldSplitIS
 IS for split Schur block. More...
 
SmartPetscObj< Mat > A
 
Range volEnts
 
Range subEnts
 

Detailed Description

Examples
plastic.cpp.

Definition at line 9 of file SetUpSchurImpl.cpp.

Constructor & Destructor Documentation

◆ SetUpSchurImpl() [1/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field,
SmartPetscObj< Mat >  m,
SmartPetscObj< Mat >  p,
EshelbianCore *  ep_core_ptr 
)
inline

Definition at line 11 of file SetUpSchurImpl.cpp.

13  : SetUpSchur(), mField(m_field), M(m), P(p), epCorePtr(ep_core_ptr) {
14  if (S) {
17  "Is expected that schur matrix is not allocated. This is "
18  "possible only is PC is set up twice");
19  }
20  }

◆ ~SetUpSchurImpl() [1/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 21 of file SetUpSchurImpl.cpp.

21 { S.reset(); }

◆ SetUpSchurImpl() [2/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field,
SmartPetscObj< DM >  sub_dm,
SmartPetscObj< IS >  field_split_is,
SmartPetscObj< AO >  ao_up 
)
inline

Definition at line 1273 of file plastic.cpp.

1275  : SetUpSchur(), mField(m_field), subDM(sub_dm),
1276  fieldSplitIS(field_split_is), aoUp(ao_up) {
1277  if (S) {
1280  "Is expected that schur matrix is not allocated. This is "
1281  "possible only is PC is set up twice");
1282  }
1283  }

◆ ~SetUpSchurImpl() [2/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 1284 of file plastic.cpp.

1284  {
1285 #ifdef ADD_CONTACT
1286  A.reset();
1287  P.reset();
1288 #endif // ADD_CONTACT
1289  S.reset();
1290  }

◆ SetUpSchurImpl() [3/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 995 of file contact.cpp.

995 : SetUpSchur(), mField(m_field) {}

◆ ~SetUpSchurImpl() [3/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 997 of file contact.cpp.

997  {
998  A.reset();
999  P.reset();
1000  S.reset();
1001  }

◆ SetUpSchurImpl() [4/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 935 of file incompressible_elasticity.cpp.

935 : SetUpSchur(), mField(m_field) {}

◆ ~SetUpSchurImpl() [4/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 936 of file incompressible_elasticity.cpp.

936 { S.reset(); }

◆ SetUpSchurImpl() [5/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 883 of file elastic.cpp.

883  : SetUpSchur(), mField(m_field) {
884  if (S) {
887  "Is expected that schur matrix is not allocated. This is "
888  "possible only is PC is set up twice");
889  }
890  }

◆ ~SetUpSchurImpl() [5/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 891 of file elastic.cpp.

891 { S.reset(); }

Member Function Documentation

◆ createSubDM() [1/2]

SmartPetscObj<DM> SetUpSchurImpl::createSubDM ( )
private

◆ createSubDM() [2/2]

SmartPetscObj< DM > SetUpSchurImpl::createSubDM ( )
private

Definition at line 1085 of file contact.cpp.

1085  {
1086  auto simple = mField.getInterface<Simple>();
1087  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1088  auto set_up = [&]() {
1090  CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
1091  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1092  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1093  CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
1094  CHKERR DMSetUp(sub_dm);
1096  };
1097  CHK_THROW_MESSAGE(set_up(), "sey up dm");
1098  return sub_dm;
1099 }

◆ postProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::postProc ( )

Definition at line 213 of file SetUpSchurImpl.cpp.

213  {
215  if (S) {
216  MOFEM_LOG("EP", Sev::verbose) << "Assemble Schur";
217  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
218  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
219  }
220  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
221  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
222  CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
223  CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
224  CHKERR MatAXPY(P, 1, M, SAME_NONZERO_PATTERN);
225 
227 }

◆ postProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::postProc ( )

◆ preProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::preProc ( )

Definition at line 202 of file SetUpSchurImpl.cpp.

202  {
204  MOFEM_LOG("EP", Sev::verbose) << "Zero Schur";
205  if (SetUpSchurImpl::S) {
206  CHKERR MatZeroEntries(S);
207  }
208  CHKERR MatZeroEntries(M);
209  CHKERR MatZeroEntries(P);
211 }

◆ preProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::preProc ( )

◆ setEntities() [1/2]

MoFEMErrorCode SetUpSchurImpl::setEntities ( )
private

◆ setEntities() [2/2]

MoFEMErrorCode SetUpSchurImpl::setEntities ( )
private

Definition at line 942 of file elastic.cpp.

942  {
944  auto simple = mField.getInterface<Simple>();
945  CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
946  SPACE_DIM, volEnts);
947  CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
948  subEnts);
949  subEnts = subtract(subEnts, volEnts);
951 };

◆ setOperator() [1/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

◆ setOperator() [2/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

◆ setOperator() [3/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

Definition at line 1101 of file contact.cpp.

1101  {
1103 
1104  double eps_stab = 1e-4;
1105  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1106  PETSC_NULL);
1107 
1108  using B =
1110  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1111 
1112  auto pip = mField.getInterface<PipelineManager>();
1113  // Boundary
1114  auto dm_is = getDMSubData(subDM)->getSmartRowIs();
1115  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
1116 
1117  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1118  pip->getOpBoundaryLhsPipeline().push_back(
1119  new OpMassStab("SIGMA", "SIGMA",
1120  [eps_stab](double, double, double) { return eps_stab; }));
1121  pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1122  {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1123 
1124  // Domain
1125  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1126  pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1127  {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1128 
1129  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1130  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1131 
1132  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1134  CHKERR MatZeroEntries(A);
1135  CHKERR MatZeroEntries(P);
1136  CHKERR MatZeroEntries(S);
1137  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1139  };
1140 
1141  post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1142  post_proc_schur_lhs_ptr]() {
1144  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1145 
1146  *post_proc_schur_lhs_ptr->matAssembleSwitch = false;
1147 
1148  auto print_mat_norm = [this](auto a, std::string prefix) {
1150  double nrm;
1151  CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1152  MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1154  };
1155 
1156  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1157  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1159  mField, post_proc_schur_lhs_ptr, 1, A)();
1160 
1161  CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1162  CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1163  CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1164 
1165  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1166  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1167 
1169  mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1170 
1171 #ifndef NDEBUG
1172  CHKERR print_mat_norm(A, "A");
1173  CHKERR print_mat_norm(P, "P");
1174  CHKERR print_mat_norm(S, "S");
1175 #endif // NDEBUG
1176 
1177  MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1179  };
1180 
1181  auto simple = mField.getInterface<Simple>();
1182  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1183  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1184  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1185 
1187 }

◆ setPC() [1/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

◆ setPC() [2/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

◆ setPC() [3/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

Definition at line 1189 of file contact.cpp.

1189  {
1191  auto simple = mField.getInterface<Simple>();
1192  SmartPetscObj<IS> is;
1193  mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1194  simple->getProblemName(), ROW, "SIGMA", 0, SPACE_DIM, is);
1195  CHKERR PCFieldSplitSetIS(pc, NULL, is);
1196  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1198 }

◆ setUp() [1/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( KSP  solver)
Examples
plastic.cpp.

Definition at line 39 of file SetUpSchurImpl.cpp.

39  {
41 
42  auto get_ents_by_dim = [&](int dim) {
43  Range ents;
44  CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents);
45  return ents;
46  };
47 
48  std::vector<std::string> schur_field_list{
49 
50  epCorePtr->stretchTensor,
51 
52  epCorePtr->bubbleField,
53 
54  epCorePtr->rotAxis,
55 
56  epCorePtr->spatialL2Disp,
57 
58  epCorePtr->piolaStress
59 
60  };
61  std::vector<std::string> field_list{epCorePtr->piolaStress,
62  epCorePtr->contactDisp};
63 
64  auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
66 
67  dm_sub = createDM(mField.get_comm(), "DMMOFEM");
68  CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
69  CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
70  CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
71  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
72  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->naturalBcElement);
73  CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->contactElement);
74 
75  auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM - 1));
76  std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr, nullptr};
77 
78  if (field_list.size() != dm_range_list.size()) {
79  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
80  "Both ranges should have the same size");
81  }
82 
83  int r_idx = 0;
84  for (auto f : field_list) {
85  MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
86  CHKERR DMMoFEMAddSubFieldRow(dm_sub, f, dm_range_list[r_idx]);
87  CHKERR DMMoFEMAddSubFieldCol(dm_sub, f, dm_range_list[r_idx]);
88  ++r_idx;
89  }
90  CHKERR DMSetUp(dm_sub);
92  };
93 
94  auto create_a00_is = [&](SmartPetscObj<IS> &is_a00) {
96  auto vols = get_ents_by_dim(SPACE_DIM);
97  std::vector<SmartPetscObj<IS>> is_vec;
98  std::vector<Range *> range_list_ptr(schur_field_list.size(), nullptr);
99  range_list_ptr.back() = &vols;
100 
101  int r_idx = 0;
102  for (auto f : schur_field_list) {
103  SmartPetscObj<IS> is;
104  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
105  "ELASTIC_PROBLEM", ROW, f, 0, MAX_DOFS_ON_ENTITY, is,
106  range_list_ptr[r_idx]);
107  is_vec.push_back(is);
108  ++r_idx;
109  }
110  if (is_vec.size()) {
111  is_a00 = is_vec[0];
112  for (auto i = 1; i < is_vec.size(); ++i) {
113  IS is_union_raw;
114  CHKERR ISExpand(is_a00, is_vec[i], &is_union_raw);
115  is_a00 = SmartPetscObj<IS>(is_union_raw);
116  }
117  CHKERR ISSort(is_a00);
118  }
120  };
121 
122  PC pc;
123  CHKERR KSPSetFromOptions(solver);
124  CHKERR KSPGetPC(solver, &pc);
125  PetscBool is_pcfs = PETSC_FALSE;
126  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
127  if (is_pcfs) {
128 
129  SmartPetscObj<DM> schur_dm, a00_dm;
130  CHKERR create_schur_dm(schur_dm);
131  CHKERR create_a00_is(isA00);
132 
133  if (S) {
136  "Is expected that schur matrix is not allocated. This is "
137  "possible only is PC is set up twice");
138  }
139 
140  S = createDMMatrix(schur_dm);
141 
142  auto set_ops = [&]() {
144  auto range_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM));
145  std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
146  nullptr);
147  ranges_list.back() = range_ptr;
148  std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
149  SmartPetscObj<AO>());
150  std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
151  SmartPetscObj<Mat>());
152  std::vector<bool> symm_list(schur_field_list.size(), false);
153 
154  auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
155  aoUp = createAOMappingIS(dm_is, PETSC_NULL);
156  ao_list.back() = aoUp;
157  mat_list.back() = S;
158 
159  const bool symmetric_system =
160  EshelbianCore::gradApperoximator <= MODERATE_ROT &&
161  EshelbianCore::rotSelector == SMALL_ROT;
162  epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
163  new OpSchurAssembleBegin());
164  epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
165  new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
166  ao_list, mat_list, symm_list,
167  symmetric_system));
168  epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
169  new OpSchurAssembleBegin());
170  epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
171  new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
172  ao_list, mat_list, symm_list,
173  symmetric_system));
174 
176  };
177 
178  auto set_pc = [&]() {
180  CHKERR PCFieldSplitSetIS(pc, NULL, isA00);
181  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
183  };
184 
185  CHKERR set_ops();
186  CHKERR set_pc();
187 
188  } else {
189  epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
190  new OpSchurAssembleBegin());
191  epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
192  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
193  epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
194  new OpSchurAssembleBegin());
195  epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
196  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
197  }
198 
200 }

◆ setUp() [2/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< KSP >  solver)

Definition at line 910 of file elastic.cpp.

910  {
912  auto pip = mField.getInterface<PipelineManager>();
913  PC pc;
914  CHKERR KSPGetPC(solver, &pc);
915  PetscBool is_pcfs = PETSC_FALSE;
916  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
917  if (is_pcfs) {
918  if (S) {
921  "Is expected that schur matrix is not allocated. This is "
922  "possible only is PC is set up twice");
923  }
925  CHKERR setUpSubDM();
927  CHKERR MatSetBlockSize(S, SPACE_DIM);
928  CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
930  CHKERR setPC(pc);
931  } else {
932  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
933  pip->getOpBoundaryLhsPipeline().push_back(
934  new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
935  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
936  pip->getOpDomainLhsPipeline().push_back(
937  new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
938  }
940 }

◆ setUp() [3/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver)

◆ setUp() [4/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver)

Definition at line 1021 of file contact.cpp.

1021  {
1023  auto simple = mField.getInterface<Simple>();
1024  auto pip = mField.getInterface<PipelineManager>();
1025  auto dm = simple->getDM();
1026 
1027  SNES snes;
1028  CHKERR TSGetSNES(solver, &snes);
1029  KSP ksp;
1030  CHKERR SNESGetKSP(snes, &ksp);
1031  CHKERR KSPSetFromOptions(ksp);
1032 
1033  PC pc;
1034  CHKERR KSPGetPC(ksp, &pc);
1035 
1036  PetscBool is_pcfs = PETSC_FALSE;
1037  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1038  if (is_pcfs) {
1039 
1040  MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
1041 
1042  if (A || P || S) {
1045  "It is expected that Schur matrix is not allocated. This is "
1046  "possible only if PC is set up twice");
1047  }
1048 
1049  A = createDMMatrix(dm);
1050  P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1051  subDM = createSubDM();
1052  S = createDMMatrix(subDM);
1053  CHKERR MatSetBlockSize(S, SPACE_DIM);
1054 
1055  auto ts_ctx_ptr = getDMTsCtx(dm);
1056  CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1057 
1058  CHKERR setOperator();
1059  CHKERR setPC(pc);
1060 
1061  } else {
1062  MOFEM_LOG("CONTACT", Sev::inform) << "No Schur pc";
1063 
1064  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1065  pip->getOpBoundaryLhsPipeline().push_back(
1066  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1067  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1068  pip->getOpDomainLhsPipeline().push_back(
1069  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1070 
1071  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1072  post_proc_schur_lhs_ptr->postProcessHook = [this,
1073  post_proc_schur_lhs_ptr]() {
1076  mField, post_proc_schur_lhs_ptr, 1.)();
1078  };
1079  auto ts_ctx_ptr = getDMTsCtx(dm);
1080  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1081  }
1083 }

◆ setUp() [5/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( TS  solver)

Definition at line 1309 of file plastic.cpp.

1309  {
1311  auto simple = mField.getInterface<Simple>();
1312  auto pip_mng = mField.getInterface<PipelineManager>();
1313 
1314  SNES snes;
1315  CHKERR TSGetSNES(solver, &snes);
1316  KSP ksp;
1317  CHKERR SNESGetKSP(snes, &ksp);
1318  CHKERR KSPSetFromOptions(ksp);
1319 
1320  PC pc;
1321  CHKERR KSPSetFromOptions(ksp);
1322  CHKERR KSPGetPC(ksp, &pc);
1323  PetscBool is_pcfs = PETSC_FALSE;
1324  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1325  if (is_pcfs) {
1326  if (S) {
1329  "Is expected that schur matrix is not allocated. This is "
1330  "possible only is PC is set up twice");
1331  }
1332 
1333 #ifdef ADD_CONTACT
1334  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1335  A = createDMMatrix(simple->getDM());
1336  P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1337  CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1338 #endif // ADD_CONTACT
1339  S = createDMMatrix(subDM);
1340  CHKERR MatSetBlockSize(S, SPACE_DIM);
1341 
1342  auto set_ops = [&]() {
1344  auto pip_mng = mField.getInterface<PipelineManager>();
1345 
1346 #ifndef ADD_CONTACT
1347  // Boundary
1348  pip_mng->getOpBoundaryLhsPipeline().push_front(
1349  new OpSchurAssembleBegin());
1350  pip_mng->getOpBoundaryLhsPipeline().push_back(
1352 
1353  {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1354  {SmartPetscObj<Mat>(), S}, {false, false}
1355 
1356  ));
1357  // Domain
1358  pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1359  pip_mng->getOpDomainLhsPipeline().push_back(
1361 
1362  {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1363  {SmartPetscObj<Mat>(), S}, {false, false}
1364 
1365  ));
1366 #else
1367 
1368  double eps_stab = 1e-4;
1369  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1370  PETSC_NULL);
1371 
1374  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1375 
1376  // Boundary
1377  pip_mng->getOpBoundaryLhsPipeline().push_front(
1378  new OpSchurAssembleBegin());
1379  pip_mng->getOpBoundaryLhsPipeline().push_back(
1380  new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1381  return eps_stab;
1382  }));
1383  pip_mng->getOpBoundaryLhsPipeline().push_back(
1385 
1386  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1389  {false, false, false}
1390 
1391  ));
1392  // Domain
1393  pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1394  pip_mng->getOpDomainLhsPipeline().push_back(
1396 
1397  {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1400  {false, false, false}
1401 
1402  ));
1403 #endif // ADD_CONTACT
1405  };
1406 
1407  auto set_assemble_elems = [&]() {
1409  auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1410  schur_asmb_pre_proc->preProcessHook = [this]() {
1412 #ifdef ADD_CONTACT
1413  CHKERR MatZeroEntries(A);
1414  CHKERR MatZeroEntries(P);
1415 #endif // ADD_CONTACT
1416  CHKERR MatZeroEntries(S);
1417  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1419  };
1420  auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1421 
1422  schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1424  MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1425 
1426 #ifndef ADD_CONTACT
1428  mField, schur_asmb_post_proc, 1)();
1429 #else // ADD_CONTACT
1430  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1431  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1432  // Apply essential constrains to A matrix
1434  mField, schur_asmb_post_proc, 1, A)();
1435  CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1436  CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1437  CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1438 #endif // ADD_CONTACT
1439 
1440  // Apply essential constrains to Schur complement
1441  CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1442  CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1444  mField, schur_asmb_post_proc, 1, S, aoUp)();
1445 
1447  };
1448  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1449  ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1450  ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1452  };
1453 
1454  auto set_pc = [&]() {
1456  CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1457  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1459  };
1460 
1461  CHKERR set_ops();
1462  CHKERR set_pc();
1463  CHKERR set_assemble_elems();
1464 
1465  } else {
1466  pip_mng->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1467  pip_mng->getOpBoundaryLhsPipeline().push_back(
1468  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1469  pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1470  pip_mng->getOpDomainLhsPipeline().push_back(
1471  new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1472  }
1473 
1474  // we do not those anymore
1475  subDM.reset();
1476  fieldSplitIS.reset();
1477  // aoUp.reset();
1479 }

◆ setUpSubDM()

MoFEMErrorCode SetUpSchurImpl::setUpSubDM ( )
private

Definition at line 953 of file elastic.cpp.

953  {
955  auto simple = mField.getInterface<Simple>();
956  subDM = createDM(mField.get_comm(), "DMMOFEM");
957  CHKERR DMMoFEMCreateSubDM(subDM, simple->getDM(), "SUB");
958  CHKERR DMMoFEMSetSquareProblem(subDM, PETSC_TRUE);
959  CHKERR DMMoFEMAddElement(subDM, simple->getDomainFEName());
960  auto sub_ents_ptr = boost::make_shared<Range>(subEnts);
961  CHKERR DMMoFEMAddSubFieldRow(subDM, "U", sub_ents_ptr);
962  CHKERR DMMoFEMAddSubFieldCol(subDM, "U", sub_ents_ptr);
963  CHKERR DMSetUp(subDM);
965 }

Member Data Documentation

◆ A

SmartPetscObj<Mat> SetUpSchurImpl::A
private

Definition at line 1012 of file contact.cpp.

◆ aoUp

SmartPetscObj< AO > SetUpSchurImpl::aoUp
private

Definition at line 31 of file SetUpSchurImpl.cpp.

◆ epCorePtr

EshelbianCore* SetUpSchurImpl::epCorePtr
private

Definition at line 36 of file SetUpSchurImpl.cpp.

◆ fieldSplitIS

SmartPetscObj<IS> SetUpSchurImpl::fieldSplitIS
private

IS for split Schur block.

Definition at line 1305 of file plastic.cpp.

◆ isA00

SmartPetscObj<IS> SetUpSchurImpl::isA00
private

Definition at line 32 of file SetUpSchurImpl.cpp.

◆ isShift

SmartPetscObj<IS> SetUpSchurImpl::isShift
private

Definition at line 33 of file SetUpSchurImpl.cpp.

◆ M

SmartPetscObj<Mat> SetUpSchurImpl::M
private

Definition at line 28 of file SetUpSchurImpl.cpp.

◆ mField

MoFEM::Interface & SetUpSchurImpl::mField
private

Definition at line 35 of file SetUpSchurImpl.cpp.

◆ P

SmartPetscObj< Mat > SetUpSchurImpl::P
private

Definition at line 29 of file SetUpSchurImpl.cpp.

◆ S

SmartPetscObj< Mat > SetUpSchurImpl::S
private

Definition at line 30 of file SetUpSchurImpl.cpp.

◆ subDM

SmartPetscObj< DM > SetUpSchurImpl::subDM
private

field split sub dm

Definition at line 1304 of file plastic.cpp.

◆ subEnts

Range SetUpSchurImpl::subEnts
private

Definition at line 907 of file elastic.cpp.

◆ volEnts

Range SetUpSchurImpl::volEnts
private

Definition at line 906 of file elastic.cpp.


The documentation for this struct was generated from the following files:
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:302
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:284
SetUpSchurImpl::createSubDM
SmartPetscObj< DM > createSubDM()
Definition: contact.cpp:1085
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:130
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
SetUpSchurImpl::epCorePtr
EshelbianCore * epCorePtr
Definition: SetUpSchurImpl.cpp:36
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:460
SetUpSchurImpl::M
SmartPetscObj< Mat > M
Definition: SetUpSchurImpl.cpp:28
SetUpSchurImpl::setUpSubDM
MoFEMErrorCode setUpSubDM()
Definition: elastic.cpp:953
SetUpSchurImpl::isA00
SmartPetscObj< IS > isA00
Definition: SetUpSchurImpl.cpp:32
SetUpSchurImpl::setPC
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1189
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::OpSchurAssembleEnd< SCHUR_DGESV >
Definition: Schur.hpp:114
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1076
SetUpSchurImpl::subDM
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1003
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
SPACE_DIM
constexpr int SPACE_DIM
Definition: contact.cpp:53
SetUpSchurImpl::volEnts
Range volEnts
Definition: elastic.cpp:906
ROW
@ ROW
Definition: definitions.h:123
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
SetUpSchurImpl::S
SmartPetscObj< Mat > S
Definition: SetUpSchurImpl.cpp:30
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::OpSchurAssembleBegin
Clear Schur complement internal data.
Definition: Schur.hpp:22
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:137
a
constexpr double a
Definition: approx_sphere.cpp:30
SetUpSchurImpl::fieldSplitIS
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1305
SPACE_DIM
constexpr int SPACE_DIM
Definition: plastic.cpp:40
SetUpSchurImpl::setEntities
MoFEMErrorCode setEntities()
Definition: elastic.cpp:942
MoFEM::OpSchurAssembleEnd< SCHUR_DSYSV >
Definition: Schur.hpp:107
SetUpSchurImpl::subEnts
Range subEnts
Definition: elastic.cpp:907
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
SetUpSchurImpl::A
SmartPetscObj< Mat > A
Definition: contact.cpp:1012
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:302
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: SetUpSchurImpl.cpp:35
MoFEM::TsSetIJacobian
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.
Definition: TsCtx.cpp:165
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
SetUpSchurImpl::setOperator
MoFEMErrorCode setOperator()
Definition: contact.cpp:1101
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:230
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BiLinearForm
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
Range
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
SetUpSchur
[Push operators to pipeline]
Definition: plastic.cpp:762
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
MoFEM::SmartPetscObj< IS >
SetUpSchurImpl::P
SmartPetscObj< Mat > P
Definition: SetUpSchurImpl.cpp:29
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
SetUpSchurImpl::aoUp
SmartPetscObj< AO > aoUp
Definition: SetUpSchurImpl.cpp:31
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1060
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:104
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127