v0.16.3
Loading...
Searching...
No Matches
ConstrainMatrixCtx.cpp
Go to the documentation of this file.
1/** \file ConstrainMatrixCtx.cpp
2 * \brief Implementation of constraint projection matrices
3 */
4
5#include <ConstrainMatrixCtx.hpp>
6
7#include <memory>
8
9namespace MoFEM {
10
11namespace {
12
13struct ProjectionEvents {
14 ProjectionEvents() {
15 PetscErrorCode ierr;
16 ierr = PetscLogEventRegister("ProjectionInit", 0, &projectionInit);
17 CHKERRABORT(PETSC_COMM_WORLD, ierr);
18 ierr = PetscLogEventRegister("ProjectionQ", 0, &projectionQ);
19 CHKERRABORT(PETSC_COMM_WORLD, ierr);
20 ierr = PetscLogEventRegister("ProjectionP", 0, &projectionP);
21 CHKERRABORT(PETSC_COMM_WORLD, ierr);
22 ierr = PetscLogEventRegister("ProjectionR", 0, &projectionR);
23 CHKERRABORT(PETSC_COMM_WORLD, ierr);
24 ierr = PetscLogEventRegister("ProjectionRT", 0, &projectionRT);
25 CHKERRABORT(PETSC_COMM_WORLD, ierr);
26 ierr =
27 PetscLogEventRegister("ProjectionCTC_QTKQ", 0, &projectionCTCPlusQTKQ);
28 CHKERRABORT(PETSC_COMM_WORLD, ierr);
29 }
30
31 PetscLogEvent projectionInit;
32 PetscLogEvent projectionQ;
33 PetscLogEvent projectionP;
34 PetscLogEvent projectionR;
35 PetscLogEvent projectionRT;
36 PetscLogEvent projectionCTCPlusQTKQ;
37};
38
39const ProjectionEvents &getProjectionEvents() {
40 static const ProjectionEvents events;
41 return events;
42}
43
44struct LogEventScope {
45 LogEventScope(PetscLogEvent event, MPI_Comm comm)
46 : event(event), comm(comm) {}
47
48 ~LogEventScope() {
49 if (active) {
50 const auto ierr = PetscLogEventEnd(event, 0, 0, 0, 0);
51 CHKERRABORT(comm, ierr);
52 }
53 }
54
55 MoFEMErrorCode end() {
56 MoFEMFunctionBeginHot;
57 active = false;
58 CHKERR PetscLogEventEnd(event, 0, 0, 0, 0);
59 MoFEMFunctionReturnHot(0);
60 }
61
62private:
63 PetscLogEvent event;
64 MPI_Comm comm;
65 bool active = true;
66};
67
68using ContextOwner = boost::shared_ptr<ConstrainMatrixCtx>;
69constexpr char inverseCCTConstraintMatrixKey[] =
70 "_MoFEM_ConstrainMatrixCtx_inverse_CCT_constraint_matrix";
71
72struct CCTInverseKSPCtx {
73 CCTInverseKSPCtx(MPI_Comm comm, Mat constrain_matrix, KSP external_ksp,
74 const std::string &options_prefix)
75 : comm(comm), C(constrain_matrix, true), ksp(external_ksp, true),
76 optionsPrefix(options_prefix), configureFromOptions(!external_ksp) {}
77
78 MoFEMErrorCode update(Mat inverse_mat);
79 MoFEMErrorCode apply(Mat inverse_mat, Vec x, Vec f, bool transpose);
80
81 MPI_Comm comm;
82 SmartPetscObj<Mat> C;
83 SmartPetscObj<Mat> CT;
84 SmartPetscObj<Mat> CCT;
85 SmartPetscObj<KSP> ksp;
86 std::string optionsPrefix;
87 PetscObjectState matrixState = 0;
88 PetscObjectState nonzeroState = 0;
89 bool configureFromOptions;
90 bool initialized = false;
91};
92
93enum class ProjectionOperator { Q, P, R, RT, CTCPlusQTKQ };
94
95MoFEMErrorCode projectionMatrixMultQ(Mat mat, Vec x, Vec f);
96MoFEMErrorCode projectionMatrixMultP(Mat mat, Vec x, Vec f);
97MoFEMErrorCode projectionMatrixMultR(Mat mat, Vec x, Vec f);
98MoFEMErrorCode projectionMatrixMultRT(Mat mat, Vec x, Vec f);
99MoFEMErrorCode projectionMatrixMultCTCPlusQTKQ(Mat mat, Vec x, Vec f);
100MoFEMErrorCode cctInverseKSPMult(Mat mat, Vec x, Vec f);
101MoFEMErrorCode cctInverseKSPMultTranspose(Mat mat, Vec x, Vec f);
102
103PetscErrorCode destroyProjectionMatrixContext(void *void_ctx) {
104 delete static_cast<ContextOwner *>(void_ctx);
105 return 0;
106}
107
108PetscErrorCode destroyCCTInverseKSPContext(void *void_ctx) {
109 delete static_cast<CCTInverseKSPCtx *>(void_ctx);
110 return 0;
111}
112
113} // namespace
114
115struct ConstrainMatrixCtxImpl : public ConstrainMatrixCtx {
116
117 ConstrainMatrixCtxImpl(MPI_Comm comm, DM x_dm, DM y_dm);
118 ConstrainMatrixCtxImpl(MPI_Comm comm, VecScatter scatter);
119
120 ~ConstrainMatrixCtxImpl() override = default;
121
122 MoFEMErrorCode setConstrainMatrix(Mat constrain_matrix,
123 Mat inverse_cct) override;
124 Mat getConstrainMatrix() const override;
125 Mat getInverseCCT() const override;
126 MoFEMErrorCode setSystemMatrix(Mat system_matrix) override;
127 Mat getSystemMatrix() const override;
128
129 MoFEMErrorCode initializeQorP(Vec x) override;
130 MoFEMErrorCode initializeQTKQ() override;
131 MoFEMErrorCode recalculateCT() override;
132 MoFEMErrorCode recalculateCTC() override;
133
134 MoFEMErrorCode applyQ(Vec x, Vec f) override;
135 MoFEMErrorCode applyP(Vec x, Vec f) override;
136 MoFEMErrorCode applyR(Vec x, Vec f) override;
137 MoFEMErrorCode applyRT(Vec x, Vec f) override;
138 MoFEMErrorCode applyCTC_QTKQ(Vec x, Vec f) override;
139
140private:
141 MoFEMErrorCode invalidateProjectionData();
142 MoFEMErrorCode invalidateSystemData();
143 MoFEMErrorCode validatePetscObject(PetscObject object,
144 const char *object_name) const;
145 MoFEMErrorCode validateVectorLayout(Vec vec, PetscInt local_size,
146 PetscInt global_size,
147 const char *vector_name) const;
148 MoFEMErrorCode validateInverseCCTLayout(Mat constrain_matrix,
149 Mat inverse_cct) const;
150 MoFEMErrorCode validateSystemMatrixLayout() const;
151 MoFEMErrorCode applyInverseCCT(Vec x, Vec f, bool transpose);
152 MoFEMErrorCode applyQImpl(Vec x, Vec f, bool transpose_inverse);
153
154 MPI_Comm comm;
155 SmartPetscObj<DM> xDm;
156 SmartPetscObj<DM> yDm;
157
158 SmartPetscObj<Mat> C;
159 SmartPetscObj<Mat> CT;
160 SmartPetscObj<Mat> inverseCCT;
161 SmartPetscObj<Mat> CTC;
162 SmartPetscObj<Mat> K;
163 SmartPetscObj<Vec> Cx;
164 SmartPetscObj<Vec> CCTm1Cx;
165 SmartPetscObj<Vec> CTCCTm1Cx;
166 SmartPetscObj<Vec> CTCx;
167 SmartPetscObj<Vec> X;
168 SmartPetscObj<Vec> Qx;
169 SmartPetscObj<Vec> KQx;
170 SmartPetscObj<VecScatter> scatter;
171
172 bool projectionInitialized = false;
173 bool qtkqInitialized = false;
174 bool buildScatterFromDms = true;
175
176 PetscInt fullLocalSize = 0;
177 PetscInt fullGlobalSize = 0;
178};
179
180ConstrainMatrixCtxImpl::ConstrainMatrixCtxImpl(MPI_Comm comm, DM x_dm, DM y_dm)
181 : comm(comm), xDm(x_dm, true), yDm(y_dm, true) {
182 getProjectionEvents();
183}
184
185ConstrainMatrixCtxImpl::ConstrainMatrixCtxImpl(MPI_Comm comm,
186 VecScatter scatter)
187 : comm(comm), scatter(scatter, true), buildScatterFromDms(false) {
188 getProjectionEvents();
189}
190
191MoFEMErrorCode
192ConstrainMatrixCtxImpl::validatePetscObject(PetscObject object,
193 const char *object_name) const {
194 MoFEMFunctionBeginHot;
195 if (!object) {
196 MOFEM_LOG("SELF", Sev::error)
197 << "Null PETSc object supplied for " << object_name;
198 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY, "Null PETSc object");
199 }
200
201 MPI_Comm object_comm;
202 CHKERR PetscObjectGetComm(object, &object_comm);
203 int relation;
204 CHKERR MPI_Comm_compare(comm, object_comm, &relation);
205 if (relation > MPI_CONGRUENT) {
206 MOFEM_LOG("SELF", Sev::error)
207 << object_name << " uses a communicator incompatible with the context";
208 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
209 "PETSc object uses an incompatible communicator");
210 }
211 MoFEMFunctionReturnHot(0);
212}
213
214MoFEMErrorCode
215ConstrainMatrixCtxImpl::validateVectorLayout(Vec vec, PetscInt local_size,
216 PetscInt global_size,
217 const char *vector_name) const {
218 MoFEMFunctionBeginHot;
219 CHKERR validatePetscObject(getPetscObject(vec), vector_name);
220 PetscInt actual_local_size;
221 PetscInt actual_global_size;
222 CHKERR VecGetLocalSize(vec, &actual_local_size);
223 CHKERR VecGetSize(vec, &actual_global_size);
224 if (actual_local_size != local_size || actual_global_size != global_size) {
225 MOFEM_LOG("SELF", Sev::error)
226 << vector_name << " layout is (" << actual_local_size << ", "
227 << actual_global_size << "), expected (" << local_size << ", "
228 << global_size << ")";
229 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
230 "Vector layout is incompatible with the projection context");
231 }
232 MoFEMFunctionReturnHot(0);
233}
234
235MoFEMErrorCode
236ConstrainMatrixCtxImpl::validateInverseCCTLayout(Mat constrain_matrix,
237 Mat inverse_cct) const {
238 MoFEMFunctionBegin;
239 CHKERR validatePetscObject(getPetscObject(constrain_matrix),
240 "constraint matrix C");
241 CHKERR validatePetscObject(getPetscObject(inverse_cct),
242 "inverse C C^T action");
243
244 PetscInt c_local_rows;
245 PetscInt c_local_columns;
246 PetscInt c_global_rows;
247 PetscInt c_global_columns;
248 PetscInt inverse_local_rows;
249 PetscInt inverse_local_columns;
250 PetscInt inverse_global_rows;
251 PetscInt inverse_global_columns;
252 CHKERR MatGetLocalSize(constrain_matrix, &c_local_rows, &c_local_columns);
253 CHKERR MatGetSize(constrain_matrix, &c_global_rows, &c_global_columns);
254 CHKERR MatGetLocalSize(inverse_cct, &inverse_local_rows,
255 &inverse_local_columns);
256 CHKERR MatGetSize(inverse_cct, &inverse_global_rows, &inverse_global_columns);
257 if (inverse_local_rows != c_local_rows ||
258 inverse_local_columns != c_local_rows ||
259 inverse_global_rows != c_global_rows ||
260 inverse_global_columns != c_global_rows) {
261 MOFEM_LOG("SELF", Sev::error)
262 << "Inverse C C^T layout, local (" << inverse_local_rows << ", "
263 << inverse_local_columns << "), global (" << inverse_global_rows << ", "
264 << inverse_global_columns
265 << "), does not match the row layout of C, local " << c_local_rows
266 << ", global " << c_global_rows;
267 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
268 "Inverse C C^T action has an incompatible layout");
269 }
270
271 PetscBool has_mult;
272 CHKERR MatHasOperation(inverse_cct, MATOP_MULT, &has_mult);
273 if (!has_mult)
274 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
275 "Inverse C C^T action does not implement MatMult");
276
277 PetscObject associated_c = nullptr;
278 CHKERR PetscObjectQuery(getPetscObject(inverse_cct),
279 inverseCCTConstraintMatrixKey, &associated_c);
280 if (associated_c && associated_c != getPetscObject(constrain_matrix)) {
281 MOFEM_LOG("SELF", Sev::error)
282 << "KSP-backed inverse C C^T action was created for a different "
283 "constraint matrix";
284 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
285 "Inverse C C^T action is associated with a different C");
286 }
287 MoFEMFunctionReturn(0);
288}
289
290MoFEMErrorCode ConstrainMatrixCtxImpl::validateSystemMatrixLayout() const {
291 MoFEMFunctionBegin;
292 CHKERR validatePetscObject(getPetscObject(K.get()), "system matrix K");
293 PetscInt local_rows;
294 PetscInt local_columns;
295 PetscInt global_rows;
296 PetscInt global_columns;
297 CHKERR MatGetLocalSize(K.get(), &local_rows, &local_columns);
298 CHKERR MatGetSize(K.get(), &global_rows, &global_columns);
299 if (local_rows != local_columns || global_rows != global_columns) {
300 MOFEM_LOG("SELF", Sev::error)
301 << "System matrix K must be square, local size (" << local_rows << ", "
302 << local_columns << "), global size (" << global_rows << ", "
303 << global_columns << ")";
304 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
305 "System matrix K must have matching row and column layouts");
306 }
307 if (projectionInitialized &&
308 (local_rows != fullLocalSize || global_rows != fullGlobalSize)) {
309 MOFEM_LOG("SELF", Sev::error)
310 << "System matrix K layout (" << local_rows << ", " << global_rows
311 << ") does not match the full vector layout (" << fullLocalSize << ", "
312 << fullGlobalSize << ")";
313 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
314 "System matrix K is incompatible with the full vector layout");
315 }
316 MoFEMFunctionReturn(0);
317}
318
319MoFEMErrorCode ConstrainMatrixCtxImpl::invalidateSystemData() {
320 MoFEMFunctionBeginHot;
321 qtkqInitialized = false;
322 CTC.reset();
323 CTCx.reset();
324 Qx.reset();
325 KQx.reset();
326 MoFEMFunctionReturnHot(0);
327}
328
329MoFEMErrorCode ConstrainMatrixCtxImpl::invalidateProjectionData() {
330 MoFEMFunctionBeginHot;
331 projectionInitialized = false;
332 CT.reset();
333 Cx.reset();
334 CCTm1Cx.reset();
335 CTCCTm1Cx.reset();
336 X.reset();
337 fullLocalSize = 0;
338 fullGlobalSize = 0;
339 if (buildScatterFromDms)
340 scatter.reset();
341 CHKERR invalidateSystemData();
342 MoFEMFunctionReturnHot(0);
343}
344
345MoFEMErrorCode ConstrainMatrixCtxImpl::setConstrainMatrix(Mat constrain_matrix,
346 Mat inverse_cct) {
347 MoFEMFunctionBegin;
348 CHKERR validateInverseCCTLayout(constrain_matrix, inverse_cct);
349 C = SmartPetscObj<Mat>(constrain_matrix, true);
350 inverseCCT = SmartPetscObj<Mat>(inverse_cct, true);
351 CHKERR invalidateProjectionData();
352 MoFEMFunctionReturn(0);
353}
354
355Mat ConstrainMatrixCtxImpl::getConstrainMatrix() const { return C.get(); }
356
357Mat ConstrainMatrixCtxImpl::getInverseCCT() const { return inverseCCT.get(); }
358
359MoFEMErrorCode ConstrainMatrixCtxImpl::setSystemMatrix(Mat system_matrix) {
360 MoFEMFunctionBegin;
361 CHKERR validatePetscObject(getPetscObject(system_matrix), "system matrix K");
362 K = SmartPetscObj<Mat>(system_matrix, true);
363 CHKERR invalidateSystemData();
364 MoFEMFunctionReturn(0);
365}
366
367Mat ConstrainMatrixCtxImpl::getSystemMatrix() const { return K.get(); }
368
369MoFEMErrorCode ConstrainMatrixCtxImpl::initializeQorP(Vec x) {
370 MoFEMFunctionBegin;
371 if (projectionInitialized)
372 MoFEMFunctionReturnHot(0);
373
374 CHKERR validateInverseCCTLayout(C, inverseCCT);
375 CHKERR validatePetscObject(getPetscObject(x), "full vector x");
376 if (buildScatterFromDms) {
377 CHKERR validatePetscObject(getPetscObject(xDm.get()), "full problem DM");
378 CHKERR validatePetscObject(getPetscObject(yDm.get()), "constraint DM");
379 } else {
380 CHKERR validatePetscObject(getPetscObject(scatter.get()),
381 "constraint scatter");
382 }
383
384 const auto &events = getProjectionEvents();
385 CHKERR PetscLogEventBegin(events.projectionInit, 0, 0, 0, 0);
386 LogEventScope event_scope(events.projectionInit, comm);
387
388 Mat raw_ct;
389 CHKERR MatTranspose(C, MAT_INITIAL_MATRIX, &raw_ct);
390 SmartPetscObj<Mat> new_ct(raw_ct);
391
392 Vec raw_x;
393 Vec raw_cx;
394 CHKERR MatCreateVecs(C, &raw_x, &raw_cx);
395 SmartPetscObj<Vec> new_x(raw_x);
396 SmartPetscObj<Vec> new_cx(raw_cx);
397
398 Vec raw_cctm1_cx;
399 CHKERR VecDuplicate(new_cx, &raw_cctm1_cx);
400 SmartPetscObj<Vec> new_cctm1_cx(raw_cctm1_cx);
401
402 Vec raw_ct_cctm1_cx;
403 CHKERR VecDuplicate(new_x, &raw_ct_cctm1_cx);
404 SmartPetscObj<Vec> new_ct_cctm1_cx(raw_ct_cctm1_cx);
405
406 SmartPetscObj<VecScatter> new_scatter = scatter;
407 if (buildScatterFromDms) {
408 MoFEM::Interface *x_m_field_ptr;
409 MoFEM::Interface *y_m_field_ptr;
410 CHKERR DMoFEMGetInterfacePtr(xDm, &x_m_field_ptr);
411 CHKERR DMoFEMGetInterfacePtr(yDm, &y_m_field_ptr);
412 if (x_m_field_ptr != y_m_field_ptr) {
413 MOFEM_LOG("SELF", Sev::error)
414 << "The full problem DM and constraint DM belong to different MoFEM "
415 "interfaces";
416 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
417 "DMs must belong to the same MoFEM interface");
418 }
419
420 const MoFEM::Problem *x_problem_ptr;
421 const MoFEM::Problem *y_problem_ptr;
422 CHKERR DMMoFEMGetProblemPtr(xDm, &x_problem_ptr);
423 CHKERR DMMoFEMGetProblemPtr(yDm, &y_problem_ptr);
424 CHKERR x_m_field_ptr->getInterface<VecManager>()->vecScatterCreate(
425 x, x_problem_ptr->getName(), ROW, new_x, y_problem_ptr->getName(), COL,
426 new_scatter);
427 }
428
429 PetscInt new_full_local_size;
430 PetscInt new_full_global_size;
431 CHKERR VecGetLocalSize(x, &new_full_local_size);
432 CHKERR VecGetSize(x, &new_full_global_size);
433
434 PetscInt c_rows;
435 PetscInt c_columns;
436 PetscInt inverse_cct_rows;
437 PetscInt inverse_cct_columns;
438 CHKERR MatGetSize(C, &c_rows, &c_columns);
439 CHKERR MatGetSize(inverseCCT, &inverse_cct_rows, &inverse_cct_columns);
440 MOFEM_LOG("WORLD", Sev::verbose)
441 << "Constraint projection initialized: C size (" << c_rows << ", "
442 << c_columns << "), inverse C C^T action size (" << inverse_cct_rows
443 << ", " << inverse_cct_columns << "), full vector size "
444 << new_full_global_size;
445
446 CT = new_ct;
447 X = new_x;
448 Cx = new_cx;
449 CCTm1Cx = new_cctm1_cx;
450 CTCCTm1Cx = new_ct_cctm1_cx;
451 scatter = new_scatter;
452 fullLocalSize = new_full_local_size;
453 fullGlobalSize = new_full_global_size;
454 projectionInitialized = true;
455
456 CHKERR event_scope.end();
457 MoFEMFunctionReturn(0);
458}
459
460MoFEMErrorCode ConstrainMatrixCtxImpl::recalculateCT() {
461 MoFEMFunctionBegin;
462 if (!projectionInitialized)
463 MoFEMFunctionReturnHot(0);
464
465 CHKERR validateInverseCCTLayout(C, inverseCCT);
466 Mat raw_mat = CT.get();
467 CHKERR MatTranspose(C, MAT_REUSE_MATRIX, &raw_mat);
468 MOFEM_LOG("WORLD", Sev::verbose) << "Recalculated C^T";
469 MoFEMFunctionReturn(0);
470}
471
472MoFEMErrorCode ConstrainMatrixCtxImpl::initializeQTKQ() {
473 MoFEMFunctionBegin;
474 if (qtkqInitialized)
475 MoFEMFunctionReturnHot(0);
476
477 CHKERR validatePetscObject(getPetscObject(C.get()), "constraint matrix C");
478 CHKERR validateSystemMatrixLayout();
479
480 const auto &events = getProjectionEvents();
481 CHKERR PetscLogEventBegin(events.projectionInit, 0, 0, 0, 0);
482 LogEventScope event_scope(events.projectionInit, comm);
483
484 Mat raw_ctc;
485 CHKERR MatTransposeMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &raw_ctc);
486 SmartPetscObj<Mat> new_ctc(raw_ctc);
487
488 Vec raw_qx;
489 Vec raw_kqx;
490 CHKERR MatCreateVecs(K, &raw_qx, &raw_kqx);
491 SmartPetscObj<Vec> new_qx(raw_qx);
492 SmartPetscObj<Vec> new_kqx(raw_kqx);
493
494 Vec raw_ctcx;
495 CHKERR MatCreateVecs(new_ctc, nullptr, &raw_ctcx);
496 SmartPetscObj<Vec> new_ctcx(raw_ctcx);
497
498 PetscInt ctc_rows;
499 PetscInt ctc_columns;
500 CHKERR MatGetSize(new_ctc, &ctc_rows, &ctc_columns);
501 MOFEM_LOG("WORLD", Sev::verbose)
502 << "Initialized C^T C + Q^T K Q work data: C^T C size (" << ctc_rows
503 << ", " << ctc_columns << ")";
504
505 CTC = new_ctc;
506 Qx = new_qx;
507 KQx = new_kqx;
508 CTCx = new_ctcx;
509 qtkqInitialized = true;
510
511 CHKERR event_scope.end();
512 MoFEMFunctionReturn(0);
513}
514
515MoFEMErrorCode ConstrainMatrixCtxImpl::recalculateCTC() {
516 MoFEMFunctionBegin;
517 if (!qtkqInitialized)
518 MoFEMFunctionReturnHot(0);
519 Mat raw_mat = CTC.get();
520 CHKERR MatTransposeMatMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &raw_mat);
521 MOFEM_LOG("WORLD", Sev::verbose) << "Recalculated C^T C";
522 MoFEMFunctionReturn(0);
523}
524
525MoFEMErrorCode ConstrainMatrixCtxImpl::applyInverseCCT(Vec x, Vec f,
526 bool transpose) {
527 MoFEMFunctionBeginHot;
528 if (transpose) {
529 PetscBool has_mult_transpose;
530 CHKERR MatHasOperation(inverseCCT, MATOP_MULT_TRANSPOSE,
531 &has_mult_transpose);
532 if (!has_mult_transpose) {
533 MOFEM_LOG("SELF", Sev::error)
534 << "Inverse C C^T action does not implement MatMultTranspose";
535 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
536 "Inverse C C^T action does not implement MatMultTranspose");
537 }
538 CHKERR MatMultTranspose(inverseCCT, x, f);
539 } else {
540 CHKERR MatMult(inverseCCT, x, f);
541 }
542 MoFEMFunctionReturnHot(0);
543}
544
545MoFEMErrorCode ConstrainMatrixCtxImpl::applyQImpl(Vec x, Vec f,
546 bool transpose_inverse) {
547 MoFEMFunctionBegin;
548 CHKERR initializeQorP(x);
549 CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "Q input");
550 CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "Q output");
551 CHKERR VecCopy(x, f);
552 CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
553 CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
554 CHKERR MatMult(C, X, Cx);
555 CHKERR applyInverseCCT(Cx, CCTm1Cx, transpose_inverse);
556 CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
557 CHKERR VecScale(CTCCTm1Cx, -1);
558 CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, ADD_VALUES, SCATTER_REVERSE);
559 CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, ADD_VALUES, SCATTER_REVERSE);
560 MoFEMFunctionReturn(0);
561}
562
563MoFEMErrorCode ConstrainMatrixCtxImpl::applyQ(Vec x, Vec f) {
564 MoFEMFunctionBegin;
565 const auto &events = getProjectionEvents();
566 CHKERR PetscLogEventBegin(events.projectionQ, 0, 0, 0, 0);
567 LogEventScope event_scope(events.projectionQ, comm);
568
569 CHKERR applyQImpl(x, f, false);
570
571 CHKERR event_scope.end();
572 MoFEMFunctionReturn(0);
573}
574
575MoFEMErrorCode ConstrainMatrixCtxImpl::applyP(Vec x, Vec f) {
576 MoFEMFunctionBegin;
577 const auto &events = getProjectionEvents();
578 CHKERR PetscLogEventBegin(events.projectionP, 0, 0, 0, 0);
579 LogEventScope event_scope(events.projectionP, comm);
580
581 CHKERR initializeQorP(x);
582 CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "P input");
583 CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "P output");
584 CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
585 CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
586 CHKERR MatMult(C, X, Cx);
587 CHKERR applyInverseCCT(Cx, CCTm1Cx, false);
588 CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
589 CHKERR VecZeroEntries(f);
590 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
591 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
592 CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
593 CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
594
595 CHKERR event_scope.end();
596 MoFEMFunctionReturn(0);
597}
598
599MoFEMErrorCode ConstrainMatrixCtxImpl::applyR(Vec x, Vec f) {
600 MoFEMFunctionBegin;
601 const auto &events = getProjectionEvents();
602 CHKERR PetscLogEventBegin(events.projectionR, 0, 0, 0, 0);
603 LogEventScope event_scope(events.projectionR, comm);
604
605 if (!projectionInitialized) {
606 MOFEM_LOG("SELF", Sev::error)
607 << "R requires initializeQorP() or a previously created Q/P matrix";
608 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
609 "Projection context is not initialized");
610 }
611 PetscInt constraint_local_size;
612 PetscInt constraint_global_size;
613 CHKERR VecGetLocalSize(Cx, &constraint_local_size);
614 CHKERR VecGetSize(Cx, &constraint_global_size);
615 CHKERR validateVectorLayout(x, constraint_local_size, constraint_global_size,
616 "R input");
617 CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize, "R output");
618 CHKERR applyInverseCCT(x, CCTm1Cx, false);
619 CHKERR MatMult(CT, CCTm1Cx, CTCCTm1Cx);
620 CHKERR VecZeroEntries(f);
621 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
622 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
623 CHKERR VecScatterBegin(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
624 CHKERR VecScatterEnd(scatter, CTCCTm1Cx, f, INSERT_VALUES, SCATTER_REVERSE);
625
626 CHKERR event_scope.end();
627 MoFEMFunctionReturn(0);
628}
629
630MoFEMErrorCode ConstrainMatrixCtxImpl::applyRT(Vec x, Vec f) {
631 MoFEMFunctionBegin;
632 const auto &events = getProjectionEvents();
633 CHKERR PetscLogEventBegin(events.projectionRT, 0, 0, 0, 0);
634 LogEventScope event_scope(events.projectionRT, comm);
635
636 if (!projectionInitialized) {
637 MOFEM_LOG("SELF", Sev::error)
638 << "R^T requires initializeQorP() or a previously created Q/P matrix";
639 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY,
640 "Projection context is not initialized");
641 }
642 PetscInt constraint_local_size;
643 PetscInt constraint_global_size;
644 CHKERR VecGetLocalSize(Cx, &constraint_local_size);
645 CHKERR VecGetSize(Cx, &constraint_global_size);
646 CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize, "R^T input");
647 CHKERR validateVectorLayout(f, constraint_local_size, constraint_global_size,
648 "R^T output");
649 CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
650 CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
651 CHKERR MatMult(C, X, Cx);
652 CHKERR applyInverseCCT(Cx, f, true);
653
654 CHKERR event_scope.end();
655 MoFEMFunctionReturn(0);
656}
657
658MoFEMErrorCode ConstrainMatrixCtxImpl::applyCTC_QTKQ(Vec x, Vec f) {
659 MoFEMFunctionBegin;
660 const auto &events = getProjectionEvents();
661 CHKERR PetscLogEventBegin(events.projectionCTCPlusQTKQ, 0, 0, 0, 0);
662 LogEventScope event_scope(events.projectionCTCPlusQTKQ, comm);
663
664 CHKERR initializeQorP(x);
665 CHKERR initializeQTKQ();
666 CHKERR validateSystemMatrixLayout();
667 CHKERR validateVectorLayout(x, fullLocalSize, fullGlobalSize,
668 "C^T C + Q^T K Q input");
669 CHKERR validateVectorLayout(f, fullLocalSize, fullGlobalSize,
670 "C^T C + Q^T K Q output");
671
672 // Consume x into CTCx before Qx/f below can overwrite it, so that x and f
673 // may safely alias the same Vec.
674 CHKERR VecScatterBegin(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
675 CHKERR VecScatterEnd(scatter, x, X, INSERT_VALUES, SCATTER_FORWARD);
676 CHKERR MatMult(CTC, X, CTCx);
677
678 CHKERR applyQImpl(x, Qx, false);
679 CHKERR MatMult(K, Qx, KQx);
680 CHKERR applyQImpl(KQx, f, true);
681 CHKERR VecScatterBegin(scatter, CTCx, f, ADD_VALUES, SCATTER_REVERSE);
682 CHKERR VecScatterEnd(scatter, CTCx, f, ADD_VALUES, SCATTER_REVERSE);
683
684 CHKERR event_scope.end();
685 MoFEMFunctionReturn(0);
686}
687
688namespace {
689
690MoFEMErrorCode CCTInverseKSPCtx::update(Mat inverse_mat) {
691 MoFEMFunctionBegin;
692 PetscObjectState current_matrix_state;
693 PetscObjectState current_nonzero_state;
694 CHKERR PetscObjectStateGet(getPetscObject(C.get()), &current_matrix_state);
695 CHKERR MatGetNonzeroState(C, &current_nonzero_state);
696
697 PetscBool local_values_changed =
698 (!initialized || current_matrix_state != matrixState) ? PETSC_TRUE
699 : PETSC_FALSE;
700 PetscBool local_structure_changed =
701 (!initialized || current_nonzero_state != nonzeroState) ? PETSC_TRUE
702 : PETSC_FALSE;
703 PetscBool values_changed;
704 PetscBool structure_changed;
705 CHKERR MPIU_Allreduce(&local_values_changed, &values_changed, 1, MPIU_BOOL,
706 MPI_LOR, comm);
707 CHKERR MPIU_Allreduce(&local_structure_changed, &structure_changed, 1,
708 MPIU_BOOL, MPI_LOR, comm);
709 if (!values_changed)
710 MoFEMFunctionReturnHot(0);
711
712 SmartPetscObj<Mat> updated_ct = CT;
713 SmartPetscObj<Mat> updated_cct = CCT;
714 if (structure_changed) {
715 Mat raw_ct;
716 CHKERR MatTranspose(C, MAT_INITIAL_MATRIX, &raw_ct);
717 updated_ct = SmartPetscObj<Mat>(raw_ct);
718
719 Mat raw_cct;
720 CHKERR MatTransposeMatMult(updated_ct, updated_ct, MAT_INITIAL_MATRIX,
721 PETSC_DEFAULT, &raw_cct);
722 updated_cct = SmartPetscObj<Mat>(raw_cct);
723 } else {
724 Mat raw_mat = updated_ct.get();
725 CHKERR MatTranspose(C, MAT_REUSE_MATRIX, &raw_mat);
726 raw_mat = updated_cct.get();
727 CHKERR MatTransposeMatMult(updated_ct, updated_ct, MAT_REUSE_MATRIX,
728 PETSC_DEFAULT, &raw_mat);
729 }
730
731 SmartPetscObj<KSP> updated_ksp = ksp;
732 const bool create_solver = !updated_ksp;
733 if (create_solver)
734 updated_ksp = MoFEM::createKSP(comm);
735 CHKERR KSPSetOperators(updated_ksp, updated_cct, updated_cct);
736 CHKERR KSPSetReusePreconditioner(updated_ksp, PETSC_FALSE);
737 if (create_solver && configureFromOptions) {
738 CHKERR KSPSetOptionsPrefix(updated_ksp, optionsPrefix.c_str());
739 CHKERR KSPSetFromOptions(updated_ksp);
740 }
741 CHKERR KSPSetUp(updated_ksp);
742
743 CT = updated_ct;
744 CCT = updated_cct;
745 ksp = updated_ksp;
746 matrixState = current_matrix_state;
747 nonzeroState = current_nonzero_state;
748 initialized = true;
749 CHKERR PetscObjectStateIncrease(getPetscObject(inverse_mat));
750
751 PetscInt rows;
752 PetscInt columns;
753 CHKERR MatGetSize(CCT, &rows, &columns);
754 MOFEM_LOG("WORLD", Sev::verbose)
755 << "Updated KSP-backed inverse C C^T action, operator size (" << rows
756 << ", " << columns << ")";
757 MoFEMFunctionReturn(0);
758}
759
760MoFEMErrorCode CCTInverseKSPCtx::apply(Mat inverse_mat, Vec x, Vec f,
761 bool transpose) {
762 MoFEMFunctionBegin;
763 PetscInt constraint_size;
764 CHKERR VecGetSize(x, &constraint_size);
765 if (!constraint_size) {
766 CHKERR VecZeroEntries(f);
767 MoFEMFunctionReturnHot(0);
768 }
769
770 CHKERR update(inverse_mat);
771 CHKERR KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
772 CHKERR KSPSetInitialGuessKnoll(ksp, PETSC_FALSE);
773 if (transpose)
774 CHKERR KSPSolveTranspose(ksp, x, f);
775 else
776 CHKERR KSPSolve(ksp, x, f);
777
778 KSPConvergedReason reason;
779 PetscInt iterations;
780 PetscReal residual_norm;
781 CHKERR KSPGetConvergedReason(ksp, &reason);
782 CHKERR KSPGetIterationNumber(ksp, &iterations);
783 CHKERR KSPGetResidualNorm(ksp, &residual_norm);
784 MOFEM_LOG("WORLD", Sev::verbose)
785 << "KSP-backed inverse C C^T action: reason "
786 << KSPConvergedReasons[reason] << ", iterations " << iterations
787 << ", residual " << residual_norm;
788 if (reason <= 0) {
789 MOFEM_LOG("SELF", Sev::error) << "KSP-backed inverse C C^T action failed: "
790 << KSPConvergedReasons[reason];
791 SETERRQ(comm, MOFEM_OPERATION_UNSUCCESSFUL,
792 "KSP-backed inverse C C^T action did not converge");
793 }
794 MoFEMFunctionReturn(0);
795}
796
797MoFEMErrorCode getCCTInverseKSPContext(Mat mat, CCTInverseKSPCtx **ctx_ptr) {
798 MoFEMFunctionBeginHot;
799 void *void_ctx = nullptr;
800 CHKERR MatShellGetContext(mat, &void_ctx);
801 auto *inverse_ctx_ptr = static_cast<CCTInverseKSPCtx *>(void_ctx);
802 if (!inverse_ctx_ptr)
803 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
804 "KSP-backed inverse C C^T shell has no context");
805 *ctx_ptr = inverse_ctx_ptr;
806 MoFEMFunctionReturnHot(0);
807}
808
809MoFEMErrorCode cctInverseKSPMult(Mat mat, Vec x, Vec f) {
810 MoFEMFunctionBegin;
811 CCTInverseKSPCtx *ctx_ptr;
812 CHKERR getCCTInverseKSPContext(mat, &ctx_ptr);
813 CHKERR ctx_ptr->apply(mat, x, f, false);
814 MoFEMFunctionReturn(0);
815}
816
817MoFEMErrorCode cctInverseKSPMultTranspose(Mat mat, Vec x, Vec f) {
818 MoFEMFunctionBegin;
819 CCTInverseKSPCtx *ctx_ptr;
820 CHKERR getCCTInverseKSPContext(mat, &ctx_ptr);
821 CHKERR ctx_ptr->apply(mat, x, f, true);
822 MoFEMFunctionReturn(0);
823}
824
825MoFEMErrorCode getProjectionContext(Mat mat, ConstrainMatrixCtx **ctx_ptr) {
826 MoFEMFunctionBeginHot;
827 void *void_ctx = nullptr;
828 CHKERR MatShellGetContext(mat, &void_ctx);
829 auto *owner_ptr = static_cast<ContextOwner *>(void_ctx);
830 if (!owner_ptr || !*owner_ptr)
831 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
832 "Projection shell matrix has no context");
833 *ctx_ptr = owner_ptr->get();
834 MoFEMFunctionReturnHot(0);
835}
836
837MoFEMErrorCode projectionMatrixMultQ(Mat mat, Vec x, Vec f) {
838 MoFEMFunctionBegin;
839 ConstrainMatrixCtx *ctx_ptr;
840 CHKERR getProjectionContext(mat, &ctx_ptr);
841 CHKERR ctx_ptr->applyQ(x, f);
842 MoFEMFunctionReturn(0);
843}
844
845MoFEMErrorCode projectionMatrixMultP(Mat mat, Vec x, Vec f) {
846 MoFEMFunctionBegin;
847 ConstrainMatrixCtx *ctx_ptr;
848 CHKERR getProjectionContext(mat, &ctx_ptr);
849 CHKERR ctx_ptr->applyP(x, f);
850 MoFEMFunctionReturn(0);
851}
852
853MoFEMErrorCode projectionMatrixMultR(Mat mat, Vec x, Vec f) {
854 MoFEMFunctionBegin;
855 ConstrainMatrixCtx *ctx_ptr;
856 CHKERR getProjectionContext(mat, &ctx_ptr);
857 CHKERR ctx_ptr->applyR(x, f);
858 MoFEMFunctionReturn(0);
859}
860
861MoFEMErrorCode projectionMatrixMultRT(Mat mat, Vec x, Vec f) {
862 MoFEMFunctionBegin;
863 ConstrainMatrixCtx *ctx_ptr;
864 CHKERR getProjectionContext(mat, &ctx_ptr);
865 CHKERR ctx_ptr->applyRT(x, f);
866 MoFEMFunctionReturn(0);
867}
868
869MoFEMErrorCode projectionMatrixMultCTCPlusQTKQ(Mat mat, Vec x, Vec f) {
870 MoFEMFunctionBegin;
871 ConstrainMatrixCtx *ctx_ptr;
872 CHKERR getProjectionContext(mat, &ctx_ptr);
873 CHKERR ctx_ptr->applyCTC_QTKQ(x, f);
874 MoFEMFunctionReturn(0);
875}
876
877MoFEMErrorCode
878createProjectionShellMatrix(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
879 Vec x, ProjectionOperator projection_operator,
880 SmartPetscObj<Mat> &shell_matrix) {
881 MoFEMFunctionBegin;
882 if (!ctx)
883 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
884 "Cannot create a projection matrix from a null context");
885 if (!x)
886 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
887 "Cannot create a projection matrix from a null vector");
888
889 CHKERR ctx->initializeQorP(x);
890 if (projection_operator == ProjectionOperator::CTCPlusQTKQ)
891 CHKERR ctx->initializeQTKQ();
892
893 PetscInt full_local_size;
894 PetscInt full_global_size;
895 CHKERR VecGetLocalSize(x, &full_local_size);
896 CHKERR VecGetSize(x, &full_global_size);
897
898 Mat C = ctx->getConstrainMatrix();
899 PetscInt constraint_local_rows;
900 PetscInt constraint_local_columns;
901 PetscInt constraint_global_rows;
902 PetscInt constraint_global_columns;
903 CHKERR MatGetLocalSize(C, &constraint_local_rows, &constraint_local_columns);
904 CHKERR MatGetSize(C, &constraint_global_rows, &constraint_global_columns);
905
906 PetscInt local_rows = full_local_size;
907 PetscInt local_columns = full_local_size;
908 PetscInt global_rows = full_global_size;
909 PetscInt global_columns = full_global_size;
910 void (*multiply_operation)(void) = nullptr;
911 const char *matrix_name = nullptr;
912
913 switch (projection_operator) {
914 case ProjectionOperator::Q:
915 multiply_operation =
916 reinterpret_cast<void (*)(void)>(projectionMatrixMultQ);
917 matrix_name = "ProjectionQ";
918 break;
919 case ProjectionOperator::P:
920 multiply_operation =
921 reinterpret_cast<void (*)(void)>(projectionMatrixMultP);
922 matrix_name = "ProjectionP";
923 break;
924 case ProjectionOperator::R:
925 local_columns = constraint_local_rows;
926 global_columns = constraint_global_rows;
927 multiply_operation =
928 reinterpret_cast<void (*)(void)>(projectionMatrixMultR);
929 matrix_name = "ProjectionR";
930 break;
931 case ProjectionOperator::RT:
932 local_rows = constraint_local_rows;
933 global_rows = constraint_global_rows;
934 multiply_operation =
935 reinterpret_cast<void (*)(void)>(projectionMatrixMultRT);
936 matrix_name = "ProjectionRT";
937 break;
938 case ProjectionOperator::CTCPlusQTKQ:
939 multiply_operation =
940 reinterpret_cast<void (*)(void)>(projectionMatrixMultCTCPlusQTKQ);
941 matrix_name = "ProjectionCTCPlusQTKQ";
942 break;
943 }
944
945 auto context_owner = std::make_unique<ContextOwner>(ctx);
946 Mat raw_shell_matrix;
947 const auto matrix_comm = getCommFromPetscObject(getPetscObject(x));
948 CHKERR MatCreateShell(matrix_comm, local_rows, local_columns, global_rows,
949 global_columns, context_owner.get(), &raw_shell_matrix);
950 shell_matrix = SmartPetscObj<Mat>(raw_shell_matrix);
951 CHKERR MatShellSetContextDestroy(shell_matrix,
952 destroyProjectionMatrixContext);
953 context_owner.release();
954 CHKERR MatShellSetOperation(shell_matrix, MATOP_MULT, multiply_operation);
955 CHKERR PetscObjectSetName(getPetscObject(shell_matrix.get()), matrix_name);
956
957 MOFEM_LOG("WORLD", Sev::verbose)
958 << "Created " << matrix_name << " shell matrix, global size ("
959 << global_rows << ", " << global_columns << ")";
960 MoFEMFunctionReturn(0);
961}
962
963void validateFactoryCommunicators(MPI_Comm context_comm, PetscObject object,
964 const char *object_name) {
965 const auto object_comm = getCommFromPetscObject(object);
966 int relation;
967 CHK_THROW_MESSAGE(MPI_Comm_compare(context_comm, object_comm, &relation),
968 "Failed to compare PETSc communicators");
969 if (relation > MPI_CONGRUENT)
970 THROW_MESSAGE(std::string(object_name) +
971 " uses a communicator incompatible with the context");
972}
973
974} // namespace
975
976boost::shared_ptr<ConstrainMatrixCtx>
977createConstrainMatrixCtx(DM x_dm, DM y_dm, Mat constrain_matrix,
978 Mat inverse_cct) {
979 if (!x_dm || !y_dm || !constrain_matrix || !inverse_cct)
980 THROW_MESSAGE(
981 "Both DMs, C, and its inverse C C^T action are required to create a "
982 "constraint-matrix context");
983
984 const auto comm = getCommFromPetscObject(getPetscObject(x_dm));
985 validateFactoryCommunicators(comm, getPetscObject(y_dm), "constraint DM");
986 validateFactoryCommunicators(comm, getPetscObject(constrain_matrix),
987 "constraint matrix C");
988 validateFactoryCommunicators(comm, getPetscObject(inverse_cct),
989 "inverse C C^T action");
990
991 MoFEM::Interface *x_m_field_ptr;
992 MoFEM::Interface *y_m_field_ptr;
993 CHK_THROW_MESSAGE(DMoFEMGetInterfacePtr(x_dm, &x_m_field_ptr),
994 "Failed to get the MoFEM interface from the full DM");
995 CHK_THROW_MESSAGE(DMoFEMGetInterfacePtr(y_dm, &y_m_field_ptr),
996 "Failed to get the MoFEM interface from the constraint DM");
997 if (x_m_field_ptr != y_m_field_ptr)
998 THROW_MESSAGE("DMs must belong to the same MoFEM interface");
999
1000 auto ctx = boost::make_shared<ConstrainMatrixCtxImpl>(comm, x_dm, y_dm);
1001 CHK_THROW_MESSAGE(ctx->setConstrainMatrix(constrain_matrix, inverse_cct),
1002 "Set C and its inverse C C^T action");
1003 return ctx;
1004}
1005
1006boost::shared_ptr<ConstrainMatrixCtx>
1007createConstrainMatrixCtx(VecScatter scatter, Mat constrain_matrix,
1008 Mat inverse_cct) {
1009 if (!scatter || !constrain_matrix || !inverse_cct)
1010 THROW_MESSAGE("A scatter, C, and its inverse C C^T action are required to "
1011 "create a constraint-matrix context");
1012 const auto comm = getCommFromPetscObject(getPetscObject(scatter));
1013 validateFactoryCommunicators(comm, getPetscObject(constrain_matrix),
1014 "constraint matrix C");
1015 validateFactoryCommunicators(comm, getPetscObject(inverse_cct),
1016 "inverse C C^T action");
1017 auto ctx = boost::make_shared<ConstrainMatrixCtxImpl>(comm, scatter);
1018 CHK_THROW_MESSAGE(ctx->setConstrainMatrix(constrain_matrix, inverse_cct),
1019 "Set C and its inverse C C^T action");
1020 return ctx;
1021}
1022
1023SmartPetscObj<Mat> createCCTInverseKSP(Mat constrain_matrix, KSP external_ksp,
1024 const std::string &options_prefix) {
1025 if (!constrain_matrix)
1026 THROW_MESSAGE(
1027 "Constraint matrix C is required to create its inverse action");
1028 if (!external_ksp && !options_prefix.empty() && options_prefix.front() == '-')
1029 THROW_MESSAGE("A PETSc options prefix must not start with '-'");
1030
1031 const auto comm = getCommFromPetscObject(getPetscObject(constrain_matrix));
1032 if (external_ksp)
1033 validateFactoryCommunicators(comm, getPetscObject(external_ksp),
1034 "inverse C C^T KSP");
1035
1036 PetscInt local_rows;
1037 PetscInt local_columns;
1038 PetscInt global_rows;
1039 PetscInt global_columns;
1040 CHK_THROW_MESSAGE(
1041 MatGetLocalSize(constrain_matrix, &local_rows, &local_columns),
1042 "Get local size of C");
1043 CHK_THROW_MESSAGE(MatGetSize(constrain_matrix, &global_rows, &global_columns),
1044 "Get global size of C");
1045
1046 auto inverse_ctx = std::make_unique<CCTInverseKSPCtx>(
1047 comm, constrain_matrix, external_ksp, options_prefix);
1048 Mat raw_inverse;
1049 CHK_THROW_MESSAGE(MatCreateShell(comm, local_rows, local_rows, global_rows,
1050 global_rows, inverse_ctx.get(),
1051 &raw_inverse),
1052 "Create KSP-backed inverse C C^T shell");
1053 SmartPetscObj<Mat> inverse(raw_inverse);
1054 CHK_THROW_MESSAGE(PetscObjectCompose(getPetscObject(inverse.get()),
1055 inverseCCTConstraintMatrixKey,
1056 getPetscObject(constrain_matrix)),
1057 "Associate KSP-backed inverse C C^T shell with C");
1058 CHK_THROW_MESSAGE(
1059 MatShellSetContextDestroy(inverse, destroyCCTInverseKSPContext),
1060 "Set KSP-backed inverse C C^T context destroy operation");
1061 inverse_ctx.release();
1062 CHK_THROW_MESSAGE(
1063 MatShellSetOperation(inverse, MATOP_MULT,
1064 reinterpret_cast<void (*)(void)>(cctInverseKSPMult)),
1065 "Set KSP-backed inverse C C^T multiply operation");
1066 CHK_THROW_MESSAGE(
1067 MatShellSetOperation(
1068 inverse, MATOP_MULT_TRANSPOSE,
1069 reinterpret_cast<void (*)(void)>(cctInverseKSPMultTranspose)),
1070 "Set KSP-backed inverse C C^T transpose-multiply operation");
1071 CHK_THROW_MESSAGE(
1072 PetscObjectSetName(getPetscObject(inverse.get()), "CCTInverseKSP"),
1073 "Name KSP-backed inverse C C^T shell");
1074
1075 MOFEM_LOG("WORLD", Sev::verbose)
1076 << "Created KSP-backed inverse C C^T shell, global size (" << global_rows
1077 << ", " << global_rows << ")";
1078 return inverse;
1079}
1080
1081SmartPetscObj<Mat>
1082createProjectionMatrixQ(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
1083 Vec x) {
1084 SmartPetscObj<Mat> matrix;
1085 CHK_THROW_MESSAGE(
1086 createProjectionShellMatrix(ctx, x, ProjectionOperator::Q, matrix),
1087 "Failed to create Q projection matrix");
1088 return matrix;
1089}
1090
1091SmartPetscObj<Mat>
1092createProjectionMatrixP(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
1093 Vec x) {
1094 SmartPetscObj<Mat> matrix;
1095 CHK_THROW_MESSAGE(
1096 createProjectionShellMatrix(ctx, x, ProjectionOperator::P, matrix),
1097 "Failed to create P projection matrix");
1098 return matrix;
1099}
1100
1101SmartPetscObj<Mat>
1102createProjectionMatrixR(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
1103 Vec x) {
1104 SmartPetscObj<Mat> matrix;
1105 CHK_THROW_MESSAGE(
1106 createProjectionShellMatrix(ctx, x, ProjectionOperator::R, matrix),
1107 "Failed to create R projection matrix");
1108 return matrix;
1109}
1110
1111SmartPetscObj<Mat>
1112createProjectionMatrixRT(const boost::shared_ptr<ConstrainMatrixCtx> &ctx,
1113 Vec x) {
1114 SmartPetscObj<Mat> matrix;
1115 CHK_THROW_MESSAGE(
1116 createProjectionShellMatrix(ctx, x, ProjectionOperator::RT, matrix),
1117 "Failed to create R^T projection matrix");
1118 return matrix;
1119}
1120
1121SmartPetscObj<Mat> createProjectionMatrixCTCPlusQTKQ(
1122 const boost::shared_ptr<ConstrainMatrixCtx> &ctx, Vec x) {
1123 SmartPetscObj<Mat> matrix;
1124 CHK_THROW_MESSAGE(createProjectionShellMatrix(
1125 ctx, x, ProjectionOperator::CTCPlusQTKQ, matrix),
1126 "Failed to create C^T C + Q^T K Q projection matrix");
1127 return matrix;
1128}
1129
1130} // namespace MoFEM
active
bool active
Definition ConstrainMatrixCtx.cpp:65
projectionRT
PetscLogEvent projectionRT
Definition ConstrainMatrixCtx.cpp:35
projectionQ
PetscLogEvent projectionQ
Definition ConstrainMatrixCtx.cpp:32
event
PetscLogEvent event
Definition ConstrainMatrixCtx.cpp:63
C
SmartPetscObj< Mat > C
Definition ConstrainMatrixCtx.cpp:82
optionsPrefix
std::string optionsPrefix
Definition ConstrainMatrixCtx.cpp:86
projectionInit
PetscLogEvent projectionInit
Definition ConstrainMatrixCtx.cpp:31
nonzeroState
PetscObjectState nonzeroState
Definition ConstrainMatrixCtx.cpp:88
CT
SmartPetscObj< Mat > CT
Definition ConstrainMatrixCtx.cpp:83
CCT
SmartPetscObj< Mat > CCT
Definition ConstrainMatrixCtx.cpp:84
projectionP
PetscLogEvent projectionP
Definition ConstrainMatrixCtx.cpp:33
projectionR
PetscLogEvent projectionR
Definition ConstrainMatrixCtx.cpp:34
comm
MPI_Comm comm
Definition ConstrainMatrixCtx.cpp:64
projectionCTCPlusQTKQ
PetscLogEvent projectionCTCPlusQTKQ
Definition ConstrainMatrixCtx.cpp:36
configureFromOptions
bool configureFromOptions
Definition ConstrainMatrixCtx.cpp:89
matrixState
PetscObjectState matrixState
Definition ConstrainMatrixCtx.cpp:87
ksp
SmartPetscObj< KSP > ksp
Definition ConstrainMatrixCtx.cpp:85
initialized
bool initialized
Definition ConstrainMatrixCtx.cpp:90
inverse
SmartPetscObj< Mat > inverse
Definition TaoNullSpace.cpp:121
rows
SmartPetscObj< IS > rows
Definition TaoNullSpace.cpp:119
COL
@ COL
Definition definitions.h:136
ROW
@ ROW
Definition definitions.h:136
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition definitions.h:577
R
@ R
Definition free_surface.cpp:630
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:451
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:439
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::getCommFromPetscObject
MPI_Comm getCommFromPetscObject(PetscObject obj)
Get the Comm From Petsc Object object.
Definition PetscSmartObj.hpp:162
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:269
MoFEM::createCCTInverseKSP
SmartPetscObj< Mat > createCCTInverseKSP(Mat constrain_matrix, KSP ksp=nullptr, const std::string &options_prefix="constraint_matrix_")
Definition ConstrainMatrixCtx.cpp:1023
MoFEM::createProjectionMatrixP
SmartPetscObj< Mat > createProjectionMatrixP(const boost::shared_ptr< ConstrainMatrixCtx > &ctx, Vec x)
Definition ConstrainMatrixCtx.cpp:1092
MoFEM::createConstrainMatrixCtx
boost::shared_ptr< ConstrainMatrixCtx > createConstrainMatrixCtx(DM x_dm, DM y_dm, Mat constrain_matrix, Mat inverse_cct)
Definition ConstrainMatrixCtx.cpp:977
MoFEM::createProjectionMatrixR
SmartPetscObj< Mat > createProjectionMatrixR(const boost::shared_ptr< ConstrainMatrixCtx > &ctx, Vec x)
Definition ConstrainMatrixCtx.cpp:1102
MoFEM::createProjectionMatrixCTCPlusQTKQ
SmartPetscObj< Mat > createProjectionMatrixCTCPlusQTKQ(const boost::shared_ptr< ConstrainMatrixCtx > &ctx, Vec x)
Definition ConstrainMatrixCtx.cpp:1121
MoFEM::createProjectionMatrixQ
SmartPetscObj< Mat > createProjectionMatrixQ(const boost::shared_ptr< ConstrainMatrixCtx > &ctx, Vec x)
Definition ConstrainMatrixCtx.cpp:1082
MoFEM::getPetscObject
PetscObject getPetscObject(T obj)
Definition PetscSmartObj.hpp:10
MoFEM::createProjectionMatrixRT
SmartPetscObj< Mat > createProjectionMatrixRT(const boost::shared_ptr< ConstrainMatrixCtx > &ctx, Vec x)
Definition ConstrainMatrixCtx.cpp:1112
MoFEM::ConstrainMatrixCtxImpl::validateInverseCCTLayout
MoFEMErrorCode validateInverseCCTLayout(Mat constrain_matrix, Mat inverse_cct) const
Definition ConstrainMatrixCtx.cpp:236
MoFEM::ConstrainMatrixCtxImpl::applyQImpl
MoFEMErrorCode applyQImpl(Vec x, Vec f, bool transpose_inverse)
Definition ConstrainMatrixCtx.cpp:545
MoFEM::ConstrainMatrixCtxImpl::initializeQTKQ
MoFEMErrorCode initializeQTKQ() override
Definition ConstrainMatrixCtx.cpp:472
MoFEM::ConstrainMatrixCtxImpl::ConstrainMatrixCtxImpl
ConstrainMatrixCtxImpl(MPI_Comm comm, DM x_dm, DM y_dm)
Definition ConstrainMatrixCtx.cpp:180
MoFEM::ConstrainMatrixCtxImpl::getConstrainMatrix
Mat getConstrainMatrix() const override
Definition ConstrainMatrixCtx.cpp:355
MoFEM::ConstrainMatrixCtxImpl::applyR
MoFEMErrorCode applyR(Vec x, Vec f) override
Definition ConstrainMatrixCtx.cpp:599
MoFEM::ConstrainMatrixCtxImpl::applyCTC_QTKQ
MoFEMErrorCode applyCTC_QTKQ(Vec x, Vec f) override
Definition ConstrainMatrixCtx.cpp:658
MoFEM::ConstrainMatrixCtxImpl::validateVectorLayout
MoFEMErrorCode validateVectorLayout(Vec vec, PetscInt local_size, PetscInt global_size, const char *vector_name) const
Definition ConstrainMatrixCtx.cpp:215
MoFEM::ConstrainMatrixCtxImpl::validateSystemMatrixLayout
MoFEMErrorCode validateSystemMatrixLayout() const
Definition ConstrainMatrixCtx.cpp:290
MoFEM::ConstrainMatrixCtxImpl::setSystemMatrix
MoFEMErrorCode setSystemMatrix(Mat system_matrix) override
Definition ConstrainMatrixCtx.cpp:359
MoFEM::ConstrainMatrixCtxImpl::applyRT
MoFEMErrorCode applyRT(Vec x, Vec f) override
Definition ConstrainMatrixCtx.cpp:630
MoFEM::ConstrainMatrixCtxImpl::applyQ
MoFEMErrorCode applyQ(Vec x, Vec f) override
Definition ConstrainMatrixCtx.cpp:563
MoFEM::ConstrainMatrixCtxImpl::applyP
MoFEMErrorCode applyP(Vec x, Vec f) override
Definition ConstrainMatrixCtx.cpp:575
MoFEM::ConstrainMatrixCtxImpl::applyInverseCCT
MoFEMErrorCode applyInverseCCT(Vec x, Vec f, bool transpose)
Definition ConstrainMatrixCtx.cpp:525
MoFEM::ConstrainMatrixCtxImpl::setConstrainMatrix
MoFEMErrorCode setConstrainMatrix(Mat constrain_matrix, Mat inverse_cct) override
Definition ConstrainMatrixCtx.cpp:345
MoFEM::ConstrainMatrixCtxImpl::recalculateCT
MoFEMErrorCode recalculateCT() override
Definition ConstrainMatrixCtx.cpp:460
MoFEM::ConstrainMatrixCtxImpl::validatePetscObject
MoFEMErrorCode validatePetscObject(PetscObject object, const char *object_name) const
Definition ConstrainMatrixCtx.cpp:192
MoFEM::ConstrainMatrixCtxImpl::initializeQorP
MoFEMErrorCode initializeQorP(Vec x) override
Definition ConstrainMatrixCtx.cpp:369
MoFEM::ConstrainMatrixCtxImpl::~ConstrainMatrixCtxImpl
~ConstrainMatrixCtxImpl() override=default
MoFEM::ConstrainMatrixCtxImpl::recalculateCTC
MoFEMErrorCode recalculateCTC() override
Definition ConstrainMatrixCtx.cpp:515
MoFEM::ConstrainMatrixCtxImpl::scatter
SmartPetscObj< VecScatter > scatter
Definition ConstrainMatrixCtx.cpp:170
MoFEM::ConstrainMatrixCtx
API for matrices that project onto, or lift from, a constraint space.
Definition ConstrainMatrixCtx.hpp:16
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:103
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition VecManager.hpp:23
Generated by Doxygen 1.9.8 and hosted at University of Glasgow