v0.14.0
SetUpSchurImpl.cpp
Go to the documentation of this file.
1 /** @file SetUpSchurImpl
2  * @brief
3  * @date 2023-05-13
4  *
5  * @license{This project is released under the MIT License.}
6  *
7  */
8 
9 struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
10 
11  SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj<Mat> m,
12  SmartPetscObj<Mat> p, EshelbianCore *ep_core_ptr)
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  }
21  virtual ~SetUpSchurImpl() { S.reset(); }
22 
23  MoFEMErrorCode setUp(KSP solver);
26 
27 private:
28  SmartPetscObj<Mat> M;
29  SmartPetscObj<Mat> P;
30  SmartPetscObj<Mat> S;
31  SmartPetscObj<AO> aoUp;
32  SmartPetscObj<IS> isA00;
33  SmartPetscObj<IS> isShift;
34 
37 };
38 
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 }
201 
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 }
212 
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 }
228 
229 boost::shared_ptr<EshelbianCore::SetUpSchur>
230 EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
231  SmartPetscObj<Mat> m,
232  SmartPetscObj<Mat> p,
233  EshelbianCore *ep_core_ptr) {
234  return boost::shared_ptr<SetUpSchur>(
235  new SetUpSchurImpl(m_field, m, p, ep_core_ptr));
236 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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:318
SetUpSchurImpl
Definition: SetUpSchurImpl.cpp:9
SetUpSchurImpl::SetUpSchurImpl
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
Definition: SetUpSchurImpl.cpp:11
SetUpSchurImpl::setUp
MoFEMErrorCode setUp(KSP solver)
Definition: SetUpSchurImpl.cpp:39
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
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:456
SetUpSchurImpl::M
SmartPetscObj< Mat > M
Definition: SetUpSchurImpl.cpp:28
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
SetUpSchurImpl::~SetUpSchurImpl
virtual ~SetUpSchurImpl()
Definition: SetUpSchurImpl.cpp:21
SetUpSchurImpl::isA00
SmartPetscObj< IS > isA00
Definition: SetUpSchurImpl.cpp:32
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1125
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1037
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
SetUpSchurImpl::isShift
SmartPetscObj< IS > isShift
Definition: SetUpSchurImpl.cpp:33
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
SetUpSchurImpl::preProc
MoFEMErrorCode preProc()
Definition: SetUpSchurImpl.cpp:202
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:215
SetUpSchurImpl::mField
MoFEM::Interface & mField
Definition: SetUpSchurImpl.cpp:35
EshelbianPlasticity::EshelbianCore
enum RotSelector EshelbianCore
Definition: EshelbianPlasticity.cpp:84
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
EshelbianPlasticity::EshelbianCore
Definition: EshelbianPlasticity.hpp:862
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:711
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
SetUpSchurImpl::postProc
MoFEMErrorCode postProc()
Definition: SetUpSchurImpl.cpp:213
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
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:1123
SetUpSchurImpl::P
SmartPetscObj< Mat > P
Definition: SetUpSchurImpl.cpp:29