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