1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/
2 #include <petsc/private/pcbddcimpl.h>
3 #include <petsc/private/pcbddcprivateimpl.h>
4 #include <petscdm.h>
5
6 static PetscBool cited = PETSC_FALSE;
7 static PetscBool cited2 = PETSC_FALSE;
8 static const char citation[] = "@article{ZampiniPCBDDC,\n"
9 "author = {Stefano Zampini},\n"
10 "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
11 "journal = {SIAM Journal on Scientific Computing},\n"
12 "volume = {38},\n"
13 "number = {5},\n"
14 "pages = {S282-S306},\n"
15 "year = {2016},\n"
16 "doi = {10.1137/15M1025785},\n"
17 "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
18 "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
19 "}\n"
20 "@article{ZampiniDualPrimal,\n"
21 "author = {Stefano Zampini},\n"
22 "title = {{D}ual-{P}rimal methods for the cardiac {B}idomain model},\n"
23 "volume = {24},\n"
24 "number = {04},\n"
25 "pages = {667-696},\n"
26 "year = {2014},\n"
27 "doi = {10.1142/S0218202513500632},\n"
28 "URL = {https://www.worldscientific.com/doi/abs/10.1142/S0218202513500632},\n"
29 "eprint = {https://www.worldscientific.com/doi/pdf/10.1142/S0218202513500632}\n"
30 "}\n";
31 static const char citation2[] = "@article{li2013nonoverlapping,\n"
32 "title={A nonoverlapping domain decomposition method for incompressible Stokes equations with continuous pressures},\n"
33 "author={Li, Jing and Tu, Xuemin},\n"
34 "journal={SIAM Journal on Numerical Analysis},\n"
35 "volume={51},\n"
36 "number={2},\n"
37 "pages={1235--1253},\n"
38 "year={2013},\n"
39 "publisher={Society for Industrial and Applied Mathematics}\n"
40 "}\n";
41
42 /*
43 This file implements the FETI-DP method in PETSc as part of KSP.
44 */
45 typedef struct {
46 KSP parentksp;
47 } KSP_FETIDPMon;
48
49 typedef struct {
50 KSP innerksp; /* the KSP for the Lagrange multipliers */
51 PC innerbddc; /* the inner BDDC object */
52 PetscBool fully_redundant; /* true for using a fully redundant set of multipliers */
53 PetscBool userbddc; /* true if the user provided the PCBDDC object */
54 PetscBool saddlepoint; /* support for saddle point problems */
55 IS pP; /* index set for pressure variables */
56 Vec rhs_flip; /* see KSPFETIDPSetUpOperators */
57 KSP_FETIDPMon *monctx; /* monitor context, used to pass user defined monitors
58 in the physical space */
59 PetscObjectState matstate; /* these are needed just in the saddle point case */
60 PetscObjectState matnnzstate; /* where we are going to use MatZeroRows on pmat */
61 PetscBool statechanged;
62 PetscBool check;
63 } KSP_FETIDP;
64
KSPFETIDPSetPressureOperator_FETIDP(KSP ksp,Mat P)65 static PetscErrorCode KSPFETIDPSetPressureOperator_FETIDP(KSP ksp, Mat P)
66 {
67 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
68
69 PetscFunctionBegin;
70 if (P) fetidp->saddlepoint = PETSC_TRUE;
71 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject)P));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
75 /*@
76 KSPFETIDPSetPressureOperator - Sets the operator used to set up the pressure preconditioner for the saddle point `KSPFETIDP` solver,
77
78 Collective
79
80 Input Parameters:
81 + ksp - the `KSPFETIDP` solver
82 - P - the linear operator to be preconditioned, usually the mass matrix.
83
84 Level: advanced
85
86 Notes:
87 The operator can be either passed in
88 .vb
89 a) monolithic global ordering,
90 b) pressure-only global ordering, or
91 c) interface pressure ordering (if `-ksp_fetidp_pressure_all false`).
92 .ve
93 In cases b) and c), the pressure ordering of dofs needs to satisfy
94 pid_1 < pid_2 iff gid_1 < gid_2
95 where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding
96 id in the monolithic global ordering.
97
98 .seealso: [](ch_ksp), `KSPFETIDP`, `MATIS`, `PCBDDC`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`, `KSPSetOperators()`
99 @*/
KSPFETIDPSetPressureOperator(KSP ksp,Mat P)100 PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P)
101 {
102 PetscFunctionBegin;
103 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
104 if (P) PetscValidHeaderSpecific(P, MAT_CLASSID, 2);
105 PetscTryMethod(ksp, "KSPFETIDPSetPressureOperator_C", (KSP, Mat), (ksp, P));
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
KSPFETIDPGetInnerKSP_FETIDP(KSP ksp,KSP * innerksp)109 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP *innerksp)
110 {
111 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
112
113 PetscFunctionBegin;
114 *innerksp = fetidp->innerksp;
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /*@
119 KSPFETIDPGetInnerKSP - Gets the `KSP` object for the Lagrange multipliers from inside a `KSPFETIDP`
120
121 Input Parameter:
122 . ksp - the `KSPFETIDP`
123
124 Output Parameter:
125 . innerksp - the `KSP` for the multipliers
126
127 Level: advanced
128
129 .seealso: [](ch_ksp), `KSPFETIDP`, `MATIS`, `PCBDDC`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerBDDC()`
130 @*/
KSPFETIDPGetInnerKSP(KSP ksp,KSP * innerksp)131 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP *innerksp)
132 {
133 PetscFunctionBegin;
134 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
135 PetscAssertPointer(innerksp, 2);
136 PetscUseMethod(ksp, "KSPFETIDPGetInnerKSP_C", (KSP, KSP *), (ksp, innerksp));
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp,PC * pc)140 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC *pc)
141 {
142 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
143
144 PetscFunctionBegin;
145 *pc = fetidp->innerbddc;
146 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
149 /*@
150 KSPFETIDPGetInnerBDDC - Gets the `PCBDDC` preconditioner used to set up the `KSPFETIDP` matrix for the Lagrange multipliers
151
152 Input Parameter:
153 . ksp - the `KSPFETIDP` Krylov solver
154
155 Output Parameter:
156 . pc - the `PCBDDC` preconditioner
157
158 Level: advanced
159
160 .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDP`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
161 @*/
KSPFETIDPGetInnerBDDC(KSP ksp,PC * pc)162 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC *pc)
163 {
164 PetscFunctionBegin;
165 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
166 PetscAssertPointer(pc, 2);
167 PetscUseMethod(ksp, "KSPFETIDPGetInnerBDDC_C", (KSP, PC *), (ksp, pc));
168 PetscFunctionReturn(PETSC_SUCCESS);
169 }
170
KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp,PC pc)171 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
172 {
173 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
174
175 PetscFunctionBegin;
176 PetscCall(PetscObjectReference((PetscObject)pc));
177 PetscCall(PCDestroy(&fetidp->innerbddc));
178 fetidp->innerbddc = pc;
179 fetidp->userbddc = PETSC_TRUE;
180 PetscFunctionReturn(PETSC_SUCCESS);
181 }
182
183 /*@
184 KSPFETIDPSetInnerBDDC - Provides the `PCBDDC` preconditioner used to set up the `KSPFETIDP` matrix for the Lagrange multipliers
185
186 Collective
187
188 Input Parameters:
189 + ksp - the `KSPFETIDP` Krylov solver
190 - pc - the `PCBDDC` preconditioner
191
192 Level: advanced
193
194 Note:
195 A `PC` is automatically created for the `KSPFETIDP` and can be accessed to change options with `KSPFETIDPGetInnerBDDC()` hence this routine is rarely needed
196
197 .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
198 @*/
KSPFETIDPSetInnerBDDC(KSP ksp,PC pc)199 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
200 {
201 PetscBool isbddc;
202
203 PetscFunctionBegin;
204 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
205 PetscValidHeaderSpecific(pc, PC_CLASSID, 2);
206 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
207 PetscCheck(isbddc, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_WRONG, "KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner");
208 PetscTryMethod(ksp, "KSPFETIDPSetInnerBDDC_C", (KSP, PC), (ksp, pc));
209 PetscFunctionReturn(PETSC_SUCCESS);
210 }
211
KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec * V)212 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp, Vec v, Vec *V)
213 {
214 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
215 Mat F;
216 Vec Xl;
217
218 PetscFunctionBegin;
219 PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
220 PetscCall(KSPBuildSolution(fetidp->innerksp, NULL, &Xl));
221 if (v) {
222 PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, v));
223 *V = v;
224 } else {
225 PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, *V));
226 }
227 PetscFunctionReturn(PETSC_SUCCESS);
228 }
229
KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,PetscCtx ctx)230 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp, PetscInt it, PetscReal rnorm, PetscCtx ctx)
231 {
232 KSP_FETIDPMon *monctx = (KSP_FETIDPMon *)ctx;
233
234 PetscFunctionBegin;
235 PetscCall(KSPMonitor(monctx->parentksp, it, rnorm));
236 PetscFunctionReturn(PETSC_SUCCESS);
237 }
238
KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal * r,PetscReal * c,PetscInt * neig)239 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
240 {
241 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
242
243 PetscFunctionBegin;
244 PetscCall(KSPComputeEigenvalues(fetidp->innerksp, nmax, r, c, neig));
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal * emax,PetscReal * emin)248 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp, PetscReal *emax, PetscReal *emin)
249 {
250 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
251
252 PetscFunctionBegin;
253 PetscCall(KSPComputeExtremeSingularValues(fetidp->innerksp, emax, emin));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
KSPFETIDPCheckOperators(KSP ksp,PetscViewer viewer)257 static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer)
258 {
259 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
260 PC_BDDC *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
261 PC_IS *pcis = (PC_IS *)fetidp->innerbddc->data;
262 Mat_IS *matis = (Mat_IS *)fetidp->innerbddc->pmat->data;
263 Mat F;
264 FETIDPMat_ctx fetidpmat_ctx;
265 Vec test_vec, test_vec_p = NULL, fetidp_global;
266 IS dirdofs, isvert;
267 MPI_Comm comm = PetscObjectComm((PetscObject)ksp);
268 PetscScalar sval, *array;
269 PetscReal val, rval;
270 const PetscInt *vertex_indices;
271 PetscInt i, n_vertices;
272 PetscBool isascii;
273
274 PetscFunctionBegin;
275 PetscCheckSameComm(ksp, 1, viewer, 2);
276 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
277 PetscCheck(isascii, comm, PETSC_ERR_SUP, "Unsupported viewer");
278 PetscCall(PetscViewerASCIIPrintf(viewer, "----------FETI-DP MAT --------------\n"));
279 PetscCall(PetscViewerASCIIAddTab(viewer, 2));
280 PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
281 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
282 PetscCall(MatView(F, viewer));
283 PetscCall(PetscViewerPopFormat(viewer));
284 PetscCall(PetscViewerASCIISubtractTab(viewer, 2));
285 PetscCall(MatShellGetContext(F, &fetidpmat_ctx));
286 PetscCall(PetscViewerASCIIPrintf(viewer, "----------FETI-DP TESTS--------------\n"));
287 PetscCall(PetscViewerASCIIPrintf(viewer, "All tests should return zero!\n"));
288 PetscCall(PetscViewerASCIIPrintf(viewer, "FETIDP MAT context in the "));
289 if (fetidp->fully_redundant) {
290 PetscCall(PetscViewerASCIIPrintf(viewer, "fully redundant case for lagrange multipliers.\n"));
291 } else {
292 PetscCall(PetscViewerASCIIPrintf(viewer, "Non-fully redundant case for lagrange multiplier.\n"));
293 }
294 PetscCall(PetscViewerFlush(viewer));
295
296 /* Get Vertices used to define the BDDC */
297 PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
298 PetscCall(ISGetLocalSize(isvert, &n_vertices));
299 PetscCall(ISGetIndices(isvert, &vertex_indices));
300
301 /******************************************************************/
302 /* TEST A/B: Test numbering of global fetidp dofs */
303 /******************************************************************/
304 PetscCall(MatCreateVecs(F, &fetidp_global, NULL));
305 PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &test_vec));
306 PetscCall(VecSet(fetidp_global, 1.0));
307 PetscCall(VecSet(test_vec, 1.));
308 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
309 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
310 if (fetidpmat_ctx->l2g_p) {
311 PetscCall(VecDuplicate(fetidpmat_ctx->vP, &test_vec_p));
312 PetscCall(VecSet(test_vec_p, 1.));
313 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidp_global, fetidpmat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
314 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidp_global, fetidpmat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
315 }
316 PetscCall(VecAXPY(test_vec, -1.0, fetidpmat_ctx->lambda_local));
317 PetscCall(VecNorm(test_vec, NORM_INFINITY, &val));
318 PetscCall(VecDestroy(&test_vec));
319 PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
320 PetscCall(PetscViewerASCIIPrintf(viewer, "A: CHECK glob to loc: % 1.14e\n", (double)rval));
321
322 if (fetidpmat_ctx->l2g_p) {
323 PetscCall(VecAXPY(test_vec_p, -1.0, fetidpmat_ctx->vP));
324 PetscCall(VecNorm(test_vec_p, NORM_INFINITY, &val));
325 PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
326 PetscCall(PetscViewerASCIIPrintf(viewer, "A: CHECK glob to loc (p): % 1.14e\n", (double)rval));
327 }
328
329 if (fetidp->fully_redundant) {
330 PetscCall(VecSet(fetidp_global, 0.0));
331 PetscCall(VecSet(fetidpmat_ctx->lambda_local, 0.5));
332 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
333 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
334 PetscCall(VecSum(fetidp_global, &sval));
335 val = PetscRealPart(sval) - fetidpmat_ctx->n_lambda;
336 PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
337 PetscCall(PetscViewerASCIIPrintf(viewer, "B: CHECK loc to glob: % 1.14e\n", (double)rval));
338 }
339
340 if (fetidpmat_ctx->l2g_p) {
341 PetscCall(VecSet(pcis->vec1_N, 1.0));
342 PetscCall(VecSet(pcis->vec1_global, 0.0));
343 PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
344 PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
345
346 PetscCall(VecSet(fetidp_global, 0.0));
347 PetscCall(VecSet(fetidpmat_ctx->vP, -1.0));
348 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
349 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
350 PetscCall(VecScatterBegin(fetidpmat_ctx->g2g_p, fetidp_global, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
351 PetscCall(VecScatterEnd(fetidpmat_ctx->g2g_p, fetidp_global, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
352 PetscCall(VecScatterBegin(fetidpmat_ctx->g2g_p, pcis->vec1_global, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
353 PetscCall(VecScatterEnd(fetidpmat_ctx->g2g_p, pcis->vec1_global, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
354 PetscCall(VecSum(fetidp_global, &sval));
355 val = PetscRealPart(sval);
356 PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
357 PetscCall(PetscViewerASCIIPrintf(viewer, "B: CHECK loc to glob (p): % 1.14e\n", (double)rval));
358 }
359
360 /******************************************************************/
361 /* TEST C: It should hold B_delta*w=0, w\in\widehat{W} */
362 /* This is the meaning of the B matrix */
363 /******************************************************************/
364
365 PetscCall(VecSetRandom(pcis->vec1_N, NULL));
366 PetscCall(VecSet(pcis->vec1_global, 0.0));
367 PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
368 PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
369 PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
370 PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
371 PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
372 PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
373 /* Action of B_delta */
374 PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
375 PetscCall(VecSet(fetidp_global, 0.0));
376 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
377 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
378 PetscCall(VecNorm(fetidp_global, NORM_INFINITY, &val));
379 PetscCall(PetscViewerASCIIPrintf(viewer, "C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n", (double)val));
380
381 /******************************************************************/
382 /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W} */
383 /* E_D = R_D^TR */
384 /* P_D = B_{D,delta}^T B_{delta} */
385 /* eq.44 Mandel Tezaur and Dohrmann 2005 */
386 /******************************************************************/
387
388 /* compute a random vector in \widetilde{W} */
389 PetscCall(VecSetRandom(pcis->vec1_N, NULL));
390 /* set zero at vertices and essential dofs */
391 PetscCall(VecGetArray(pcis->vec1_N, &array));
392 for (i = 0; i < n_vertices; i++) array[vertex_indices[i]] = 0.0;
393 PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirdofs));
394 if (dirdofs) {
395 const PetscInt *idxs;
396 PetscInt ndir;
397
398 PetscCall(ISGetLocalSize(dirdofs, &ndir));
399 PetscCall(ISGetIndices(dirdofs, &idxs));
400 for (i = 0; i < ndir; i++) array[idxs[i]] = 0.0;
401 PetscCall(ISRestoreIndices(dirdofs, &idxs));
402 }
403 PetscCall(VecRestoreArray(pcis->vec1_N, &array));
404 /* store w for final comparison */
405 PetscCall(VecDuplicate(pcis->vec1_B, &test_vec));
406 PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, test_vec, INSERT_VALUES, SCATTER_FORWARD));
407 PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, test_vec, INSERT_VALUES, SCATTER_FORWARD));
408
409 /* Jump operator P_D : results stored in pcis->vec1_B */
410 /* Action of B_delta */
411 PetscCall(MatMult(fetidpmat_ctx->B_delta, test_vec, fetidpmat_ctx->lambda_local));
412 PetscCall(VecSet(fetidp_global, 0.0));
413 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
414 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
415 /* Action of B_Ddelta^T */
416 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
417 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
418 PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));
419
420 /* Average operator E_D : results stored in pcis->vec2_B */
421 PetscCall(PCBDDCScalingExtension(fetidpmat_ctx->pc, test_vec, pcis->vec1_global));
422 PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec2_B, INSERT_VALUES, SCATTER_FORWARD));
423 PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec2_B, INSERT_VALUES, SCATTER_FORWARD));
424
425 /* test E_D=I-P_D */
426 PetscCall(VecAXPY(pcis->vec1_B, 1.0, pcis->vec2_B));
427 PetscCall(VecAXPY(pcis->vec1_B, -1.0, test_vec));
428 PetscCall(VecNorm(pcis->vec1_B, NORM_INFINITY, &val));
429 PetscCall(VecDestroy(&test_vec));
430 PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
431 PetscCall(PetscViewerASCIIPrintf(viewer, "%d: CHECK infty norm of E_D + P_D - I: %1.14e\n", PetscGlobalRank, (double)val));
432
433 /******************************************************************/
434 /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W} */
435 /* eq.48 Mandel Tezaur and Dohrmann 2005 */
436 /******************************************************************/
437
438 PetscCall(VecSetRandom(pcis->vec1_N, NULL));
439 /* set zero at vertices and essential dofs */
440 PetscCall(VecGetArray(pcis->vec1_N, &array));
441 for (i = 0; i < n_vertices; i++) array[vertex_indices[i]] = 0.0;
442 if (dirdofs) {
443 const PetscInt *idxs;
444 PetscInt ndir;
445
446 PetscCall(ISGetLocalSize(dirdofs, &ndir));
447 PetscCall(ISGetIndices(dirdofs, &idxs));
448 for (i = 0; i < ndir; i++) array[idxs[i]] = 0.0;
449 PetscCall(ISRestoreIndices(dirdofs, &idxs));
450 }
451 PetscCall(VecRestoreArray(pcis->vec1_N, &array));
452
453 /* Jump operator P_D : results stored in pcis->vec1_B */
454
455 PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
456 PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
457 /* Action of B_delta */
458 PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
459 PetscCall(VecSet(fetidp_global, 0.0));
460 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
461 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
462 /* Action of B_Ddelta^T */
463 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
464 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
465 PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));
466 /* scaling */
467 PetscCall(PCBDDCScalingExtension(fetidpmat_ctx->pc, pcis->vec1_B, pcis->vec1_global));
468 PetscCall(VecNorm(pcis->vec1_global, NORM_INFINITY, &val));
469 PetscCall(PetscViewerASCIIPrintf(viewer, "E: CHECK infty norm of R^T_D P_D: % 1.14e\n", (double)val));
470
471 if (!fetidp->fully_redundant) {
472 /******************************************************************/
473 /* TEST F: It should holds B_{delta}B^T_{D,delta}=I */
474 /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005 */
475 /******************************************************************/
476 PetscCall(VecDuplicate(fetidp_global, &test_vec));
477 PetscCall(VecSetRandom(fetidp_global, NULL));
478 if (fetidpmat_ctx->l2g_p) {
479 PetscCall(VecSet(fetidpmat_ctx->vP, 0.));
480 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
481 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
482 }
483 /* Action of B_Ddelta^T */
484 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
485 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
486 PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));
487 /* Action of B_delta */
488 PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
489 PetscCall(VecSet(test_vec, 0.0));
490 PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, test_vec, ADD_VALUES, SCATTER_FORWARD));
491 PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, test_vec, ADD_VALUES, SCATTER_FORWARD));
492 PetscCall(VecAXPY(fetidp_global, -1., test_vec));
493 PetscCall(VecNorm(fetidp_global, NORM_INFINITY, &val));
494 PetscCall(PetscViewerASCIIPrintf(viewer, "E: CHECK infty norm of P^T_D - I: % 1.14e\n", (double)val));
495 PetscCall(VecDestroy(&test_vec));
496 }
497 PetscCall(PetscViewerASCIIPrintf(viewer, "-------------------------------------\n"));
498 PetscCall(PetscViewerFlush(viewer));
499 PetscCall(VecDestroy(&test_vec_p));
500 PetscCall(ISDestroy(&dirdofs));
501 PetscCall(VecDestroy(&fetidp_global));
502 PetscCall(ISRestoreIndices(isvert, &vertex_indices));
503 PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
504 PetscFunctionReturn(PETSC_SUCCESS);
505 }
506
KSPFETIDPSetUpOperators(KSP ksp)507 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp)
508 {
509 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
510 PC_BDDC *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
511 Mat A, Ap;
512 PetscInt fidp[8] = {-1}, nfp = 8;
513 PetscMPIInt size;
514 PetscBool ismatis, pisz = PETSC_FALSE, allp = PETSC_FALSE, schp = PETSC_FALSE;
515 PetscBool flip = PETSC_FALSE; /* Usually, Stokes is written (B = -\int_\Omega \nabla \cdot u q)
516 | A B'| | v | = | f |
517 | B 0 | | p | = | g |
518 If -ksp_fetidp_saddlepoint_flip is true, the code assumes it is written as
519 | A B'| | v | = | f |
520 |-B 0 | | p | = |-g |
521 */
522 PetscObjectState matstate, matnnzstate;
523
524 PetscFunctionBegin;
525 PetscOptionsBegin(PetscObjectComm((PetscObject)ksp), ((PetscObject)ksp)->prefix, "FETI-DP options", "PC");
526 PetscCall(PetscOptionsIntArray("-ksp_fetidp_pressure_field", "Field id for pressures for saddle-point problems", NULL, fidp, &nfp, NULL));
527 PetscCall(PetscOptionsBool("-ksp_fetidp_pressure_all", "Use the whole pressure set instead of just that at the interface", NULL, allp, &allp, NULL));
528 PetscCall(PetscOptionsBool("-ksp_fetidp_saddlepoint_flip", "Flip the sign of the pressure-velocity (lower-left) block", NULL, flip, &flip, NULL));
529 PetscCall(PetscOptionsBool("-ksp_fetidp_pressure_schur", "Use a BDDC solver for pressure", NULL, schp, &schp, NULL));
530 PetscOptionsEnd();
531
532 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size));
533 fetidp->saddlepoint = (nfp > 0 ? PETSC_TRUE : fetidp->saddlepoint);
534 if (size == 1) fetidp->saddlepoint = PETSC_FALSE;
535
536 PetscCall(KSPGetOperators(ksp, &A, &Ap));
537 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
538 PetscCheck(ismatis, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Amat should be of type MATIS");
539
540 /* Quiet return if the matrix states are unchanged.
541 Needed only for the saddle point case since it uses MatZeroRows
542 on a matrix that may not have changed */
543 PetscCall(PetscObjectStateGet((PetscObject)A, &matstate));
544 PetscCall(MatGetNonzeroState(A, &matnnzstate));
545 if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(PETSC_SUCCESS);
546 fetidp->matstate = matstate;
547 fetidp->matnnzstate = matnnzstate;
548 fetidp->statechanged = fetidp->saddlepoint;
549
550 /* see if we have some fields attached */
551 if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
552 DM dm;
553 PetscContainer c;
554
555 PetscCall(KSPGetDM(ksp, &dm));
556 PetscCall(PetscObjectQuery((PetscObject)A, "_convert_nest_lfields", (PetscObject *)&c));
557 if (dm) {
558 IS *fields;
559 PetscInt nf, i;
560
561 PetscCall(DMCreateFieldDecomposition(dm, &nf, NULL, &fields, NULL));
562 PetscCall(PCBDDCSetDofsSplitting(fetidp->innerbddc, nf, fields));
563 for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
564 PetscCall(PetscFree(fields));
565 } else if (c) {
566 MatISLocalFields lf;
567
568 PetscCall(PetscContainerGetPointer(c, &lf));
569 PetscCall(PCBDDCSetDofsSplittingLocal(fetidp->innerbddc, lf->nr, lf->rf));
570 }
571 }
572
573 if (!fetidp->saddlepoint) {
574 PetscCall(PCSetOperators(fetidp->innerbddc, A, A));
575 } else {
576 Mat nA, lA, PPmat;
577 MatNullSpace nnsp;
578 IS pP;
579 PetscInt totP;
580
581 PetscCall(MatISGetLocalMat(A, &lA));
582 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", (PetscObject)lA));
583
584 pP = fetidp->pP;
585 if (!pP) { /* first time, need to compute pressure dofs */
586 PC_IS *pcis = (PC_IS *)fetidp->innerbddc->data;
587 Mat_IS *matis = (Mat_IS *)A->data;
588 ISLocalToGlobalMapping l2g;
589 IS lP = NULL, II, pII, lPall, Pall, is1, is2;
590 const PetscInt *idxs;
591 PetscInt nl, ni, *widxs;
592 PetscInt i, j, n_neigh, *neigh, *n_shared, **shared, *count;
593 PetscInt rst, ren, n;
594 PetscBool ploc;
595
596 PetscCall(MatGetLocalSize(A, &nl, NULL));
597 PetscCall(MatGetOwnershipRange(A, &rst, &ren));
598 PetscCall(MatGetLocalSize(lA, &n, NULL));
599 PetscCall(MatISGetLocalToGlobalMapping(A, &l2g, NULL));
600
601 if (!pcis->is_I_local) { /* need to compute interior dofs */
602 PetscCall(PetscCalloc1(n, &count));
603 PetscCall(ISLocalToGlobalMappingGetInfo(l2g, &n_neigh, &neigh, &n_shared, &shared));
604 for (i = 1; i < n_neigh; i++)
605 for (j = 0; j < n_shared[i]; j++) count[shared[i][j]] += 1;
606 for (i = 0, j = 0; i < n; i++)
607 if (!count[i]) count[j++] = i;
608 PetscCall(ISLocalToGlobalMappingRestoreInfo(l2g, &n_neigh, &neigh, &n_shared, &shared));
609 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, count, PETSC_OWN_POINTER, &II));
610 } else {
611 PetscCall(PetscObjectReference((PetscObject)pcis->is_I_local));
612 II = pcis->is_I_local;
613 }
614
615 /* interior dofs in layout */
616 PetscCall(PetscArrayzero(matis->sf_leafdata, n));
617 PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
618 PetscCall(ISGetLocalSize(II, &ni));
619 PetscCall(ISGetIndices(II, &idxs));
620 for (i = 0; i < ni; i++) matis->sf_leafdata[idxs[i]] = 1;
621 PetscCall(ISRestoreIndices(II, &idxs));
622 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
623 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
624 PetscCall(PetscMalloc1(PetscMax(nl, n), &widxs));
625 for (i = 0, ni = 0; i < nl; i++)
626 if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
627 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &pII));
628
629 /* pressure dofs */
630 Pall = NULL;
631 lPall = NULL;
632 ploc = PETSC_FALSE;
633 if (nfp == 0) { /* zero pressure block */
634 PetscInt np;
635
636 PetscCall(MatFindZeroDiagonals(A, &Pall));
637 PetscCall(ISGetSize(Pall, &np));
638 if (!np) { /* zero-block not found, defaults to last field (if set) */
639 nfp = 1;
640 fidp[0] = pcbddc->n_ISForDofsLocal ? pcbddc->n_ISForDofsLocal - 1 : pcbddc->n_ISForDofs - 1;
641 PetscCall(ISDestroy(&Pall));
642 } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
643 PetscCall(PCBDDCSetDofsSplitting(fetidp->innerbddc, 1, &Pall));
644 }
645 }
646 if (!Pall) { /* look for registered fields when no zero block has been found */
647 IS *tis;
648
649 PetscCall(PetscMalloc1(nfp, &tis));
650 if (pcbddc->n_ISForDofsLocal) {
651 for (PetscInt i = 0; i < nfp; i++) {
652 PetscInt fid = fidp[i];
653
654 PetscCheck(fid >= 0 && fid < pcbddc->n_ISForDofsLocal, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Invalid field id for pressure %" PetscInt_FMT ", max %" PetscInt_FMT, fid, pcbddc->n_ISForDofsLocal);
655 /* need a sequential IS */
656 PetscCall(ISOnComm(pcbddc->ISForDofsLocal[fid], PETSC_COMM_SELF, PETSC_COPY_VALUES, &tis[i]));
657 }
658 PetscCall(ISConcatenate(PETSC_COMM_SELF, nfp, tis, &lPall));
659 ploc = PETSC_TRUE;
660 } else if (pcbddc->n_ISForDofs) {
661 for (PetscInt i = 0; i < nfp; i++) {
662 PetscInt fid = fidp[i];
663
664 PetscCheck(fid >= 0 && fid < pcbddc->n_ISForDofs, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Invalid field id for pressure %" PetscInt_FMT ", max %" PetscInt_FMT, fid, pcbddc->n_ISForDofs);
665 PetscCall(PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]));
666 tis[i] = pcbddc->ISForDofs[fid];
667 }
668 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), nfp, tis, &Pall));
669 PetscCall(ISSort(Pall));
670 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal");
671 for (PetscInt i = 0; i < nfp; i++) PetscCall(ISDestroy(&tis[i]));
672 PetscCall(PetscFree(tis));
673 }
674
675 /* if the user requested the entire pressure,
676 remove the interior pressure dofs from II (or pII) */
677 if (allp) {
678 if (ploc) {
679 IS nII;
680 PetscCall(ISDifference(II, lPall, &nII));
681 PetscCall(ISDestroy(&II));
682 II = nII;
683 } else {
684 IS nII;
685 PetscCall(ISDifference(pII, Pall, &nII));
686 PetscCall(ISDestroy(&pII));
687 pII = nII;
688 }
689 }
690 if (ploc) {
691 PetscCall(ISDifference(lPall, II, &lP));
692 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP", (PetscObject)lP));
693 } else {
694 PetscCall(ISDifference(Pall, pII, &pP));
695 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", (PetscObject)pP));
696 /* need all local pressure dofs */
697 PetscCall(PetscArrayzero(matis->sf_leafdata, n));
698 PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
699 PetscCall(ISGetLocalSize(Pall, &ni));
700 PetscCall(ISGetIndices(Pall, &idxs));
701 for (i = 0; i < ni; i++) matis->sf_rootdata[idxs[i] - rst] = 1;
702 PetscCall(ISRestoreIndices(Pall, &idxs));
703 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
704 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
705 for (i = 0, ni = 0; i < n; i++)
706 if (matis->sf_leafdata[i]) widxs[ni++] = i;
707 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ni, widxs, PETSC_COPY_VALUES, &lPall));
708 }
709
710 if (!Pall) {
711 PetscCall(PetscArrayzero(matis->sf_leafdata, n));
712 PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
713 PetscCall(ISGetLocalSize(lPall, &ni));
714 PetscCall(ISGetIndices(lPall, &idxs));
715 for (i = 0; i < ni; i++) matis->sf_leafdata[idxs[i]] = 1;
716 PetscCall(ISRestoreIndices(lPall, &idxs));
717 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
718 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
719 for (i = 0, ni = 0; i < nl; i++)
720 if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
721 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &Pall));
722 }
723 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject)Pall));
724
725 if (flip) {
726 PetscInt npl;
727 PetscCall(ISGetLocalSize(Pall, &npl));
728 PetscCall(ISGetIndices(Pall, &idxs));
729 PetscCall(MatCreateVecs(A, NULL, &fetidp->rhs_flip));
730 PetscCall(VecSet(fetidp->rhs_flip, 1.));
731 PetscCall(VecSetOption(fetidp->rhs_flip, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
732 for (i = 0; i < npl; i++) PetscCall(VecSetValue(fetidp->rhs_flip, idxs[i], -1., INSERT_VALUES));
733 PetscCall(VecAssemblyBegin(fetidp->rhs_flip));
734 PetscCall(VecAssemblyEnd(fetidp->rhs_flip));
735 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_flip", (PetscObject)fetidp->rhs_flip));
736 PetscCall(ISRestoreIndices(Pall, &idxs));
737 }
738 PetscCall(ISDestroy(&Pall));
739 PetscCall(ISDestroy(&pII));
740
741 /* local selected pressures in subdomain-wise and global ordering */
742 PetscCall(PetscArrayzero(matis->sf_leafdata, n));
743 PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
744 if (!ploc) {
745 PetscInt *widxs2;
746
747 PetscCheck(pP, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Missing parallel pressure IS");
748 PetscCall(ISGetLocalSize(pP, &ni));
749 PetscCall(ISGetIndices(pP, &idxs));
750 for (i = 0; i < ni; i++) matis->sf_rootdata[idxs[i] - rst] = 1;
751 PetscCall(ISRestoreIndices(pP, &idxs));
752 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
753 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
754 for (i = 0, ni = 0; i < n; i++)
755 if (matis->sf_leafdata[i]) widxs[ni++] = i;
756 PetscCall(PetscMalloc1(ni, &widxs2));
757 PetscCall(ISLocalToGlobalMappingApply(l2g, ni, widxs, widxs2));
758 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ni, widxs, PETSC_COPY_VALUES, &lP));
759 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP", (PetscObject)lP));
760 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs2, PETSC_OWN_POINTER, &is1));
761 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_gP", (PetscObject)is1));
762 PetscCall(ISDestroy(&is1));
763 } else {
764 PetscCheck(lP, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing sequential pressure IS");
765 PetscCall(ISGetLocalSize(lP, &ni));
766 PetscCall(ISGetIndices(lP, &idxs));
767 for (i = 0; i < ni; i++)
768 if (idxs[i] >= 0 && idxs[i] < n) matis->sf_leafdata[idxs[i]] = 1;
769 PetscCall(ISRestoreIndices(lP, &idxs));
770 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
771 PetscCall(ISLocalToGlobalMappingApply(l2g, ni, idxs, widxs));
772 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &is1));
773 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_gP", (PetscObject)is1));
774 PetscCall(ISDestroy(&is1));
775 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
776 for (i = 0, ni = 0; i < nl; i++)
777 if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
778 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &pP));
779 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", (PetscObject)pP));
780 }
781 PetscCall(PetscFree(widxs));
782
783 /* If there's any "interior pressure",
784 we may want to use a discrete harmonic solver instead
785 of a Stokes harmonic for the Dirichlet preconditioner
786 Need to extract the interior velocity dofs in interior dofs ordering (iV)
787 and interior pressure dofs in local ordering (iP) */
788 if (!allp) {
789 ISLocalToGlobalMapping l2g_t;
790
791 PetscCall(ISDifference(lPall, lP, &is1));
792 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_iP", (PetscObject)is1));
793 PetscCall(ISDifference(II, is1, &is2));
794 PetscCall(ISDestroy(&is1));
795 PetscCall(ISLocalToGlobalMappingCreateIS(II, &l2g_t));
796 PetscCall(ISGlobalToLocalMappingApplyIS(l2g_t, IS_GTOLM_DROP, is2, &is1));
797 PetscCall(ISGetLocalSize(is1, &i));
798 PetscCall(ISGetLocalSize(is2, &j));
799 PetscCheck(i == j, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local sizes %" PetscInt_FMT " and %" PetscInt_FMT " for iV", i, j);
800 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_iV", (PetscObject)is1));
801 PetscCall(ISLocalToGlobalMappingDestroy(&l2g_t));
802 PetscCall(ISDestroy(&is1));
803 PetscCall(ISDestroy(&is2));
804 }
805
806 /* exclude selected pressures from the inner BDDC */
807 if (pcbddc->DirichletBoundariesLocal) {
808 IS list[2], plP, isout;
809
810 /* need a parallel IS */
811 PetscCall(ISOnComm(lP, PetscObjectComm((PetscObject)ksp), PETSC_COPY_VALUES, &plP));
812 list[0] = plP;
813 list[1] = pcbddc->DirichletBoundariesLocal;
814 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), 2, list, &isout));
815 PetscCall(ISSortRemoveDups(isout));
816 PetscCall(ISDestroy(&plP));
817 PetscCall(PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc, isout));
818 PetscCall(ISDestroy(&isout));
819 } else if (pcbddc->DirichletBoundaries) {
820 IS list[2], isout;
821
822 list[0] = pP;
823 list[1] = pcbddc->DirichletBoundaries;
824 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), 2, list, &isout));
825 PetscCall(ISSortRemoveDups(isout));
826 PetscCall(PCBDDCSetDirichletBoundaries(fetidp->innerbddc, isout));
827 PetscCall(ISDestroy(&isout));
828 } else {
829 IS plP;
830
831 /* need a parallel IS */
832 PetscCall(ISOnComm(lP, PetscObjectComm((PetscObject)ksp), PETSC_COPY_VALUES, &plP));
833 PetscCall(PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc, plP));
834 PetscCall(ISDestroy(&plP));
835 }
836
837 /* Need to zero output of interface BDDC for lP */
838 {
839 IS BB, lP_I, lP_B;
840
841 PetscCall(ISComplement(II, 0, n, &BB));
842 PetscCall(ISEmbed(lP, II, PETSC_TRUE, &lP_I));
843 PetscCall(ISEmbed(lP, BB, PETSC_TRUE, &lP_B));
844 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP_I", (PetscObject)lP_I));
845 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP_B", (PetscObject)lP_B));
846 PetscCall(ISDestroy(&BB));
847 PetscCall(ISDestroy(&lP_I));
848 PetscCall(ISDestroy(&lP_B));
849 }
850 PetscCall(ISDestroy(&II));
851
852 /* save CSR information for the pressure BDDC solver (if any) */
853 if (schp) {
854 PetscInt np, nt;
855
856 PetscCall(MatGetSize(matis->A, &nt, NULL));
857 PetscCall(ISGetLocalSize(lP, &np));
858 if (np) {
859 PetscInt *xadj = pcbddc->mat_graph->xadj;
860 PetscInt *adjn = pcbddc->mat_graph->adjncy;
861 PetscInt nv = pcbddc->mat_graph->nvtxs_csr;
862
863 if (nv && nv == nt) {
864 ISLocalToGlobalMapping pmap;
865 PetscInt *schp_csr, *schp_xadj, *schp_adjn, p;
866 PetscContainer c;
867
868 PetscCall(ISLocalToGlobalMappingCreateIS(lPall, &pmap));
869 PetscCall(ISGetIndices(lPall, &idxs));
870 for (p = 0, nv = 0; p < np; p++) {
871 PetscInt x, n = idxs[p];
872
873 PetscCall(ISGlobalToLocalMappingApply(pmap, IS_GTOLM_DROP, xadj[n + 1] - xadj[n], adjn + xadj[n], &x, NULL));
874 nv += x;
875 }
876 PetscCall(PetscMalloc1(np + 1 + nv, &schp_csr));
877 schp_xadj = schp_csr;
878 schp_adjn = schp_csr + np + 1;
879 for (p = 0, schp_xadj[0] = 0; p < np; p++) {
880 PetscInt x, n = idxs[p];
881
882 PetscCall(ISGlobalToLocalMappingApply(pmap, IS_GTOLM_DROP, xadj[n + 1] - xadj[n], adjn + xadj[n], &x, schp_adjn + schp_xadj[p]));
883 schp_xadj[p + 1] = schp_xadj[p] + x;
884 }
885 PetscCall(ISRestoreIndices(lPall, &idxs));
886 PetscCall(ISLocalToGlobalMappingDestroy(&pmap));
887 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
888 PetscCall(PetscContainerSetPointer(c, schp_csr));
889 PetscCall(PetscContainerSetCtxDestroy(c, PetscCtxDestroyDefault));
890 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pCSR", (PetscObject)c));
891 PetscCall(PetscContainerDestroy(&c));
892 }
893 }
894 }
895 PetscCall(ISDestroy(&lPall));
896 PetscCall(ISDestroy(&lP));
897 fetidp->pP = pP;
898 }
899
900 /* total number of selected pressure dofs */
901 PetscCall(ISGetSize(fetidp->pP, &totP));
902
903 /* Set operator for inner BDDC */
904 if (totP || fetidp->rhs_flip) {
905 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &nA));
906 } else {
907 PetscCall(PetscObjectReference((PetscObject)A));
908 nA = A;
909 }
910 if (fetidp->rhs_flip) {
911 PetscCall(MatDiagonalScale(nA, fetidp->rhs_flip, NULL));
912 if (totP) {
913 Mat lA2;
914
915 PetscCall(MatISGetLocalMat(nA, &lA));
916 PetscCall(MatDuplicate(lA, MAT_COPY_VALUES, &lA2));
917 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", (PetscObject)lA2));
918 PetscCall(MatDestroy(&lA2));
919 }
920 }
921
922 /* assign operator to compute inner bddc */
923 if (totP) {
924 /* in this case, we remove all the used pressure couplings */
925 PetscCall(MatSetOption(nA, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
926 PetscCall(MatZeroRowsColumnsIS(nA, fetidp->pP, 1., NULL, NULL));
927 } else {
928 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", NULL));
929 }
930 PetscCall(MatGetNearNullSpace(Ap, &nnsp));
931 if (!nnsp) PetscCall(MatGetNullSpace(Ap, &nnsp));
932 if (!nnsp) PetscCall(MatGetNearNullSpace(A, &nnsp));
933 if (!nnsp) PetscCall(MatGetNullSpace(A, &nnsp));
934 PetscCall(MatSetNearNullSpace(nA, nnsp));
935 PetscCall(PCSetOperators(fetidp->innerbddc, nA, nA));
936 PetscCall(MatDestroy(&nA));
937
938 /* non-zero rhs on interior dofs when applying the preconditioner */
939 if (totP) pcbddc->switch_static = PETSC_TRUE;
940
941 /* if there are no interface pressures, set inner bddc flag for benign saddle point */
942 if (!totP) {
943 pcbddc->benign_saddle_point = PETSC_TRUE;
944 pcbddc->compute_nonetflux = PETSC_TRUE;
945 }
946
947 /* Operators for pressure preconditioner */
948 if (totP) {
949 /* Extract pressure block if needed */
950 if (!pisz) {
951 Mat C;
952 IS nzrows = NULL;
953
954 PetscCall(MatCreateSubMatrix(A, fetidp->pP, fetidp->pP, MAT_INITIAL_MATRIX, &C));
955 PetscCall(MatFindNonzeroRows(C, &nzrows));
956 if (nzrows) {
957 PetscInt i;
958
959 PetscCall(ISGetSize(nzrows, &i));
960 PetscCall(ISDestroy(&nzrows));
961 if (!i) pisz = PETSC_TRUE;
962 }
963 if (!pisz) {
964 PetscCall(MatScale(C, -1.)); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized, Biot... */
965 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_C", (PetscObject)C));
966 }
967 PetscCall(MatDestroy(&C));
968 }
969 /* Divergence mat */
970 if (!pcbddc->divudotp) {
971 Mat B;
972 IS P;
973 IS l2l = NULL;
974 PetscBool save;
975
976 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&P));
977 if (!pisz) {
978 IS F, V, Ps;
979 PetscInt m, M;
980
981 PetscCall(MatGetOwnershipRange(A, &m, &M));
982 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), M - m, m, 1, &F));
983 PetscCall(ISDuplicate(P, &Ps));
984 PetscCall(ISSort(Ps));
985 PetscCall(ISComplement(Ps, m, M, &V));
986 PetscCall(ISDestroy(&Ps));
987 PetscCall(MatCreateSubMatrix(A, P, V, MAT_INITIAL_MATRIX, &B));
988 {
989 Mat_IS *Bmatis = (Mat_IS *)B->data;
990 PetscCall(PetscObjectReference((PetscObject)Bmatis->getsub_cis));
991 l2l = Bmatis->getsub_cis;
992 }
993 PetscCall(ISDestroy(&V));
994 PetscCall(ISDestroy(&F));
995 } else {
996 PetscCall(MatCreateSubMatrix(A, P, NULL, MAT_INITIAL_MATRIX, &B));
997 }
998 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */
999 PetscCall(PCBDDCSetDivergenceMat(fetidp->innerbddc, B, PETSC_FALSE, l2l));
1000 pcbddc->compute_nonetflux = save;
1001 PetscCall(MatDestroy(&B));
1002 PetscCall(ISDestroy(&l2l));
1003 }
1004 if (A != Ap) { /* user has provided a different Pmat, this always supersedes the setter (TODO: is it OK?) */
1005 /* use monolithic operator, we restrict later */
1006 PetscCall(KSPFETIDPSetPressureOperator(ksp, Ap));
1007 }
1008 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject *)&PPmat));
1009
1010 /* PPmat not present, use some default choice */
1011 if (!PPmat) {
1012 Mat C;
1013
1014 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_C", (PetscObject *)&C));
1015 if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */
1016 PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
1017 } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */
1018 IS P;
1019
1020 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&P));
1021 PetscCall(MatCreateSubMatrix(A, P, P, MAT_INITIAL_MATRIX, &C));
1022 PetscCall(MatScale(C, -1.));
1023 PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
1024 PetscCall(MatDestroy(&C));
1025 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */
1026 PetscInt nl;
1027
1028 PetscCall(ISGetLocalSize(fetidp->pP, &nl));
1029 PetscCall(MatCreate(PetscObjectComm((PetscObject)ksp), &C));
1030 PetscCall(MatSetSizes(C, nl, nl, totP, totP));
1031 PetscCall(MatSetType(C, MATAIJ));
1032 PetscCall(MatMPIAIJSetPreallocation(C, 1, NULL, 0, NULL));
1033 PetscCall(MatSeqAIJSetPreallocation(C, 1, NULL));
1034 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1035 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1036 PetscCall(MatShift(C, 1.));
1037 PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
1038 PetscCall(MatDestroy(&C));
1039 }
1040 }
1041
1042 /* Preconditioned operator for the pressure block */
1043 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject *)&PPmat));
1044 if (PPmat) {
1045 Mat C;
1046 IS Pall;
1047 PetscInt AM, PAM, PAN, pam, pan, am, an, pl, pIl, pAg, pIg;
1048
1049 PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&Pall));
1050 PetscCall(MatGetSize(A, &AM, NULL));
1051 PetscCall(MatGetSize(PPmat, &PAM, &PAN));
1052 PetscCall(ISGetSize(Pall, &pAg));
1053 PetscCall(ISGetSize(fetidp->pP, &pIg));
1054 PetscCall(MatGetLocalSize(PPmat, &pam, &pan));
1055 PetscCall(MatGetLocalSize(A, &am, &an));
1056 PetscCall(ISGetLocalSize(Pall, &pIl));
1057 PetscCall(ISGetLocalSize(fetidp->pP, &pl));
1058 PetscCheck(PAM == PAN, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Pressure matrix must be square, unsupported %" PetscInt_FMT " x %" PetscInt_FMT, PAM, PAN);
1059 PetscCheck(pam == pan, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Local sizes of pressure matrix must be equal, unsupported %" PetscInt_FMT " x %" PetscInt_FMT, pam, pan);
1060 PetscCheck(pam == am || pam == pl || pam == pIl, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid number of local rows %" PetscInt_FMT " for pressure matrix! Supported are %" PetscInt_FMT ", %" PetscInt_FMT " or %" PetscInt_FMT, pam, am, pl, pIl);
1061 PetscCheck(pan == an || pan == pl || pan == pIl, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid number of local columns %" PetscInt_FMT " for pressure matrix! Supported are %" PetscInt_FMT ", %" PetscInt_FMT " or %" PetscInt_FMT, pan, an, pl, pIl);
1062 if (PAM == AM) { /* monolithic ordering, restrict to pressure */
1063 if (schp) {
1064 PetscCall(MatCreateSubMatrix(PPmat, Pall, Pall, MAT_INITIAL_MATRIX, &C));
1065 PetscCall(MatNullSpacePropagateAny_Private(PPmat, Pall, C));
1066 } else {
1067 PetscCall(MatCreateSubMatrix(PPmat, fetidp->pP, fetidp->pP, MAT_INITIAL_MATRIX, &C));
1068 }
1069 } else if (pAg == PAM) { /* global ordering for pressure only */
1070 if (!allp && !schp) { /* solving for interface pressure only */
1071 IS restr;
1072
1073 PetscCall(ISRenumber(fetidp->pP, NULL, NULL, &restr));
1074 PetscCall(MatCreateSubMatrix(PPmat, restr, restr, MAT_INITIAL_MATRIX, &C));
1075 PetscCall(ISDestroy(&restr));
1076 } else {
1077 PetscCall(PetscObjectReference((PetscObject)PPmat));
1078 C = PPmat;
1079 }
1080 } else if (pIg == PAM) { /* global ordering for selected pressure only */
1081 PetscCheck(!schp, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Need the entire matrix");
1082 PetscCall(PetscObjectReference((PetscObject)PPmat));
1083 C = PPmat;
1084 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Unable to use the pressure matrix");
1085
1086 PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
1087 PetscCall(MatDestroy(&C));
1088 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Missing Pmat for pressure block");
1089 } else { /* totP == 0 */
1090 PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", NULL));
1091 }
1092 }
1093 PetscFunctionReturn(PETSC_SUCCESS);
1094 }
1095
KSPSetUp_FETIDP(KSP ksp)1096 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp)
1097 {
1098 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1099 PC_BDDC *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
1100 PetscBool flg;
1101
1102 PetscFunctionBegin;
1103 PetscCall(KSPFETIDPSetUpOperators(ksp));
1104 /* set up BDDC */
1105 PetscCall(PCSetErrorIfFailure(fetidp->innerbddc, ksp->errorifnotconverged));
1106 PetscCall(PCSetUp(fetidp->innerbddc));
1107 /* FETI-DP as it is implemented needs an exact coarse solver */
1108 if (pcbddc->coarse_ksp) {
1109 PetscCall(KSPSetTolerances(pcbddc->coarse_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, 1000));
1110 PetscCall(KSPSetNormType(pcbddc->coarse_ksp, KSP_NORM_DEFAULT));
1111 }
1112 /* FETI-DP as it is implemented needs exact local Neumann solvers */
1113 PetscCall(KSPSetTolerances(pcbddc->ksp_R, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, 1000));
1114 PetscCall(KSPSetNormType(pcbddc->ksp_R, KSP_NORM_DEFAULT));
1115
1116 /* setup FETI-DP operators
1117 If fetidp->statechanged is true, we need to update the operators
1118 needed in the saddle-point case. This should be replaced
1119 by a better logic when the FETI-DP matrix and preconditioner will
1120 have their own classes */
1121 if (pcbddc->new_primal_space || fetidp->statechanged) {
1122 Mat F; /* the FETI-DP matrix */
1123 PC D; /* the FETI-DP preconditioner */
1124 PetscCall(KSPReset(fetidp->innerksp));
1125 PetscCall(PCBDDCCreateFETIDPOperators(fetidp->innerbddc, fetidp->fully_redundant, ((PetscObject)ksp)->prefix, &F, &D));
1126 PetscCall(KSPSetOperators(fetidp->innerksp, F, F));
1127 PetscCall(KSPSetTolerances(fetidp->innerksp, ksp->rtol, ksp->abstol, ksp->divtol, ksp->max_it));
1128 PetscCall(KSPSetPC(fetidp->innerksp, D));
1129 PetscCall(PetscObjectIncrementTabLevel((PetscObject)D, (PetscObject)fetidp->innerksp, 0));
1130 PetscCall(KSPSetFromOptions(fetidp->innerksp));
1131 PetscCall(MatCreateVecs(F, &fetidp->innerksp->vec_rhs, &fetidp->innerksp->vec_sol));
1132 PetscCall(MatDestroy(&F));
1133 PetscCall(PCDestroy(&D));
1134 if (fetidp->check) {
1135 PetscViewer viewer;
1136
1137 if (!pcbddc->dbg_viewer) {
1138 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
1139 } else {
1140 viewer = pcbddc->dbg_viewer;
1141 }
1142 PetscCall(KSPFETIDPCheckOperators(ksp, viewer));
1143 }
1144 }
1145 fetidp->statechanged = PETSC_FALSE;
1146 pcbddc->new_primal_space = PETSC_FALSE;
1147
1148 /* propagate settings to the inner solve */
1149 PetscCall(KSPGetComputeSingularValues(ksp, &flg));
1150 PetscCall(KSPSetComputeSingularValues(fetidp->innerksp, flg));
1151 if (ksp->res_hist) PetscCall(KSPSetResidualHistory(fetidp->innerksp, ksp->res_hist, ksp->res_hist_max, ksp->res_hist_reset));
1152 PetscCall(KSPSetErrorIfNotConverged(fetidp->innerksp, ksp->errorifnotconverged));
1153 PetscCall(KSPSetUp(fetidp->innerksp));
1154 PetscFunctionReturn(PETSC_SUCCESS);
1155 }
1156
KSPSolve_FETIDP(KSP ksp)1157 static PetscErrorCode KSPSolve_FETIDP(KSP ksp)
1158 {
1159 Mat F, A;
1160 MatNullSpace nsp;
1161 Vec X, B, Xl, Bl;
1162 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1163 PC_BDDC *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
1164 KSPConvergedReason reason;
1165 PC pc;
1166 PCFailedReason pcreason;
1167 PetscInt hist_len;
1168 int flg;
1169
1170 PetscFunctionBegin;
1171 PetscCall(PetscCitationsRegister(citation, &cited));
1172 if (fetidp->saddlepoint) PetscCall(PetscCitationsRegister(citation2, &cited2));
1173 PetscCall(KSPGetOperators(ksp, &A, NULL));
1174 PetscCall(KSPGetRhs(ksp, &B));
1175 PetscCall(KSPGetSolution(ksp, &X));
1176 PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
1177 PetscCall(KSPGetRhs(fetidp->innerksp, &Bl));
1178 PetscCall(KSPGetSolution(fetidp->innerksp, &Xl));
1179 PetscCall(PCBDDCMatFETIDPGetRHS(F, B, Bl));
1180 if (ksp->transpose_solve) {
1181 PetscCall(KSPSolveTranspose(fetidp->innerksp, Bl, Xl));
1182 } else {
1183 PetscCall(KSPSolve(fetidp->innerksp, Bl, Xl));
1184 }
1185 PetscCall(KSPGetConvergedReason(fetidp->innerksp, &reason));
1186 PetscCall(KSPGetPC(fetidp->innerksp, &pc));
1187 PetscCall(PCGetFailedReason(pc, &pcreason));
1188 flg = (reason < 0 && reason != KSP_DIVERGED_ITS) || pcreason;
1189 PetscCall(VecFlag(Xl, flg));
1190 if (flg) {
1191 PetscInt its;
1192
1193 PetscCall(KSPGetIterationNumber(fetidp->innerksp, &its));
1194 ksp->reason = KSP_DIVERGED_PC_FAILED;
1195 PetscCall(PetscInfo(ksp, "Inner KSP solve failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
1196 }
1197 PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, X));
1198 PetscCall(MatGetNullSpace(A, &nsp));
1199 if (nsp) PetscCall(MatNullSpaceRemove(nsp, X));
1200 /* update ksp with stats from inner ksp */
1201 PetscCall(KSPGetConvergedReason(fetidp->innerksp, &ksp->reason));
1202 PetscCall(KSPGetIterationNumber(fetidp->innerksp, &ksp->its));
1203 ksp->totalits += ksp->its;
1204 PetscCall(KSPGetResidualHistory(fetidp->innerksp, NULL, &hist_len));
1205 ksp->res_hist_len = (size_t)hist_len;
1206 /* restore defaults for inner BDDC (Pre/PostSolve flags) */
1207 pcbddc->temp_solution_used = PETSC_FALSE;
1208 pcbddc->rhs_change = PETSC_FALSE;
1209 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1210 PetscFunctionReturn(PETSC_SUCCESS);
1211 }
1212
KSPReset_FETIDP(KSP ksp)1213 static PetscErrorCode KSPReset_FETIDP(KSP ksp)
1214 {
1215 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1216 PC_BDDC *pcbddc;
1217
1218 PetscFunctionBegin;
1219 PetscCall(ISDestroy(&fetidp->pP));
1220 PetscCall(VecDestroy(&fetidp->rhs_flip));
1221 /* avoid PCReset that does not take into account ref counting */
1222 PetscCall(PCDestroy(&fetidp->innerbddc));
1223 PetscCall(PCCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerbddc));
1224 PetscCall(PCSetType(fetidp->innerbddc, PCBDDC));
1225 pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
1226 pcbddc->symmetric_primal = PETSC_FALSE;
1227 PetscCall(KSPDestroy(&fetidp->innerksp));
1228 fetidp->saddlepoint = PETSC_FALSE;
1229 fetidp->matstate = -1;
1230 fetidp->matnnzstate = -1;
1231 fetidp->statechanged = PETSC_TRUE;
1232 PetscFunctionReturn(PETSC_SUCCESS);
1233 }
1234
KSPDestroy_FETIDP(KSP ksp)1235 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
1236 {
1237 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1238
1239 PetscFunctionBegin;
1240 PetscCall(KSPReset_FETIDP(ksp));
1241 PetscCall(PCDestroy(&fetidp->innerbddc));
1242 PetscCall(KSPDestroy(&fetidp->innerksp));
1243 PetscCall(PetscFree(fetidp->monctx));
1244 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetInnerBDDC_C", NULL));
1245 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerBDDC_C", NULL));
1246 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerKSP_C", NULL));
1247 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetPressureOperator_C", NULL));
1248 PetscCall(PetscFree(ksp->data));
1249 PetscFunctionReturn(PETSC_SUCCESS);
1250 }
1251
KSPView_FETIDP(KSP ksp,PetscViewer viewer)1252 static PetscErrorCode KSPView_FETIDP(KSP ksp, PetscViewer viewer)
1253 {
1254 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1255 PetscBool isascii;
1256
1257 PetscFunctionBegin;
1258 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1259 if (isascii) {
1260 PetscCall(PetscViewerASCIIPrintf(viewer, " fully redundant: %d\n", fetidp->fully_redundant));
1261 PetscCall(PetscViewerASCIIPrintf(viewer, " saddle point: %d\n", fetidp->saddlepoint));
1262 PetscCall(PetscViewerASCIIPrintf(viewer, "Inner KSP solver details\n"));
1263 }
1264 PetscCall(PetscViewerASCIIPushTab(viewer));
1265 PetscCall(KSPView(fetidp->innerksp, viewer));
1266 PetscCall(PetscViewerASCIIPopTab(viewer));
1267 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Inner BDDC solver details\n"));
1268 PetscCall(PetscViewerASCIIPushTab(viewer));
1269 PetscCall(PCView(fetidp->innerbddc, viewer));
1270 PetscCall(PetscViewerASCIIPopTab(viewer));
1271 PetscFunctionReturn(PETSC_SUCCESS);
1272 }
1273
KSPSetFromOptions_FETIDP(KSP ksp,PetscOptionItems PetscOptionsObject)1274 static PetscErrorCode KSPSetFromOptions_FETIDP(KSP ksp, PetscOptionItems PetscOptionsObject)
1275 {
1276 KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
1277
1278 PetscFunctionBegin;
1279 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
1280 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp, ((PetscObject)ksp)->prefix));
1281 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp, "fetidp_"));
1282 if (!fetidp->userbddc) {
1283 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc, ((PetscObject)ksp)->prefix));
1284 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc, "fetidp_bddc_"));
1285 }
1286 PetscOptionsHeadBegin(PetscOptionsObject, "KSP FETIDP options");
1287 PetscCall(PetscOptionsBool("-ksp_fetidp_fullyredundant", "Use fully redundant multipliers", "none", fetidp->fully_redundant, &fetidp->fully_redundant, NULL));
1288 PetscCall(PetscOptionsBool("-ksp_fetidp_saddlepoint", "Activates support for saddle-point problems", NULL, fetidp->saddlepoint, &fetidp->saddlepoint, NULL));
1289 PetscCall(PetscOptionsBool("-ksp_fetidp_check", "Activates verbose debugging output FETI-DP operators", NULL, fetidp->check, &fetidp->check, NULL));
1290 PetscOptionsHeadEnd();
1291 PetscCall(PCSetFromOptions(fetidp->innerbddc));
1292 PetscFunctionReturn(PETSC_SUCCESS);
1293 }
1294
1295 /*MC
1296 KSPFETIDP - The FETI-DP method {cite}`farhat2001feti`
1297
1298 Options Database Keys:
1299 + -ksp_fetidp_fullyredundant <false> - use a fully redundant set of Lagrange multipliers
1300 . -ksp_fetidp_saddlepoint <false> - activates support for saddle point problems, see {cite}`tu2015feti`
1301 . -ksp_fetidp_saddlepoint_flip <false> - see note below
1302 . -ksp_fetidp_pressure_field <-1> - activates support for saddle point problems, and identifies the pressure field id.
1303 If this information is not provided, the pressure field is detected by using `MatFindZeroDiagonals()`.
1304 - -ksp_fetidp_pressure_all <false> - if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field.
1305
1306 Level: Advanced
1307
1308 Notes:
1309 The matrix for the `KSP` must be of type `MATIS`.
1310
1311 Usually, an incompressible Stokes problem is written as
1312 .vb
1313 | A B^T | | v | = | f |
1314 | B 0 | | p | = | g |
1315 .ve
1316 with B representing $ -\int_\Omega \nabla \cdot u q $. If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as
1317 .vb
1318 | A B^T | | v | = | f |
1319 |-B 0 | | p | = |-g |
1320 .ve
1321
1322 The FETI-DP linear system (automatically generated constructing an internal `PCBDDC` object) is solved using an internal `KSP` object.
1323
1324 Options for the inner `KSP` and for the customization of the `PCBDDC` object can be specified at command line by using the prefixes `-fetidp_` and `-fetidp_bddc_`. E.g.,
1325 .vb
1326 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false
1327 .ve
1328 will use `KSPGMRES` for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric `PCBDDC`.
1329
1330 For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with `KSPFETIDPSetPressureOperator()`.
1331 Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix).
1332 If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver.
1333 Options for the pressure solver can be prefixed with `-fetidp_fielsplit_p_`, E.g.
1334 .vb
1335 -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps
1336 .ve
1337 In order to use the deluxe version of FETI-DP, you must customize the inner `PCBDDC` operator with -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_deluxe_singlemat and use
1338 non-redundant multipliers, i.e. `-ksp_fetidp_fullyredundant false`. Options for the scaling solver are prefixed by `-fetidp_bddelta_`, E.g.
1339 .vb
1340 -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu
1341 .ve
1342
1343 Some of the basic options such as the maximum number of iterations and tolerances are automatically passed from this `KSP` to the inner `KSP` that actually performs the iterations.
1344
1345 The converged reason and number of iterations computed are passed from the inner `KSP` to this `KSP` at the end of the solution.
1346
1347 Developer Note:
1348 Even though this method does not directly use any norms, the user is allowed to set the `KSPNormType` to any value.
1349 This is so users do not have to change `KSPNormType` options when they switch from other `KSP` methods to this one.
1350
1351 .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
1352 M*/
KSPCreate_FETIDP(KSP ksp)1353 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
1354 {
1355 KSP_FETIDP *fetidp;
1356 KSP_FETIDPMon *monctx;
1357 PC_BDDC *pcbddc;
1358 PC pc;
1359
1360 PetscFunctionBegin;
1361 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 3));
1362 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 2));
1363 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1364 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_RIGHT, 2));
1365 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
1366 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
1367 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1368
1369 PetscCall(PetscNew(&fetidp));
1370 fetidp->matstate = -1;
1371 fetidp->matnnzstate = -1;
1372 fetidp->statechanged = PETSC_TRUE;
1373
1374 ksp->data = (void *)fetidp;
1375 ksp->ops->setup = KSPSetUp_FETIDP;
1376 ksp->ops->solve = KSPSolve_FETIDP;
1377 ksp->ops->destroy = KSPDestroy_FETIDP;
1378 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP;
1379 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP;
1380 ksp->ops->view = KSPView_FETIDP;
1381 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP;
1382 ksp->ops->buildsolution = KSPBuildSolution_FETIDP;
1383 ksp->ops->buildresidual = KSPBuildResidualDefault;
1384 PetscCall(KSPGetPC(ksp, &pc));
1385 PetscCall(PCSetType(pc, PCNONE));
1386 /* create the inner KSP for the Lagrange multipliers */
1387 PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerksp));
1388 PetscCall(KSPGetPC(fetidp->innerksp, &pc));
1389 PetscCall(PCSetType(pc, PCNONE));
1390 /* monitor */
1391 PetscCall(PetscNew(&monctx));
1392 monctx->parentksp = ksp;
1393 fetidp->monctx = monctx;
1394 PetscCall(KSPMonitorSet(fetidp->innerksp, KSPMonitor_FETIDP, fetidp->monctx, NULL));
1395 /* create the inner BDDC */
1396 PetscCall(PCCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerbddc));
1397 PetscCall(PCSetType(fetidp->innerbddc, PCBDDC));
1398 /* make sure we always obtain a consistent FETI-DP matrix
1399 for symmetric problems, the user can always customize it through the command line */
1400 pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
1401 pcbddc->symmetric_primal = PETSC_FALSE;
1402 /* composed functions */
1403 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetInnerBDDC_C", KSPFETIDPSetInnerBDDC_FETIDP));
1404 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerBDDC_C", KSPFETIDPGetInnerBDDC_FETIDP));
1405 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerKSP_C", KSPFETIDPGetInnerKSP_FETIDP));
1406 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetPressureOperator_C", KSPFETIDPSetPressureOperator_FETIDP));
1407 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
1408 ksp->setupnewmatrix = PETSC_TRUE;
1409 PetscFunctionReturn(PETSC_SUCCESS);
1410 }
1411