v0.15.0
Loading...
Searching...
No Matches
Namespaces | Macros | Functions | Variables
ConstrainMatrixCtx.cpp File Reference

Implementation of projection matrix. More...

#include <ConstrainMatrixCtx.hpp>

Go to the source code of this file.

Namespaces

namespace  MoFEM
 implementation of Data Operators for Forces and Sources
 

Macros

#define INIT_DATA_CONSTRAINMATRIXCTX
 

Functions

MoFEMErrorCode MoFEM::ProjectionMatrixMultOpQ (Mat Q, Vec x, Vec f)
 Multiplication operator for Q = I-CTC(CCT)^-1C.
 
MoFEMErrorCode MoFEM::ConstrainMatrixMultOpP (Mat P, Vec x, Vec f)
 Multiplication operator for P = CT(CCT)^-1C.
 
MoFEMErrorCode MoFEM::ConstrainMatrixMultOpR (Mat R, Vec x, Vec f)
 Multiplication operator for R = CT(CCT)^-1.
 
MoFEMErrorCode MoFEM::ConstrainMatrixMultOpRT (Mat RT, Vec x, Vec f)
 Multiplication operator for RT = (CCT)^-TC.
 
MoFEMErrorCode MoFEM::ConstrainMatrixMultOpCTC_QTKQ (Mat CTC_QTKQ, Vec x, Vec f)
 Multiplication operator for RT = (CCT)^-TC.
 
MoFEMErrorCode MoFEM::ConstrainMatrixDestroyOpPorQ (Mat Q)
 Destroy shell matrix Q.
 
MoFEMErrorCode MoFEM::ConstrainMatrixDestroyOpQTKQ (Mat QTKQ)
 Destroy shell matrix.
 

Variables

static const bool MoFEM::debug = false
 

Detailed Description

Implementation of projection matrix.

Definition in file ConstrainMatrixCtx.cpp.

Macro Definition Documentation

◆ INIT_DATA_CONSTRAINMATRIXCTX

#define INIT_DATA_CONSTRAINMATRIXCTX
Value:
C(PETSC_NULLPTR), CT(PETSC_NULLPTR), CCT(PETSC_NULLPTR), CTC(PETSC_NULLPTR), \
K(PETSC_NULLPTR), Cx(PETSC_NULLPTR), CCTm1_Cx(PETSC_NULLPTR), \
CT_CCTm1_Cx(PETSC_NULLPTR), CTCx(PETSC_NULLPTR), Qx(PETSC_NULLPTR), \
KQx(PETSC_NULLPTR), initQorP(true), initQTKQ(true), createKSP(create_ksp), \
createScatter(true), cancelKSPMonitor(true), \
ownConstrainMatrix(own_contrain_matrix)

Definition at line 11 of file ConstrainMatrixCtx.cpp.

22 : mField(m_field), INIT_DATA_CONSTRAINMATRIXCTX, xProblem(x_problem),
23 yProblem(y_problem) {
24 PetscLogEventRegister("ProjectionInit", 0, &MOFEM_EVENT_projInit);
25 PetscLogEventRegister("ProjectionQ", 0, &MOFEM_EVENT_projQ);
26 PetscLogEventRegister("ProjectionP", 0, &MOFEM_EVENT_projP);
27 PetscLogEventRegister("ProjectionR", 0, &MOFEM_EVENT_projR);
28 PetscLogEventRegister("ProjectionRT", 0, &MOFEM_EVENT_projRT);
29 PetscLogEventRegister("ProjectionCTC_QTKQ", 0, &MOFEM_EVENT_projCTC_QTKQ);
30}
31
32ConstrainMatrixCtx::ConstrainMatrixCtx(MoFEM::Interface &m_field,
33 VecScatter scatter, bool create_ksp,
34 bool own_contrain_matrix)
35 : mField(m_field), INIT_DATA_CONSTRAINMATRIXCTX, sCatter(scatter) {
36 PetscLogEventRegister("ProjectionInit", 0, &MOFEM_EVENT_projInit);
37 PetscLogEventRegister("ProjectionQ", 0, &MOFEM_EVENT_projQ);
38 PetscLogEventRegister("ProjectionP", 0, &MOFEM_EVENT_projP);
39 PetscLogEventRegister("ProjectionR", 0, &MOFEM_EVENT_projR);
40 PetscLogEventRegister("ProjectionRT", 0, &MOFEM_EVENT_projRT);
41 PetscLogEventRegister("ProjectionCTC_QTKQ", 0, &MOFEM_EVENT_projCTC_QTKQ);
42}
43
44MoFEMErrorCode ConstrainMatrixCtx::initializeQorP(Vec x) {
45 MoFEMFunctionBegin;
46 if (initQorP) {
47 initQorP = false;
48
49 PetscLogEventBegin(MOFEM_EVENT_projInit, 0, 0, 0, 0);
50 CHKERR MatTranspose(C, MAT_INITIAL_MATRIX, &CT);
51 // need to be calculated when C is changed
52 CHKERR MatTransposeMatMult(CT, CT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCT);
53 if (createKSP) {
54 CHKERR KSPCreate(mField.get_comm(), &kSP);
55 // need to be recalculated when C is changed
56 CHKERR KSPSetOperators(kSP, CCT, CCT);
57 CHKERR KSPSetFromOptions(kSP);
58 CHKERR KSPSetInitialGuessKnoll(kSP, PETSC_TRUE);
59 CHKERR KSPGetTolerances(kSP, &rTol, &absTol, &dTol, &maxIts);
60 CHKERR KSPSetUp(kSP);
61 if (cancelKSPMonitor) {
62 CHKERR KSPMonitorCancel(kSP);
63 }
64 }
65#if PETSC_VERSION_GE(3, 5, 3)
66 CHKERR MatCreateVecs(C, &X, PETSC_NULLPTR);
67 CHKERR MatCreateVecs(C, PETSC_NULLPTR, &Cx);
68 CHKERR MatCreateVecs(CCT, PETSC_NULLPTR, &CCTm1_Cx);
69#else
70 CHKERR MatGetVecs(C, &X, PETSC_NULLPTR);
71 CHKERR MatGetVecs(C, PETSC_NULLPTR, &Cx);
72 CHKERR MatGetVecs(CCT, PETSC_NULLPTR, &CCTm1_Cx);
73#endif
74 CHKERR VecDuplicate(X, &CT_CCTm1_Cx);
75 if (createScatter) {
76 CHKERR mField.getInterface<VecManager>()->vecScatterCreate(
77 x, xProblem, ROW, X, yProblem, COL, &sCatter);
78 }
79 PetscLogEventEnd(MOFEM_EVENT_projInit, 0, 0, 0, 0);
80 }
81 MoFEMFunctionReturn(0);
82}
83
84MoFEMErrorCode ConstrainMatrixCtx::recalculateCTandCCT() {
85 MoFEMFunctionBegin;
86 if (initQorP)
87 MoFEMFunctionReturnHot(0);
88 CHKERR MatTranspose(C, MAT_REUSE_MATRIX, &CT);
89 CHKERR MatTransposeMatMult(CT, CT, MAT_REUSE_MATRIX, PETSC_DEFAULT, &CCT);
90 MoFEMFunctionReturn(0);
91}
92
93MoFEMErrorCode ConstrainMatrixCtx::destroyQorP() {
94 MoFEMFunctionBegin;
95 if (initQorP)
96 MoFEMFunctionReturnHot(0);
97 CHKERR MatDestroy(&CT);
98 CHKERR MatDestroy(&CCT);
99 if (createKSP) {
100 CHKERR KSPDestroy(&kSP);
101 }
102 CHKERR VecDestroy(&X);
103 CHKERR VecDestroy(&Cx);
104 CHKERR VecDestroy(&CCTm1_Cx);
105 CHKERR VecDestroy(&CT_CCTm1_Cx);
106 if (createScatter) {
107 CHKERR VecScatterDestroy(&sCatter);
108 }
109 initQorP = true;
110 MoFEMFunctionReturn(0);
111}
112
113MoFEMErrorCode ConstrainMatrixCtx::initializeQTKQ() {
114 MoFEMFunctionBegin;
115 if (initQTKQ) {
116 initQTKQ = false;
117 PetscLogEventBegin(MOFEM_EVENT_projInit, 0, 0, 0, 0);
118 // need to be recalculated when C is changed
119 CHKERR MatTransposeMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CTC);
120 if (debug) {
121 // MatView(CCT,PETSC_VIEWER_DRAW_WORLD);
122 int m, n;
123 MatGetSize(CCT, &m, &n);
124 PetscPrintf(mField.get_comm(), "CCT size (%d,%d)\n", m, n);
125 // std::string wait;
126 // std::cin >> wait;
127 }
128#if PETSC_VERSION_GE(3, 5, 3)
129 CHKERR MatCreateVecs(K, &Qx, PETSC_NULLPTR);
130 CHKERR MatCreateVecs(K, PETSC_NULLPTR, &KQx);
131 CHKERR MatCreateVecs(CTC, PETSC_NULLPTR, &CTCx);
132#else
133 CHKERR MatGetVecs(K, &Qx, PETSC_NULLPTR);
134 CHKERR MatGetVecs(K, PETSC_NULLPTR, &KQx);
135 CHKERR MatGetVecs(CTC, PETSC_NULLPTR, &CTCx);
136#endif
137 PetscLogEventEnd(MOFEM_EVENT_projInit, 0, 0, 0, 0);
138 }
139 MoFEMFunctionReturn(0);
140}
141
142MoFEMErrorCode ConstrainMatrixCtx::recalculateCTC() {
143 MoFEMFunctionBegin;
144 if (initQTKQ)
145 MoFEMFunctionReturnHot(0);
146 CHKERR MatTransposeMatMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &CTC);
147 MoFEMFunctionReturn(0);
148}
149
150MoFEMErrorCode ConstrainMatrixCtx::destroyQTKQ() {
151 MoFEMFunctionBegin;
152 if (initQTKQ)
153 MoFEMFunctionReturnHot(0);
154 CHKERR MatDestroy(&CTC);
155 CHKERR VecDestroy(&Qx);
156 CHKERR VecDestroy(&KQx);
157 CHKERR VecDestroy(&CTCx);
158 initQTKQ = true;
159 MoFEMFunctionReturn(0);
160}
161
162MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f) {
163 MoFEMFunctionBegin;
164 void *void_ctx;
165 CHKERR MatShellGetContext(Q, &void_ctx);
166 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
167 PetscLogEventBegin(ctx->MOFEM_EVENT_projQ, 0, 0, 0, 0);
168 CHKERR ctx->initializeQorP(x);
169 CHKERR VecCopy(x, f);
170 CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
171 SCATTER_FORWARD);
172 CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
173 if (debug) {
174 // CHKERR VecView(ctx->X,PETSC_VIEWER_STDOUT_WORLD);
175 CHKERR VecScatterBegin(ctx->sCatter, ctx->X, f, INSERT_VALUES,
176 SCATTER_REVERSE);
177 CHKERR VecScatterEnd(ctx->sCatter, ctx->X, f, INSERT_VALUES,
178 SCATTER_REVERSE);
179 PetscBool flg;
180 CHKERR VecEqual(x, f, &flg);
181 if (flg == PETSC_FALSE)
182 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "scatter is not working");
183 }
184 CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
185 CHKERR KSPSolve(ctx->kSP, ctx->Cx, ctx->CCTm1_Cx);
186 CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
187 CHKERR VecScale(ctx->CT_CCTm1_Cx, -1);
188 CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, ADD_VALUES,
189 SCATTER_REVERSE);
190 CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, ADD_VALUES,
191 SCATTER_REVERSE);
192 PetscLogEventEnd(ctx->MOFEM_EVENT_projQ, 0, 0, 0, 0);
193 MoFEMFunctionReturn(0);
194}
195
196MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f) {
197 MoFEMFunctionBegin;
198 void *void_ctx;
199 CHKERR MatShellGetContext(P, &void_ctx);
200 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
201 PetscLogEventBegin(ctx->MOFEM_EVENT_projP, 0, 0, 0, 0);
202 CHKERR ctx->initializeQorP(x);
203 CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
204 SCATTER_FORWARD);
205 CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
206 CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
207 CHKERR KSPSolve(ctx->kSP, ctx->Cx, ctx->CCTm1_Cx);
208 CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
209 CHKERR VecZeroEntries(f);
210 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
211 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
212 CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
213 SCATTER_REVERSE);
214 CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
215 SCATTER_REVERSE);
216 PetscLogEventEnd(ctx->MOFEM_EVENT_projP, 0, 0, 0, 0);
217 MoFEMFunctionReturn(0);
218}
219
220MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f) {
221 MoFEMFunctionBegin;
222 void *void_ctx;
223 CHKERR MatShellGetContext(R, &void_ctx);
224 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
225 PetscLogEventBegin(ctx->MOFEM_EVENT_projR, 0, 0, 0, 0);
226 if (ctx->initQorP)
227 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
228 "you have to call first initQorP or use Q matrix");
229 CHKERR KSPSolve(ctx->kSP, x, ctx->CCTm1_Cx);
230 CHKERR MatMult(ctx->CT, ctx->CCTm1_Cx, ctx->CT_CCTm1_Cx);
231 CHKERR VecZeroEntries(f);
232 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
233 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
234 CHKERR VecScatterBegin(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
235 SCATTER_REVERSE);
236 CHKERR VecScatterEnd(ctx->sCatter, ctx->CT_CCTm1_Cx, f, INSERT_VALUES,
237 SCATTER_REVERSE);
238 PetscLogEventEnd(ctx->MOFEM_EVENT_projR, 0, 0, 0, 0);
239 MoFEMFunctionReturn(0);
240}
241
242MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f) {
243 MoFEMFunctionBegin;
244 void *void_ctx;
245 CHKERR MatShellGetContext(RT, &void_ctx);
246 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
247 PetscLogEventBegin(ctx->MOFEM_EVENT_projRT, 0, 0, 0, 0);
248 if (ctx->initQorP)
249 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
250 "you have to call first initQorP or use Q matrix");
251 CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
252 SCATTER_FORWARD);
253 CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
254 CHKERR MatMult(ctx->C, ctx->X, ctx->Cx);
255 CHKERR KSPSolve(ctx->kSP, ctx->Cx, f);
256 PetscLogEventEnd(ctx->MOFEM_EVENT_projRT, 0, 0, 0, 0);
257 MoFEMFunctionReturn(0);
258}
259
260MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x, Vec f) {
261 MoFEMFunctionBegin;
262 void *void_ctx;
263 CHKERR MatShellGetContext(CTC_QTKQ, &void_ctx);
264 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
265 PetscLogEventBegin(ctx->MOFEM_EVENT_projCTC_QTKQ, 0, 0, 0, 0);
266 Mat Q;
267 int M, N, m, n;
268 CHKERR MatGetSize(ctx->K, &M, &N);
269 CHKERR MatGetLocalSize(ctx->K, &m, &n);
270 CHKERR MatCreateShell(ctx->mField.get_comm(), m, n, M, N, ctx, &Q);
271 CHKERR MatShellSetOperation(Q, MATOP_MULT,
272 (void (*)(void))ProjectionMatrixMultOpQ);
273 CHKERR ctx->initializeQTKQ();
274 CHKERR MatMult(Q, x, ctx->Qx);
275 CHKERR MatMult(ctx->K, ctx->Qx, ctx->KQx);
276 CHKERR MatMult(Q, ctx->KQx, f);
277 CHKERR VecScatterBegin(ctx->sCatter, x, ctx->X, INSERT_VALUES,
278 SCATTER_FORWARD);
279 CHKERR VecScatterEnd(ctx->sCatter, x, ctx->X, INSERT_VALUES, SCATTER_FORWARD);
280 CHKERR MatMult(ctx->CTC, ctx->X, ctx->CTCx);
281 CHKERR VecScatterBegin(ctx->sCatter, ctx->CTCx, f, ADD_VALUES,
282 SCATTER_REVERSE);
283 CHKERR VecScatterEnd(ctx->sCatter, ctx->CTCx, f, ADD_VALUES, SCATTER_REVERSE);
284 CHKERR MatDestroy(&Q);
285 PetscLogEventEnd(ctx->MOFEM_EVENT_projCTC_QTKQ, 0, 0, 0, 0);
286 MoFEMFunctionReturn(0);
287}
288
289MoFEMErrorCode ConstrainMatrixDestroyOpPorQ(Mat Q) {
290 MoFEMFunctionBegin;
291 void *void_ctx;
292 CHKERR MatShellGetContext(Q, &void_ctx);
293 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
294 CHKERR ctx->destroyQorP();
295 MoFEMFunctionReturn(0);
296}
297MoFEMErrorCode ConstrainMatrixDestroyOpQTKQ(Mat QTKQ) {
298 MoFEMFunctionBegin;
299 void *void_ctx;
300 CHKERR MatShellGetContext(QTKQ, &void_ctx);
301 ConstrainMatrixCtx *ctx = (ConstrainMatrixCtx *)void_ctx;
302 CHKERR ctx->destroyQTKQ();
303 MoFEMFunctionReturn(0);
304}
305
306}
INIT_DATA_CONSTRAINMATRIXCTX
#define INIT_DATA_CONSTRAINMATRIXCTX
Definition ConstrainMatrixCtx.cpp:11
COL
@ COL
Definition definitions.h:136
ROW
@ ROW
Definition definitions.h:136
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_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
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
debug
static const bool debug
Definition dm_create_subdm.cpp:12
R
@ R
Definition free_surface.cpp:624
MoFEM::ConstrainMatrixMultOpCTC_QTKQ
MoFEMErrorCode ConstrainMatrixMultOpCTC_QTKQ(Mat CTC_QTKQ, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
Definition ConstrainMatrixCtx.cpp:261
MoFEM::ConstrainMatrixMultOpP
MoFEMErrorCode ConstrainMatrixMultOpP(Mat P, Vec x, Vec f)
Multiplication operator for P = CT(CCT)^-1C.
Definition ConstrainMatrixCtx.cpp:197
MoFEM::ConstrainMatrixDestroyOpPorQ
MoFEMErrorCode ConstrainMatrixDestroyOpPorQ(Mat Q)
Destroy shell matrix Q.
Definition ConstrainMatrixCtx.cpp:290
MoFEM::ProjectionMatrixMultOpQ
MoFEMErrorCode ProjectionMatrixMultOpQ(Mat Q, Vec x, Vec f)
Multiplication operator for Q = I-CTC(CCT)^-1C.
Definition ConstrainMatrixCtx.cpp:163
MoFEM::ConstrainMatrixMultOpR
MoFEMErrorCode ConstrainMatrixMultOpR(Mat R, Vec x, Vec f)
Multiplication operator for R = CT(CCT)^-1.
Definition ConstrainMatrixCtx.cpp:221
MoFEM::ConstrainMatrixDestroyOpQTKQ
MoFEMErrorCode ConstrainMatrixDestroyOpQTKQ(Mat QTKQ)
Destroy shell matrix.
Definition ConstrainMatrixCtx.cpp:298
MoFEM::ConstrainMatrixMultOpRT
MoFEMErrorCode ConstrainMatrixMultOpRT(Mat RT, Vec x, Vec f)
Multiplication operator for RT = (CCT)^-TC.
Definition ConstrainMatrixCtx.cpp:243
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition PlasticOps.hpp:117
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
N
const int N
Definition speed_test.cpp:3
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
Generated by Doxygen 1.9.8 and hosted at University of Glasgow