Actual source code: fetidp.c

  1: #include <petsc/private/kspimpl.h>
  2: #include <../src/ksp/pc/impls/bddc/bddc.h>
  3: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
  4: #include <petscdm.h>

  6: static PetscBool  cited  = PETSC_FALSE;
  7: static PetscBool  cited2 = PETSC_FALSE;
  8: static const char citation[] =
  9: "@article{ZampiniPCBDDC,\n"
 10: "author = {Stefano Zampini},\n"
 11: "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
 12: "journal = {SIAM Journal on Scientific Computing},\n"
 13: "volume = {38},\n"
 14: "number = {5},\n"
 15: "pages = {S282-S306},\n"
 16: "year = {2016},\n"
 17: "doi = {10.1137/15M1025785},\n"
 18: "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
 19: "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
 20: "}\n"
 21: "@article{ZampiniDualPrimal,\n"
 22: "author = {Stefano Zampini},\n"
 23: "title = {{D}ual-{P}rimal methods for the cardiac {B}idomain model},\n"
 24: "volume = {24},\n"
 25: "number = {04},\n"
 26: "pages = {667-696},\n"
 27: "year = {2014},\n"
 28: "doi = {10.1142/S0218202513500632},\n"
 29: "URL = {https://www.worldscientific.com/doi/abs/10.1142/S0218202513500632},\n"
 30: "eprint = {https://www.worldscientific.com/doi/pdf/10.1142/S0218202513500632}\n"
 31: "}\n";
 32: static const char citation2[] =
 33: "@article{li2013nonoverlapping,\n"
 34: "title={A nonoverlapping domain decomposition method for incompressible Stokes equations with continuous pressures},\n"
 35: "author={Li, Jing and Tu, Xuemin},\n"
 36: "journal={SIAM Journal on Numerical Analysis},\n"
 37: "volume={51},\n"
 38: "number={2},\n"
 39: "pages={1235--1253},\n"
 40: "year={2013},\n"
 41: "publisher={Society for Industrial and Applied Mathematics}\n"
 42: "}\n";

 44: /*
 45:     This file implements the FETI-DP method in PETSc as part of KSP.
 46: */
 47: typedef struct {
 48:   KSP parentksp;
 49: } KSP_FETIDPMon;

 51: typedef struct {
 52:   KSP              innerksp;         /* the KSP for the Lagrange multipliers */
 53:   PC               innerbddc;        /* the inner BDDC object */
 54:   PetscBool        fully_redundant;  /* true for using a fully redundant set of multipliers */
 55:   PetscBool        userbddc;         /* true if the user provided the PCBDDC object */
 56:   PetscBool        saddlepoint;      /* support for saddle point problems */
 57:   IS               pP;               /* index set for pressure variables */
 58:   Vec              rhs_flip;         /* see KSPFETIDPSetUpOperators */
 59:   KSP_FETIDPMon    *monctx;          /* monitor context, used to pass user defined monitors
 60:                                         in the physical space */
 61:   PetscObjectState matstate;         /* these are needed just in the saddle point case */
 62:   PetscObjectState matnnzstate;      /* where we are going to use MatZeroRows on pmat */
 63:   PetscBool        statechanged;
 64:   PetscBool        check;
 65: } KSP_FETIDP;

 67: static PetscErrorCode KSPFETIDPSetPressureOperator_FETIDP(KSP ksp, Mat P)
 68: {
 69:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

 71:   if (P) fetidp->saddlepoint = PETSC_TRUE;
 72:   PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)P);
 73:   return 0;
 74: }

 76: /*@
 77:  KSPFETIDPSetPressureOperator - Sets the operator used to setup the pressure preconditioner for saddle point FETI-DP.

 79:    Collective on ksp

 81:    Input Parameters:
 82: +  ksp - the FETI-DP Krylov solver
 83: -  P - the linear operator to be preconditioned, usually the mass matrix.

 85:    Level: advanced

 87:    Notes:
 88:     The operator can be either passed in a) monolithic global ordering, b) pressure-only global ordering
 89:           or c) interface pressure ordering (if -ksp_fetidp_pressure_all false).
 90:           In cases b) and c), the pressure ordering of dofs needs to satisfy
 91:              pid_1 < pid_2  iff  gid_1 < gid_2
 92:           where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding
 93:           id in the monolithic global ordering.

 95: .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators
 96: @*/
 97: PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P)
 98: {
101:   PetscTryMethod(ksp,"KSPFETIDPSetPressureOperator_C",(KSP,Mat),(ksp,P));
102:   return 0;
103: }

105: static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp)
106: {
107:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

109:   *innerksp = fetidp->innerksp;
110:   return 0;
111: }

113: /*@
114:  KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers

116:    Input Parameters:
117: +  ksp - the FETI-DP KSP
118: -  innerksp - the KSP for the multipliers

120:    Level: advanced

122:    Notes:

124: .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC
125: @*/
126: PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp)
127: {
130:   PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));
131:   return 0;
132: }

134: static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc)
135: {
136:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

138:   *pc = fetidp->innerbddc;
139:   return 0;
140: }

142: /*@
143:  KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers

145:    Input Parameters:
146: +  ksp - the FETI-DP Krylov solver
147: -  pc - the BDDC preconditioner

149:    Level: advanced

151:    Notes:

153: .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP
154: @*/
155: PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc)
156: {
159:   PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));
160:   return 0;
161: }

163: static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
164: {
165:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

167:   PetscObjectReference((PetscObject)pc);
168:   PCDestroy(&fetidp->innerbddc);
169:   fetidp->innerbddc = pc;
170:   fetidp->userbddc  = PETSC_TRUE;
171:   return 0;
172: }

174: /*@
175:  KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers

177:    Collective on ksp

179:    Input Parameters:
180: +  ksp - the FETI-DP Krylov solver
181: -  pc - the BDDC preconditioner

183:    Level: advanced

185:    Notes:

187: .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP
188: @*/
189: PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
190: {
191:   PetscBool      isbddc;

195:   PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);
197:   PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));
198:   return 0;
199: }

201: static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V)
202: {
203:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
204:   Mat            F;
205:   Vec            Xl;

207:   KSPGetOperators(fetidp->innerksp,&F,NULL);
208:   KSPBuildSolution(fetidp->innerksp,NULL,&Xl);
209:   if (v) {
210:     PCBDDCMatFETIDPGetSolution(F,Xl,v);
211:     *V   = v;
212:   } else {
213:     PCBDDCMatFETIDPGetSolution(F,Xl,*V);
214:   }
215:   return 0;
216: }

218: static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx)
219: {
220:   KSP_FETIDPMon  *monctx = (KSP_FETIDPMon*)ctx;

222:   KSPMonitor(monctx->parentksp,it,rnorm);
223:   return 0;
224: }

226: static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
227: {
228:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

230:   KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);
231:   return 0;
232: }

234: static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin)
235: {
236:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

238:   KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);
239:   return 0;
240: }

242: static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer)
243: {
244:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
245:   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
246:   PC_IS          *pcis = (PC_IS*)fetidp->innerbddc->data;
247:   Mat_IS         *matis = (Mat_IS*)fetidp->innerbddc->pmat->data;
248:   Mat            F;
249:   FETIDPMat_ctx  fetidpmat_ctx;
250:   Vec            test_vec,test_vec_p = NULL,fetidp_global;
251:   IS             dirdofs,isvert;
252:   MPI_Comm       comm = PetscObjectComm((PetscObject)ksp);
253:   PetscScalar    sval,*array;
254:   PetscReal      val,rval;
255:   const PetscInt *vertex_indices;
256:   PetscInt       i,n_vertices;
257:   PetscBool      isascii;

260:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
262:   PetscViewerASCIIPrintf(viewer,"----------FETI-DP MAT  --------------\n");
263:   PetscViewerASCIIAddTab(viewer,2);
264:   KSPGetOperators(fetidp->innerksp,&F,NULL);
265:   PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
266:   MatView(F,viewer);
267:   PetscViewerPopFormat(viewer);
268:   PetscViewerASCIISubtractTab(viewer,2);
269:   MatShellGetContext(F,&fetidpmat_ctx);
270:   PetscViewerASCIIPrintf(viewer,"----------FETI-DP TESTS--------------\n");
271:   PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");
272:   PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");
273:   if (fetidp->fully_redundant) {
274:     PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");
275:   } else {
276:     PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");
277:   }
278:   PetscViewerFlush(viewer);

280:   /* Get Vertices used to define the BDDC */
281:   PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);
282:   ISGetLocalSize(isvert,&n_vertices);
283:   ISGetIndices(isvert,&vertex_indices);

285:   /******************************************************************/
286:   /* TEST A/B: Test numbering of global fetidp dofs                 */
287:   /******************************************************************/
288:   MatCreateVecs(F,&fetidp_global,NULL);
289:   VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);
290:   VecSet(fetidp_global,1.0);
291:   VecSet(test_vec,1.);
292:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
293:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
294:   if (fetidpmat_ctx->l2g_p) {
295:     VecDuplicate(fetidpmat_ctx->vP,&test_vec_p);
296:     VecSet(test_vec_p,1.);
297:     VecScatterBegin(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);
298:     VecScatterEnd(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);
299:   }
300:   VecAXPY(test_vec,-1.0,fetidpmat_ctx->lambda_local);
301:   VecNorm(test_vec,NORM_INFINITY,&val);
302:   VecDestroy(&test_vec);
303:   MPI_Reduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,0,comm);
304:   PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc: % 1.14e\n",rval);

306:   if (fetidpmat_ctx->l2g_p) {
307:     VecAXPY(test_vec_p,-1.0,fetidpmat_ctx->vP);
308:     VecNorm(test_vec_p,NORM_INFINITY,&val);
309:     MPI_Reduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,0,comm);
310:     PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc (p): % 1.14e\n",rval);
311:   }

313:   if (fetidp->fully_redundant) {
314:     VecSet(fetidp_global,0.0);
315:     VecSet(fetidpmat_ctx->lambda_local,0.5);
316:     VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
317:     VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
318:     VecSum(fetidp_global,&sval);
319:     val  = PetscRealPart(sval)-fetidpmat_ctx->n_lambda;
320:     MPI_Reduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,0,comm);
321:     PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob: % 1.14e\n",rval);
322:   }

324:   if (fetidpmat_ctx->l2g_p) {
325:     VecSet(pcis->vec1_N,1.0);
326:     VecSet(pcis->vec1_global,0.0);
327:     VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
328:     VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);

330:     VecSet(fetidp_global,0.0);
331:     VecSet(fetidpmat_ctx->vP,-1.0);
332:     VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
333:     VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
334:     VecScatterBegin(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
335:     VecScatterEnd(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
336:     VecScatterBegin(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);
337:     VecScatterEnd(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);
338:     VecSum(fetidp_global,&sval);
339:     val  = PetscRealPart(sval);
340:     MPI_Reduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,0,comm);
341:     PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob (p): % 1.14e\n",rval);
342:   }

344:   /******************************************************************/
345:   /* TEST C: It should hold B_delta*w=0, w\in\widehat{W}            */
346:   /* This is the meaning of the B matrix                            */
347:   /******************************************************************/

349:   VecSetRandom(pcis->vec1_N,NULL);
350:   VecSet(pcis->vec1_global,0.0);
351:   VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
352:   VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
353:   VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
354:   VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
355:   VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
356:   VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
357:   /* Action of B_delta */
358:   MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);
359:   VecSet(fetidp_global,0.0);
360:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
361:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
362:   VecNorm(fetidp_global,NORM_INFINITY,&val);
363:   PetscViewerASCIIPrintf(viewer,"C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",val);

365:   /******************************************************************/
366:   /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W}       */
367:   /* E_D = R_D^TR                                                   */
368:   /* P_D = B_{D,delta}^T B_{delta}                                  */
369:   /* eq.44 Mandel Tezaur and Dohrmann 2005                          */
370:   /******************************************************************/

372:   /* compute a random vector in \widetilde{W} */
373:   VecSetRandom(pcis->vec1_N,NULL);
374:   /* set zero at vertices and essential dofs */
375:   VecGetArray(pcis->vec1_N,&array);
376:   for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
377:   PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirdofs);
378:   if (dirdofs) {
379:     const PetscInt *idxs;
380:     PetscInt       ndir;

382:     ISGetLocalSize(dirdofs,&ndir);
383:     ISGetIndices(dirdofs,&idxs);
384:     for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
385:     ISRestoreIndices(dirdofs,&idxs);
386:   }
387:   VecRestoreArray(pcis->vec1_N,&array);
388:   /* store w for final comparison */
389:   VecDuplicate(pcis->vec1_B,&test_vec);
390:   VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);
391:   VecScatterEnd(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);

393:   /* Jump operator P_D : results stored in pcis->vec1_B */
394:   /* Action of B_delta */
395:   MatMult(fetidpmat_ctx->B_delta,test_vec,fetidpmat_ctx->lambda_local);
396:   VecSet(fetidp_global,0.0);
397:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
398:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
399:   /* Action of B_Ddelta^T */
400:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
401:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
402:   MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);

404:   /* Average operator E_D : results stored in pcis->vec2_B */
405:   PCBDDCScalingExtension(fetidpmat_ctx->pc,test_vec,pcis->vec1_global);
406:   VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);
407:   VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);

409:   /* test E_D=I-P_D */
410:   VecAXPY(pcis->vec1_B,1.0,pcis->vec2_B);
411:   VecAXPY(pcis->vec1_B,-1.0,test_vec);
412:   VecNorm(pcis->vec1_B,NORM_INFINITY,&val);
413:   VecDestroy(&test_vec);
414:   MPI_Reduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,0,comm);
415:   PetscViewerASCIIPrintf(viewer,"D: CHECK infty norm of E_D + P_D - I: % 1.14e\n",PetscGlobalRank,val);

417:   /******************************************************************/
418:   /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W}           */
419:   /* eq.48 Mandel Tezaur and Dohrmann 2005                          */
420:   /******************************************************************/

422:   VecSetRandom(pcis->vec1_N,NULL);
423:   /* set zero at vertices and essential dofs */
424:   VecGetArray(pcis->vec1_N,&array);
425:   for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
426:   if (dirdofs) {
427:     const PetscInt *idxs;
428:     PetscInt       ndir;

430:     ISGetLocalSize(dirdofs,&ndir);
431:     ISGetIndices(dirdofs,&idxs);
432:     for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
433:     ISRestoreIndices(dirdofs,&idxs);
434:   }
435:   VecRestoreArray(pcis->vec1_N,&array);

437:   /* Jump operator P_D : results stored in pcis->vec1_B */

439:   VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
440:   VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
441:   /* Action of B_delta */
442:   MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);
443:   VecSet(fetidp_global,0.0);
444:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
445:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);
446:   /* Action of B_Ddelta^T */
447:   VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
448:   VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
449:   MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);
450:   /* scaling */
451:   PCBDDCScalingExtension(fetidpmat_ctx->pc,pcis->vec1_B,pcis->vec1_global);
452:   VecNorm(pcis->vec1_global,NORM_INFINITY,&val);
453:   PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of R^T_D P_D: % 1.14e\n",val);

455:   if (!fetidp->fully_redundant) {
456:     /******************************************************************/
457:     /* TEST F: It should holds B_{delta}B^T_{D,delta}=I               */
458:     /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005               */
459:     /******************************************************************/
460:     VecDuplicate(fetidp_global,&test_vec);
461:     VecSetRandom(fetidp_global,NULL);
462:     if (fetidpmat_ctx->l2g_p) {
463:       VecSet(fetidpmat_ctx->vP,0.);
464:       VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);
465:       VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);
466:     }
467:     /* Action of B_Ddelta^T */
468:     VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
469:     VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
470:     MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);
471:     /* Action of B_delta */
472:     MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);
473:     VecSet(test_vec,0.0);
474:     VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);
475:     VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);
476:     VecAXPY(fetidp_global,-1.,test_vec);
477:     VecNorm(fetidp_global,NORM_INFINITY,&val);
478:     PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of P^T_D - I: % 1.14e\n",val);
479:     VecDestroy(&test_vec);
480:   }
481:   PetscViewerASCIIPrintf(viewer,"-------------------------------------\n");
482:   PetscViewerFlush(viewer);
483:   VecDestroy(&test_vec_p);
484:   ISDestroy(&dirdofs);
485:   VecDestroy(&fetidp_global);
486:   ISRestoreIndices(isvert,&vertex_indices);
487:   PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);
488:   return 0;
489: }

491: static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp)
492: {
493:   KSP_FETIDP       *fetidp = (KSP_FETIDP*)ksp->data;
494:   PC_BDDC          *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
495:   Mat              A,Ap;
496:   PetscInt         fid = -1;
497:   PetscMPIInt      size;
498:   PetscBool        ismatis,pisz,allp,schp;
499:   PetscBool        flip; /* Usually, Stokes is written (B = -\int_\Omega \nabla \cdot u q)
500:                            | A B'| | v | = | f |
501:                            | B 0 | | p | = | g |
502:                             If -ksp_fetidp_saddlepoint_flip is true, the code assumes it is written as
503:                            | A B'| | v | = | f |
504:                            |-B 0 | | p | = |-g |
505:                          */
506:   PetscObjectState matstate, matnnzstate;
507:   PetscErrorCode   ierr;

509:   pisz = PETSC_FALSE;
510:   flip = PETSC_FALSE;
511:   allp = PETSC_FALSE;
512:   schp = PETSC_FALSE;
513:   PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");
514:   PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);
515:   PetscOptionsBool("-ksp_fetidp_pressure_all","Use the whole pressure set instead of just that at the interface",NULL,allp,&allp,NULL);
516:   PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);
517:   PetscOptionsBool("-ksp_fetidp_pressure_schur","Use a BDDC solver for pressure",NULL,schp,&schp,NULL);
518:   PetscOptionsEnd();

520:   MPI_Comm_size(PetscObjectComm((PetscObject)ksp),&size);
521:   fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint);
522:   if (size == 1) fetidp->saddlepoint = PETSC_FALSE;

524:   KSPGetOperators(ksp,&A,&Ap);
525:   PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);

528:   /* Quiet return if the matrix states are unchanged.
529:      Needed only for the saddle point case since it uses MatZeroRows
530:      on a matrix that may not have changed */
531:   PetscObjectStateGet((PetscObject)A,&matstate);
532:   MatGetNonzeroState(A,&matnnzstate);
533:   if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) return 0;
534:   fetidp->matstate     = matstate;
535:   fetidp->matnnzstate  = matnnzstate;
536:   fetidp->statechanged = fetidp->saddlepoint;

538:   /* see if we have some fields attached */
539:   if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
540:     DM             dm;
541:     PetscContainer c;

543:     KSPGetDM(ksp,&dm);
544:     PetscObjectQuery((PetscObject)A,"_convert_nest_lfields",(PetscObject*)&c);
545:     if (dm) {
546:       IS      *fields;
547:       PetscInt nf,i;

549:       DMCreateFieldDecomposition(dm,&nf,NULL,&fields,NULL);
550:       PCBDDCSetDofsSplitting(fetidp->innerbddc,nf,fields);
551:       for (i=0;i<nf;i++) {
552:         ISDestroy(&fields[i]);
553:       }
554:       PetscFree(fields);
555:     } else if (c) {
556:       MatISLocalFields lf;

558:       PetscContainerGetPointer(c,(void**)&lf);
559:       PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);
560:     }
561:   }

563:   if (!fetidp->saddlepoint) {
564:     PCSetOperators(fetidp->innerbddc,A,A);
565:   } else {
566:     Mat          nA,lA,PPmat;
567:     MatNullSpace nnsp;
568:     IS           pP;
569:     PetscInt     totP;

571:     MatISGetLocalMat(A,&lA);
572:     PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);

574:     pP = fetidp->pP;
575:     if (!pP) { /* first time, need to compute pressure dofs */
576:       PC_IS                  *pcis = (PC_IS*)fetidp->innerbddc->data;
577:       Mat_IS                 *matis = (Mat_IS*)(A->data);
578:       ISLocalToGlobalMapping l2g;
579:       IS                     lP = NULL,II,pII,lPall,Pall,is1,is2;
580:       const PetscInt         *idxs;
581:       PetscInt               nl,ni,*widxs;
582:       PetscInt               i,j,n_neigh,*neigh,*n_shared,**shared,*count;
583:       PetscInt               rst,ren,n;
584:       PetscBool              ploc;

586:       MatGetLocalSize(A,&nl,NULL);
587:       MatGetOwnershipRange(A,&rst,&ren);
588:       MatGetLocalSize(lA,&n,NULL);
589:       MatISGetLocalToGlobalMapping(A,&l2g,NULL);

591:       if (!pcis->is_I_local) { /* need to compute interior dofs */
592:         PetscCalloc1(n,&count);
593:         ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);
594:         for (i=1;i<n_neigh;i++)
595:           for (j=0;j<n_shared[i];j++)
596:             count[shared[i][j]] += 1;
597:         for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i;
598:         ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);
599:         ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);
600:       } else {
601:         PetscObjectReference((PetscObject)pcis->is_I_local);
602:         II   = pcis->is_I_local;
603:       }

605:       /* interior dofs in layout */
606:       PetscArrayzero(matis->sf_leafdata,n);
607:       PetscArrayzero(matis->sf_rootdata,nl);
608:       ISGetLocalSize(II,&ni);
609:       ISGetIndices(II,&idxs);
610:       for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
611:       ISRestoreIndices(II,&idxs);
612:       PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
613:       PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
614:       PetscMalloc1(PetscMax(nl,n),&widxs);
615:       for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
616:       ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);

618:       /* pressure dofs */
619:       Pall  = NULL;
620:       lPall = NULL;
621:       ploc  = PETSC_FALSE;
622:       if (fid < 0) { /* zero pressure block */
623:         PetscInt np;

625:         MatFindZeroDiagonals(A,&Pall);
626:         ISGetSize(Pall,&np);
627:         if (!np) { /* zero-block not found, defaults to last field (if set) */
628:           fid  = pcbddc->n_ISForDofsLocal ? pcbddc->n_ISForDofsLocal - 1 : pcbddc->n_ISForDofs - 1;
629:           ISDestroy(&Pall);
630:         } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
631:           PCBDDCSetDofsSplitting(fetidp->innerbddc,1,&Pall);
632:         }
633:       }
634:       if (!Pall) { /* look for registered fields */
635:         if (pcbddc->n_ISForDofsLocal) {
636:           PetscInt np;

639:           /* need a sequential IS */
640:           ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);
641:           ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);
642:           ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);
643:           ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);
644:           ploc = PETSC_TRUE;
645:         } else if (pcbddc->n_ISForDofs) {
647:           PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);
648:           Pall = pcbddc->ISForDofs[fid];
649:         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal");
650:       }

652:       /* if the user requested the entire pressure,
653:          remove the interior pressure dofs from II (or pII) */
654:       if (allp) {
655:         if (ploc) {
656:           IS nII;
657:           ISDifference(II,lPall,&nII);
658:           ISDestroy(&II);
659:           II   = nII;
660:         } else {
661:           IS nII;
662:           ISDifference(pII,Pall,&nII);
663:           ISDestroy(&pII);
664:           pII  = nII;
665:         }
666:       }
667:       if (ploc) {
668:         ISDifference(lPall,II,&lP);
669:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);
670:       } else {
671:         ISDifference(Pall,pII,&pP);
672:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);
673:         /* need all local pressure dofs */
674:         PetscArrayzero(matis->sf_leafdata,n);
675:         PetscArrayzero(matis->sf_rootdata,nl);
676:         ISGetLocalSize(Pall,&ni);
677:         ISGetIndices(Pall,&idxs);
678:         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
679:         ISRestoreIndices(Pall,&idxs);
680:         PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
681:         PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
682:         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
683:         ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);
684:       }

686:       if (!Pall) {
687:         PetscArrayzero(matis->sf_leafdata,n);
688:         PetscArrayzero(matis->sf_rootdata,nl);
689:         ISGetLocalSize(lPall,&ni);
690:         ISGetIndices(lPall,&idxs);
691:         for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
692:         ISRestoreIndices(lPall,&idxs);
693:         PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
694:         PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
695:         for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
696:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);
697:       }
698:       PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject)Pall);

700:       if (flip) {
701:         PetscInt npl;
702:         ISGetLocalSize(Pall,&npl);
703:         ISGetIndices(Pall,&idxs);
704:         MatCreateVecs(A,NULL,&fetidp->rhs_flip);
705:         VecSet(fetidp->rhs_flip,1.);
706:         VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
707:         for (i=0;i<npl;i++) {
708:           VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);
709:         }
710:         VecAssemblyBegin(fetidp->rhs_flip);
711:         VecAssemblyEnd(fetidp->rhs_flip);
712:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);
713:         ISRestoreIndices(Pall,&idxs);
714:       }
715:       ISDestroy(&Pall);
716:       ISDestroy(&pII);

718:       /* local selected pressures in subdomain-wise and global ordering */
719:       PetscArrayzero(matis->sf_leafdata,n);
720:       PetscArrayzero(matis->sf_rootdata,nl);
721:       if (!ploc) {
722:         PetscInt *widxs2;

725:         ISGetLocalSize(pP,&ni);
726:         ISGetIndices(pP,&idxs);
727:         for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
728:         ISRestoreIndices(pP,&idxs);
729:         PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
730:         PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE);
731:         for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
732:         PetscMalloc1(ni,&widxs2);
733:         ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs2);
734:         ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);
735:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);
736:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs2,PETSC_OWN_POINTER,&is1);
737:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);
738:         ISDestroy(&is1);
739:       } else {
741:         ISGetLocalSize(lP,&ni);
742:         ISGetIndices(lP,&idxs);
743:         for (i=0;i<ni;i++)
744:           if (idxs[i] >=0 && idxs[i] < n)
745:             matis->sf_leafdata[idxs[i]] = 1;
746:         ISRestoreIndices(lP,&idxs);
747:         PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
748:         ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);
749:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);
750:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);
751:         ISDestroy(&is1);
752:         PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPI_REPLACE);
753:         for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
754:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);
755:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);
756:       }
757:       PetscFree(widxs);

759:       /* If there's any "interior pressure",
760:          we may want to use a discrete harmonic solver instead
761:          of a Stokes harmonic for the Dirichlet preconditioner
762:          Need to extract the interior velocity dofs in interior dofs ordering (iV)
763:          and interior pressure dofs in local ordering (iP) */
764:       if (!allp) {
765:         ISLocalToGlobalMapping l2g_t;

767:         ISDifference(lPall,lP,&is1);
768:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);
769:         ISDifference(II,is1,&is2);
770:         ISDestroy(&is1);
771:         ISLocalToGlobalMappingCreateIS(II,&l2g_t);
772:         ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);
773:         ISGetLocalSize(is1,&i);
774:         ISGetLocalSize(is2,&j);
776:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);
777:         ISLocalToGlobalMappingDestroy(&l2g_t);
778:         ISDestroy(&is1);
779:         ISDestroy(&is2);
780:       }
781:       ISDestroy(&II);

783:       /* exclude selected pressures from the inner BDDC */
784:       if (pcbddc->DirichletBoundariesLocal) {
785:         IS       list[2],plP,isout;
786:         PetscInt np;

788:         /* need a parallel IS */
789:         ISGetLocalSize(lP,&np);
790:         ISGetIndices(lP,&idxs);
791:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);
792:         list[0] = plP;
793:         list[1] = pcbddc->DirichletBoundariesLocal;
794:         ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);
795:         ISSortRemoveDups(isout);
796:         ISDestroy(&plP);
797:         ISRestoreIndices(lP,&idxs);
798:         PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);
799:         ISDestroy(&isout);
800:       } else if (pcbddc->DirichletBoundaries) {
801:         IS list[2],isout;

803:         list[0] = pP;
804:         list[1] = pcbddc->DirichletBoundaries;
805:         ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);
806:         ISSortRemoveDups(isout);
807:         PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);
808:         ISDestroy(&isout);
809:       } else {
810:         IS       plP;
811:         PetscInt np;

813:         /* need a parallel IS */
814:         ISGetLocalSize(lP,&np);
815:         ISGetIndices(lP,&idxs);
816:         ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);
817:         PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);
818:         ISDestroy(&plP);
819:         ISRestoreIndices(lP,&idxs);
820:       }

822:       /* save CSR information for the pressure BDDC solver (if any) */
823:       if (schp) {
824:         PetscInt np,nt;

826:         MatGetSize(matis->A,&nt,NULL);
827:         ISGetLocalSize(lP,&np);
828:         if (np) {
829:           PetscInt *xadj = pcbddc->mat_graph->xadj;
830:           PetscInt *adjn = pcbddc->mat_graph->adjncy;
831:           PetscInt nv = pcbddc->mat_graph->nvtxs_csr;

833:           if (nv && nv == nt) {
834:             ISLocalToGlobalMapping pmap;
835:             PetscInt               *schp_csr,*schp_xadj,*schp_adjn,p;
836:             PetscContainer         c;

838:             ISLocalToGlobalMappingCreateIS(lPall,&pmap);
839:             ISGetIndices(lPall,&idxs);
840:             for (p = 0, nv = 0; p < np; p++) {
841:               PetscInt x,n = idxs[p];

843:               ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,NULL);
844:               nv  += x;
845:             }
846:             PetscMalloc1(np + 1 + nv,&schp_csr);
847:             schp_xadj = schp_csr;
848:             schp_adjn = schp_csr + np + 1;
849:             for (p = 0, schp_xadj[0] = 0; p < np; p++) {
850:               PetscInt x,n = idxs[p];

852:               ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,schp_adjn + schp_xadj[p]);
853:               schp_xadj[p+1] = schp_xadj[p] + x;
854:             }
855:             ISRestoreIndices(lPall,&idxs);
856:             ISLocalToGlobalMappingDestroy(&pmap);
857:             PetscContainerCreate(PETSC_COMM_SELF,&c);
858:             PetscContainerSetPointer(c,schp_csr);
859:             PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);
860:             PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pCSR",(PetscObject)c);
861:             PetscContainerDestroy(&c);

863:           }
864:         }
865:       }
866:       ISDestroy(&lPall);
867:       ISDestroy(&lP);
868:       fetidp->pP = pP;
869:     }

871:     /* total number of selected pressure dofs */
872:     ISGetSize(fetidp->pP,&totP);

874:     /* Set operator for inner BDDC */
875:     if (totP || fetidp->rhs_flip) {
876:       MatDuplicate(A,MAT_COPY_VALUES,&nA);
877:     } else {
878:       PetscObjectReference((PetscObject)A);
879:       nA   = A;
880:     }
881:     if (fetidp->rhs_flip) {
882:       MatDiagonalScale(nA,fetidp->rhs_flip,NULL);
883:       if (totP) {
884:         Mat lA2;

886:         MatISGetLocalMat(nA,&lA);
887:         MatDuplicate(lA,MAT_COPY_VALUES,&lA2);
888:         PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);
889:         MatDestroy(&lA2);
890:       }
891:     }

893:     if (totP) {
894:       MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
895:       MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);
896:     } else {
897:       PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);
898:     }
899:     MatGetNearNullSpace(Ap,&nnsp);
900:     if (!nnsp) {
901:       MatGetNullSpace(Ap,&nnsp);
902:     }
903:     if (!nnsp) {
904:       MatGetNearNullSpace(A,&nnsp);
905:     }
906:     if (!nnsp) {
907:       MatGetNullSpace(A,&nnsp);
908:     }
909:     MatSetNearNullSpace(nA,nnsp);
910:     PCSetOperators(fetidp->innerbddc,nA,nA);
911:     MatDestroy(&nA);

913:     /* non-zero rhs on interior dofs when applying the preconditioner */
914:     if (totP) pcbddc->switch_static = PETSC_TRUE;

916:     /* if there are no interface pressures, set inner bddc flag for benign saddle point */
917:     if (!totP) {
918:       pcbddc->benign_saddle_point = PETSC_TRUE;
919:       pcbddc->compute_nonetflux   = PETSC_TRUE;
920:     }

922:     /* Operators for pressure preconditioner */
923:     if (totP) {
924:       /* Extract pressure block if needed */
925:       if (!pisz) {
926:         Mat C;
927:         IS  nzrows = NULL;

929:         MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);
930:         MatFindNonzeroRows(C,&nzrows);
931:         if (nzrows) {
932:           PetscInt i;

934:           ISGetSize(nzrows,&i);
935:           ISDestroy(&nzrows);
936:           if (!i) pisz = PETSC_TRUE;
937:         }
938:         if (!pisz) {
939:           MatScale(C,-1.); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized */
940:           PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);
941:         }
942:         MatDestroy(&C);
943:       }
944:       /* Divergence mat */
945:       if (!pcbddc->divudotp) {
946:         Mat       B;
947:         IS        P;
948:         IS        l2l = NULL;
949:         PetscBool save;

951:         PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);
952:         if (!pisz) {
953:           IS       F,V;
954:           PetscInt m,M;

956:           MatGetOwnershipRange(A,&m,&M);
957:           ISCreateStride(PetscObjectComm((PetscObject)A),M-m,m,1,&F);
958:           ISComplement(P,m,M,&V);
959:           MatCreateSubMatrix(A,P,V,MAT_INITIAL_MATRIX,&B);
960:           {
961:             Mat_IS *Bmatis = (Mat_IS*)B->data;
962:             PetscObjectReference((PetscObject)Bmatis->getsub_cis);
963:             l2l  = Bmatis->getsub_cis;
964:           }
965:           ISDestroy(&V);
966:           ISDestroy(&F);
967:         } else {
968:           MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);
969:         }
970:         save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */
971:         PCBDDCSetDivergenceMat(fetidp->innerbddc,B,PETSC_FALSE,l2l);
972:         pcbddc->compute_nonetflux = save;
973:         MatDestroy(&B);
974:         ISDestroy(&l2l);
975:       }
976:       if (A != Ap) { /* user has provided a different Pmat, this always superseeds the setter (TODO: is it OK?) */
977:         /* use monolithic operator, we restrict later */
978:         KSPFETIDPSetPressureOperator(ksp,Ap);
979:       }
980:       PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);

982:       /* PPmat not present, use some default choice */
983:       if (!PPmat) {
984:         Mat C;

986:         PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);
987:         if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */
988:           KSPFETIDPSetPressureOperator(ksp,C);
989:         } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */
990:           IS  P;

992:           PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);
993:           MatCreateSubMatrix(A,P,P,MAT_INITIAL_MATRIX,&C);
994:           MatScale(C,-1.);
995:           KSPFETIDPSetPressureOperator(ksp,C);
996:           MatDestroy(&C);
997:         } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */
998:           PetscInt nl;

1000:           ISGetLocalSize(fetidp->pP,&nl);
1001:           MatCreate(PetscObjectComm((PetscObject)ksp),&C);
1002:           MatSetSizes(C,nl,nl,totP,totP);
1003:           MatSetType(C,MATAIJ);
1004:           MatMPIAIJSetPreallocation(C,1,NULL,0,NULL);
1005:           MatSeqAIJSetPreallocation(C,1,NULL);
1006:           MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1007:           MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1008:           MatShift(C,1.);
1009:           KSPFETIDPSetPressureOperator(ksp,C);
1010:           MatDestroy(&C);
1011:         }
1012:       }

1014:       /* Preconditioned operator for the pressure block */
1015:       PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);
1016:       if (PPmat) {
1017:         Mat      C;
1018:         IS       Pall;
1019:         PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg;

1021:         PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);
1022:         MatGetSize(A,&AM,NULL);
1023:         MatGetSize(PPmat,&PAM,&PAN);
1024:         ISGetSize(Pall,&pAg);
1025:         ISGetSize(fetidp->pP,&pIg);
1026:         MatGetLocalSize(PPmat,&pam,&pan);
1027:         MatGetLocalSize(A,&am,&an);
1028:         ISGetLocalSize(Pall,&pIl);
1029:         ISGetLocalSize(fetidp->pP,&pl);
1034:         if (PAM == AM) { /* monolithic ordering, restrict to pressure */
1035:           if (schp) {
1036:             MatCreateSubMatrix(PPmat,Pall,Pall,MAT_INITIAL_MATRIX,&C);
1037:           } else {
1038:             MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);
1039:           }
1040:         } else if (pAg == PAM) { /* global ordering for pressure only */
1041:           if (!allp && !schp) { /* solving for interface pressure only */
1042:             IS restr;

1044:             ISRenumber(fetidp->pP,NULL,NULL,&restr);
1045:             MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);
1046:             ISDestroy(&restr);
1047:           } else {
1048:             PetscObjectReference((PetscObject)PPmat);
1049:             C    = PPmat;
1050:           }
1051:         } else if (pIg == PAM) { /* global ordering for selected pressure only */
1053:           PetscObjectReference((PetscObject)PPmat);
1054:           C    = PPmat;
1055:         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix");

1057:         KSPFETIDPSetPressureOperator(ksp,C);
1058:         MatDestroy(&C);
1059:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing Pmat for pressure block");
1060:     } else { /* totP == 0 */
1061:       PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);
1062:     }
1063:   }
1064:   return 0;
1065: }

1067: static PetscErrorCode KSPSetUp_FETIDP(KSP ksp)
1068: {
1069:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1070:   PC_BDDC        *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1071:   PetscBool      flg;

1073:   KSPFETIDPSetUpOperators(ksp);
1074:   /* set up BDDC */
1075:   PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);
1076:   PCSetUp(fetidp->innerbddc);
1077:   /* FETI-DP as it is implemented needs an exact coarse solver */
1078:   if (pcbddc->coarse_ksp) {
1079:     KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);
1080:     KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);
1081:   }
1082:   /* FETI-DP as it is implemented needs exact local Neumann solvers */
1083:   KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);
1084:   KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);

1086:   /* setup FETI-DP operators
1087:      If fetidp->statechanged is true, we need to update the operators
1088:      needed in the saddle-point case. This should be replaced
1089:      by a better logic when the FETI-DP matrix and preconditioner will
1090:      have their own classes */
1091:   if (pcbddc->new_primal_space || fetidp->statechanged) {
1092:     Mat F; /* the FETI-DP matrix */
1093:     PC  D; /* the FETI-DP preconditioner */
1094:     KSPReset(fetidp->innerksp);
1095:     PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);
1096:     KSPSetOperators(fetidp->innerksp,F,F);
1097:     KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);
1098:     KSPSetPC(fetidp->innerksp,D);
1099:     PetscObjectIncrementTabLevel((PetscObject)D,(PetscObject)fetidp->innerksp,0);
1100:     KSPSetFromOptions(fetidp->innerksp);
1101:     MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);
1102:     MatDestroy(&F);
1103:     PCDestroy(&D);
1104:     if (fetidp->check) {
1105:       PetscViewer viewer;

1107:       if (!pcbddc->dbg_viewer) {
1108:         viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
1109:       } else {
1110:         viewer = pcbddc->dbg_viewer;
1111:       }
1112:       KSPFETIDPCheckOperators(ksp,viewer);
1113:     }
1114:   }
1115:   fetidp->statechanged     = PETSC_FALSE;
1116:   pcbddc->new_primal_space = PETSC_FALSE;

1118:   /* propagate settings to the inner solve */
1119:   KSPGetComputeSingularValues(ksp,&flg);
1120:   KSPSetComputeSingularValues(fetidp->innerksp,flg);
1121:   if (ksp->res_hist) {
1122:     KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);
1123:   }
1124:   KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);
1125:   KSPSetUp(fetidp->innerksp);
1126:   return 0;
1127: }

1129: static PetscErrorCode KSPSolve_FETIDP(KSP ksp)
1130: {
1131:   Mat                F,A;
1132:   MatNullSpace       nsp;
1133:   Vec                X,B,Xl,Bl;
1134:   KSP_FETIDP         *fetidp = (KSP_FETIDP*)ksp->data;
1135:   PC_BDDC            *pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1136:   KSPConvergedReason reason;
1137:   PC                 pc;
1138:   PCFailedReason     pcreason;
1139:   PetscInt           hist_len;

1141:   PetscCitationsRegister(citation,&cited);
1142:   if (fetidp->saddlepoint) {
1143:     PetscCitationsRegister(citation2,&cited2);
1144:   }
1145:   KSPGetOperators(ksp,&A,NULL);
1146:   KSPGetRhs(ksp,&B);
1147:   KSPGetSolution(ksp,&X);
1148:   KSPGetOperators(fetidp->innerksp,&F,NULL);
1149:   KSPGetRhs(fetidp->innerksp,&Bl);
1150:   KSPGetSolution(fetidp->innerksp,&Xl);
1151:   PCBDDCMatFETIDPGetRHS(F,B,Bl);
1152:   if (ksp->transpose_solve) {
1153:     KSPSolveTranspose(fetidp->innerksp,Bl,Xl);
1154:   } else {
1155:     KSPSolve(fetidp->innerksp,Bl,Xl);
1156:   }
1157:   KSPGetConvergedReason(fetidp->innerksp,&reason);
1158:   KSPGetPC(fetidp->innerksp,&pc);
1159:   PCGetFailedReason(pc,&pcreason);
1160:   if ((reason < 0 && reason != KSP_DIVERGED_ITS) || pcreason) {
1161:     PetscInt its;
1162:     KSPGetIterationNumber(fetidp->innerksp,&its);
1163:     ksp->reason = KSP_DIVERGED_PC_FAILED;
1164:     VecSetInf(Xl);
1165:     PetscInfo(ksp,"Inner KSP solve failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
1166:   }
1167:   PCBDDCMatFETIDPGetSolution(F,Xl,X);
1168:   MatGetNullSpace(A,&nsp);
1169:   if (nsp) {
1170:     MatNullSpaceRemove(nsp,X);
1171:   }
1172:   /* update ksp with stats from inner ksp */
1173:   KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);
1174:   KSPGetIterationNumber(fetidp->innerksp,&ksp->its);
1175:   ksp->totalits += ksp->its;
1176:   KSPGetResidualHistory(fetidp->innerksp,NULL,&hist_len);
1177:   ksp->res_hist_len = (size_t) hist_len;
1178:   /* restore defaults for inner BDDC (Pre/PostSolve flags) */
1179:   pcbddc->temp_solution_used        = PETSC_FALSE;
1180:   pcbddc->rhs_change                = PETSC_FALSE;
1181:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1182:   return 0;
1183: }

1185: static PetscErrorCode KSPReset_FETIDP(KSP ksp)
1186: {
1187:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1188:   PC_BDDC        *pcbddc;

1190:   ISDestroy(&fetidp->pP);
1191:   VecDestroy(&fetidp->rhs_flip);
1192:   /* avoid PCReset that does not take into account ref counting */
1193:   PCDestroy(&fetidp->innerbddc);
1194:   PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);
1195:   PCSetType(fetidp->innerbddc,PCBDDC);
1196:   pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1197:   pcbddc->symmetric_primal = PETSC_FALSE;
1198:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);
1199:   KSPDestroy(&fetidp->innerksp);
1200:   fetidp->saddlepoint  = PETSC_FALSE;
1201:   fetidp->matstate     = -1;
1202:   fetidp->matnnzstate  = -1;
1203:   fetidp->statechanged = PETSC_TRUE;
1204:   return 0;
1205: }

1207: static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
1208: {
1209:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

1211:   KSPReset_FETIDP(ksp);
1212:   PCDestroy(&fetidp->innerbddc);
1213:   KSPDestroy(&fetidp->innerksp);
1214:   PetscFree(fetidp->monctx);
1215:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);
1216:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);
1217:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);
1218:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);
1219:   PetscFree(ksp->data);
1220:   return 0;
1221: }

1223: static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer)
1224: {
1225:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;
1226:   PetscBool      iascii;

1228:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1229:   if (iascii) {
1230:     PetscViewerASCIIPrintf(viewer,"  fully redundant: %d\n",fetidp->fully_redundant);
1231:     PetscViewerASCIIPrintf(viewer,"  saddle point:    %d\n",fetidp->saddlepoint);
1232:     PetscViewerASCIIPrintf(viewer,"Inner KSP solver details\n");
1233:   }
1234:   PetscViewerASCIIPushTab(viewer);
1235:   KSPView(fetidp->innerksp,viewer);
1236:   PetscViewerASCIIPopTab(viewer);
1237:   if (iascii) {
1238:     PetscViewerASCIIPrintf(viewer,"Inner BDDC solver details\n");
1239:   }
1240:   PetscViewerASCIIPushTab(viewer);
1241:   PCView(fetidp->innerbddc,viewer);
1242:   PetscViewerASCIIPopTab(viewer);
1243:   return 0;
1244: }

1246: static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp)
1247: {
1248:   KSP_FETIDP     *fetidp = (KSP_FETIDP*)ksp->data;

1250:   /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
1251:   PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);
1252:   PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");
1253:   if (!fetidp->userbddc) {
1254:     PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);
1255:     PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");
1256:   }
1257:   PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");
1258:   PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);
1259:   PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);
1260:   PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);
1261:   PetscOptionsTail();
1262:   PCSetFromOptions(fetidp->innerbddc);
1263:   return 0;
1264: }

1266: /*MC
1267:      KSPFETIDP - The FETI-DP method

1269:    This class implements the FETI-DP method [1].
1270:    The matrix for the KSP must be of type MATIS.
1271:    The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object.

1273:    Options Database Keys:
1274: +   -ksp_fetidp_fullyredundant <false>   - use a fully redundant set of Lagrange multipliers
1275: .   -ksp_fetidp_saddlepoint <false>      - activates support for saddle point problems, see [2]
1276: .   -ksp_fetidp_saddlepoint_flip <false> - usually, an incompressible Stokes problem is written as
1277:                                            | A B^T | | v | = | f |
1278:                                            | B 0   | | p | = | g |
1279:                                            with B representing -\int_\Omega \nabla \cdot u q.
1280:                                            If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as
1281:                                            | A B^T | | v | = | f |
1282:                                            |-B 0   | | p | = |-g |
1283: .   -ksp_fetidp_pressure_field <-1>      - activates support for saddle point problems, and identifies the pressure field id.
1284:                                            If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals().
1285: -   -ksp_fetidp_pressure_all <false>     - if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field.

1287:    Level: Advanced

1289:    Notes:
1290:     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.,
1291: .vb
1292:       -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false
1293: .ve
1294:    will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC.

1296:    For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator().
1297:    Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix).
1298:    If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver.
1299:    Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g.
1300: .vb
1301:       -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps
1302: .ve
1303:    In order to use the deluxe version of FETI-DP, you must customize the inner BDDC operator with -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_deluxe_singlemat and use
1304:    non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g.
1305: .vb
1306:       -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu
1307: .ve

1309:    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.

1311:    The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution.

1313:    Developer Notes:
1314:     Even though this method does not directly use any norms, the user is allowed to set the KSPNormType to any value.
1315:     This is so users do not have to change KSPNormType options when they switch from other KSP methods to this one.

1317:    References:
1318: +  * - 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
1319: -  * - 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

1321: .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP()
1322: M*/
1323: PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
1324: {
1325:   KSP_FETIDP     *fetidp;
1326:   KSP_FETIDPMon  *monctx;
1327:   PC_BDDC        *pcbddc;
1328:   PC             pc;

1330:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);
1331:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);
1332:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1333:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);
1334:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
1335:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
1336:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);

1338:   PetscNewLog(ksp,&fetidp);
1339:   fetidp->matstate     = -1;
1340:   fetidp->matnnzstate  = -1;
1341:   fetidp->statechanged = PETSC_TRUE;

1343:   ksp->data = (void*)fetidp;
1344:   ksp->ops->setup                        = KSPSetUp_FETIDP;
1345:   ksp->ops->solve                        = KSPSolve_FETIDP;
1346:   ksp->ops->destroy                      = KSPDestroy_FETIDP;
1347:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_FETIDP;
1348:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP;
1349:   ksp->ops->view                         = KSPView_FETIDP;
1350:   ksp->ops->setfromoptions               = KSPSetFromOptions_FETIDP;
1351:   ksp->ops->buildsolution                = KSPBuildSolution_FETIDP;
1352:   ksp->ops->buildresidual                = KSPBuildResidualDefault;
1353:   KSPGetPC(ksp,&pc);
1354:   PCSetType(pc,PCNONE);
1355:   /* create the inner KSP for the Lagrange multipliers */
1356:   KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);
1357:   KSPGetPC(fetidp->innerksp,&pc);
1358:   PCSetType(pc,PCNONE);
1359:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);
1360:   /* monitor */
1361:   PetscNew(&monctx);
1362:   monctx->parentksp = ksp;
1363:   fetidp->monctx = monctx;
1364:   KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);
1365:   /* create the inner BDDC */
1366:   PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);
1367:   PCSetType(fetidp->innerbddc,PCBDDC);
1368:   /* make sure we always obtain a consistent FETI-DP matrix
1369:      for symmetric problems, the user can always customize it through the command line */
1370:   pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1371:   pcbddc->symmetric_primal = PETSC_FALSE;
1372:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);
1373:   /* composed functions */
1374:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);
1375:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);
1376:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);
1377:   PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);
1378:   /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
1379:   ksp->setupnewmatrix = PETSC_TRUE;
1380:   return 0;
1381: }