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