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