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: }