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