v0.14.0
Loading...
Searching...
No Matches
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
9struct SetUpSchurImpl : public EshelbianCore::SetUpSchur {
10
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);
24 MoFEMErrorCode preProc();
25 MoFEMErrorCode postProc();
26
27private:
34
36 EshelbianCore *epCorePtr;
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, epCorePtr->bubbleField,
51
52 epCorePtr->rotAxis, epCorePtr->spatialL2Disp,
53
54 epCorePtr->piolaStress
55
56 };
57 std::vector<std::string> field_list{epCorePtr->piolaStress};
58
59 auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
61
62 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
63 CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
64 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
65 CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
66 CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
67 CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->naturalBcElement);
68
69 auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM - 1));
70 std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr};
71
72 if (field_list.size() != dm_range_list.size()) {
73 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
74 "Both ranges should have the same size");
75 }
76
77 int r_idx = 0;
78 for (auto f : field_list) {
79 MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
80 CHKERR DMMoFEMAddSubFieldRow(dm_sub, f.c_str(), dm_range_list[r_idx]);
81 CHKERR DMMoFEMAddSubFieldCol(dm_sub, f.c_str(), dm_range_list[r_idx]);
82 ++r_idx;
83 }
84 CHKERR DMSetUp(dm_sub);
86 };
87
88 auto create_a00_is = [&](SmartPetscObj<IS> &is_a00) {
90 auto vols = get_ents_by_dim(SPACE_DIM);
91 std::vector<SmartPetscObj<IS>> is_vec;
92 std::vector<Range *> range_list_ptr(schur_field_list.size(), &vols);
93 int r_idx = 0;
94 for (auto f : schur_field_list) {
96 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
97 "ELASTIC_PROBLEM", ROW, f, 0, MAX_DOFS_ON_ENTITY, is,
98 range_list_ptr[r_idx]);
99 is_vec.push_back(is);
100 ++r_idx;
101 }
102 if (is_vec.size()) {
103 is_a00 = is_vec[0];
104 for (auto i = 1; i < is_vec.size(); ++i) {
105 IS is_union_raw;
106 CHKERR ISExpand(is_a00, is_vec[i], &is_union_raw);
107 is_a00 = SmartPetscObj<IS>(is_union_raw);
108 }
109 CHKERR ISSort(is_a00);
110 }
112 };
113
114 PC pc;
115 CHKERR KSPSetFromOptions(solver);
116 CHKERR KSPGetPC(solver, &pc);
117 PetscBool is_pcfs = PETSC_FALSE;
118 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
119 if (is_pcfs) {
120
121 SmartPetscObj<DM> schur_dm, a00_dm;
122 CHKERR create_schur_dm(schur_dm);
123 CHKERR create_a00_is(isA00);
124
125 if (S) {
128 "Is expected that schur matrix is not allocated. This is "
129 "possible only is PC is set up twice");
130 }
131
132 S = createDMMatrix(schur_dm);
133
134 auto set_ops = [&]() {
136 auto range_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM));
137 std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
138 range_ptr);
139 std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
141 std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
143 std::vector<bool> symm_list(schur_field_list.size(), false);
144
145 auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
146 aoUp = createAOMappingIS(dm_is, PETSC_NULL);
147 ao_list.back() = aoUp;
148 mat_list.back() = S;
149
150 const bool symmetric_system =
151 EshelbianCore::gradApperoximator <= MODERATE_ROT &&
152 EshelbianCore::rotSelector == SMALL_ROT;
153 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
155 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
156 new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
157 ao_list, mat_list, symm_list,
158 symmetric_system));
159 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
161 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
162 new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
163 ao_list, mat_list, symm_list,
164 symmetric_system));
165
167 };
168
169 auto set_pc = [&]() {
171 CHKERR PCFieldSplitSetIS(pc, NULL, isA00);
172 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
174 };
175
176 CHKERR set_ops();
177 CHKERR set_pc();
178
179 } else {
180 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
182 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
183 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
184 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
186 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
187 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
188 }
189
191}
192
195 MOFEM_LOG("EP", Sev::verbose) << "Zero Schur";
196 if (SetUpSchurImpl::S) {
197 CHKERR MatZeroEntries(S);
198 }
199 CHKERR MatZeroEntries(M);
200 CHKERR MatZeroEntries(P);
202}
203
206 if (S) {
207 MOFEM_LOG("EP", Sev::verbose) << "Assemble Schur";
208 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
209 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
210 }
211 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
212 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
213 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
214 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
215 CHKERR MatAXPY(P, 1, M, SAME_NONZERO_PATTERN);
216
218}
219
220boost::shared_ptr<EshelbianCore::SetUpSchur>
221EshelbianCore::SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field,
222 SmartPetscObj<Mat> m,
223 SmartPetscObj<Mat> p,
224 EshelbianCore *ep_core_ptr) {
225 return boost::shared_ptr<SetUpSchur>(
226 new SetUpSchurImpl(m_field, m, p, ep_core_ptr));
227}
static Index< 'p', 3 > p
@ ROW
Definition: definitions.h:123
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
const int dim
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
Definition: plastic.cpp:762
EshelbianCore * epCorePtr
SmartPetscObj< IS > isA00
SetUpSchurImpl(MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
MoFEMErrorCode setUp(KSP solver)
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SmartPetscObj< Mat > M
virtual ~SetUpSchurImpl()
MoFEMErrorCode postProc()
SmartPetscObj< IS > isShift
SmartPetscObj< AO > aoUp
MoFEMErrorCode preProc()
MoFEM::Interface & mField