1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/
2 #include <../src/ksp/pc/impls/bddc/bddc.h>
3 #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
4 #include <petscdm.h>
5
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";
43
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;
50
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;
66
KSPFETIDPSetPressureOperator_FETIDP(KSP ksp,Mat P)67 static PetscErrorCode KSPFETIDPSetPressureOperator_FETIDP(KSP ksp, Mat P)
68 {
69 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
70 PetscErrorCode ierr;
71
72 PetscFunctionBegin;
73 if (P) fetidp->saddlepoint = PETSC_TRUE;
74 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)P);CHKERRQ(ierr);
75 PetscFunctionReturn(0);
76 }
77
78 /*@
79 KSPFETIDPSetPressureOperator - Sets the operator used to setup the pressure preconditioner for saddle point FETI-DP.
80
81 Collective on ksp
82
83 Input Parameters:
84 + ksp - the FETI-DP Krylov solver
85 - P - the linear operator to be preconditioned, usually the mass matrix.
86
87 Level: advanced
88
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.
96
97 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators
98 @*/
KSPFETIDPSetPressureOperator(KSP ksp,Mat P)99 PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P)
100 {
101 PetscErrorCode ierr;
102
103 PetscFunctionBegin;
104 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
105 if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,2);
106 ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperator_C",(KSP,Mat),(ksp,P));CHKERRQ(ierr);
107 PetscFunctionReturn(0);
108 }
109
KSPFETIDPGetInnerKSP_FETIDP(KSP ksp,KSP * innerksp)110 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp)
111 {
112 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
113
114 PetscFunctionBegin;
115 *innerksp = fetidp->innerksp;
116 PetscFunctionReturn(0);
117 }
118
119 /*@
120 KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers
121
122 Input Parameters:
123 + ksp - the FETI-DP KSP
124 - innerksp - the KSP for the multipliers
125
126 Level: advanced
127
128 Notes:
129
130 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC
131 @*/
KSPFETIDPGetInnerKSP(KSP ksp,KSP * innerksp)132 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp)
133 {
134 PetscErrorCode ierr;
135
136 PetscFunctionBegin;
137 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
138 PetscValidPointer(innerksp,2);
139 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr);
140 PetscFunctionReturn(0);
141 }
142
KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp,PC * pc)143 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc)
144 {
145 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
146
147 PetscFunctionBegin;
148 *pc = fetidp->innerbddc;
149 PetscFunctionReturn(0);
150 }
151
152 /*@
153 KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
154
155 Input Parameters:
156 + ksp - the FETI-DP Krylov solver
157 - pc - the BDDC preconditioner
158
159 Level: advanced
160
161 Notes:
162
163 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP
164 @*/
KSPFETIDPGetInnerBDDC(KSP ksp,PC * pc)165 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc)
166 {
167 PetscErrorCode ierr;
168
169 PetscFunctionBegin;
170 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
171 PetscValidPointer(pc,2);
172 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr);
173 PetscFunctionReturn(0);
174 }
175
KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp,PC pc)176 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
177 {
178 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
179 PetscErrorCode ierr;
180
181 PetscFunctionBegin;
182 ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr);
183 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
184 fetidp->innerbddc = pc;
185 fetidp->userbddc = PETSC_TRUE;
186 PetscFunctionReturn(0);
187 }
188
189 /*@
190 KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers
191
192 Collective on ksp
193
194 Input Parameters:
195 + ksp - the FETI-DP Krylov solver
196 - pc - the BDDC preconditioner
197
198 Level: advanced
199
200 Notes:
201
202 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP
203 @*/
KSPFETIDPSetInnerBDDC(KSP ksp,PC pc)204 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
205 {
206 PetscBool isbddc;
207 PetscErrorCode ierr;
208
209 PetscFunctionBegin;
210 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
211 PetscValidHeaderSpecific(pc,PC_CLASSID,2);
212 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr);
213 if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner");
214 ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr);
215 PetscFunctionReturn(0);
216 }
217
KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec * V)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;
223 PetscErrorCode ierr;
224
225 PetscFunctionBegin;
226 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
227 ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr);
228 if (v) {
229 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr);
230 *V = v;
231 } else {
232 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr);
233 }
234 PetscFunctionReturn(0);
235 }
236
KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void * ctx)237 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx)
238 {
239 KSP_FETIDPMon *monctx = (KSP_FETIDPMon*)ctx;
240 PetscErrorCode ierr;
241
242 PetscFunctionBegin;
243 ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr);
244 PetscFunctionReturn(0);
245 }
246
KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal * r,PetscReal * c,PetscInt * neig)247 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
248 {
249 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
250 PetscErrorCode ierr;
251
252 PetscFunctionBegin;
253 ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr);
254 PetscFunctionReturn(0);
255 }
256
KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal * emax,PetscReal * emin)257 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin)
258 {
259 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
260 PetscErrorCode ierr;
261
262 PetscFunctionBegin;
263 ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr);
264 PetscFunctionReturn(0);
265 }
266
KSPFETIDPCheckOperators(KSP ksp,PetscViewer viewer)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;
283 PetscErrorCode ierr;
284
285 PetscFunctionBegin;
286 PetscCheckSameComm(ksp,1,viewer,2);
287 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
288 if (!isascii) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported viewer");
289 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP MAT --------------\n");CHKERRQ(ierr);
290 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr);
291 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
292 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
293 ierr = MatView(F,viewer);CHKERRQ(ierr);
294 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
295 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr);
296 ierr = MatShellGetContext(F,(void**)&fetidpmat_ctx);CHKERRQ(ierr);
297 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP TESTS--------------\n");CHKERRQ(ierr);
298 ierr = PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");CHKERRQ(ierr);
299 ierr = PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");CHKERRQ(ierr);
300 if (fetidp->fully_redundant) {
301 ierr = PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");CHKERRQ(ierr);
302 } else {
303 ierr = PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");CHKERRQ(ierr);
304 }
305 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
306
307 /* Get Vertices used to define the BDDC */
308 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
309 ierr = ISGetLocalSize(isvert,&n_vertices);CHKERRQ(ierr);
310 ierr = ISGetIndices(isvert,&vertex_indices);CHKERRQ(ierr);
311
312 /******************************************************************/
313 /* TEST A/B: Test numbering of global fetidp dofs */
314 /******************************************************************/
315 ierr = MatCreateVecs(F,&fetidp_global,NULL);CHKERRQ(ierr);
316 ierr = VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);CHKERRQ(ierr);
317 ierr = VecSet(fetidp_global,1.0);CHKERRQ(ierr);
318 ierr = VecSet(test_vec,1.);CHKERRQ(ierr);
319 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
320 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
321 if (fetidpmat_ctx->l2g_p) {
322 ierr = VecDuplicate(fetidpmat_ctx->vP,&test_vec_p);CHKERRQ(ierr);
323 ierr = VecSet(test_vec_p,1.);CHKERRQ(ierr);
324 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
325 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
326 }
327 ierr = VecAXPY(test_vec,-1.0,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
328 ierr = VecNorm(test_vec,NORM_INFINITY,&val);CHKERRQ(ierr);
329 ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
330 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
331 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc: % 1.14e\n",rval);CHKERRQ(ierr);
332
333 if (fetidpmat_ctx->l2g_p) {
334 ierr = VecAXPY(test_vec_p,-1.0,fetidpmat_ctx->vP);CHKERRQ(ierr);
335 ierr = VecNorm(test_vec_p,NORM_INFINITY,&val);CHKERRQ(ierr);
336 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
337 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc (p): % 1.14e\n",rval);CHKERRQ(ierr);
338 }
339
340 if (fetidp->fully_redundant) {
341 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
342 ierr = VecSet(fetidpmat_ctx->lambda_local,0.5);CHKERRQ(ierr);
343 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
344 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
345 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr);
346 val = PetscRealPart(sval)-fetidpmat_ctx->n_lambda;
347 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
348 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob: % 1.14e\n",rval);CHKERRQ(ierr);
349 }
350
351 if (fetidpmat_ctx->l2g_p) {
352 ierr = VecSet(pcis->vec1_N,1.0);CHKERRQ(ierr);
353 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
354 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
355 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
356
357 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
358 ierr = VecSet(fetidpmat_ctx->vP,-1.0);CHKERRQ(ierr);
359 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
360 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
361 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
362 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
363 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
364 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
365 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr);
366 val = PetscRealPart(sval);
367 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
368 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob (p): % 1.14e\n",rval);CHKERRQ(ierr);
369 }
370
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 /******************************************************************/
375
376 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
377 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr);
378 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
379 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
380 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
381 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
382 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
383 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
384 /* Action of B_delta */
385 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
386 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
387 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
388 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
389 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr);
390 ierr = PetscViewerASCIIPrintf(viewer,"C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",val);CHKERRQ(ierr);
391
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 /******************************************************************/
398
399 /* compute a random vector in \widetilde{W} */
400 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
401 /* set zero at vertices and essential dofs */
402 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
403 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
404 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirdofs);CHKERRQ(ierr);
405 if (dirdofs) {
406 const PetscInt *idxs;
407 PetscInt ndir;
408
409 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr);
410 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr);
411 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
412 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr);
413 }
414 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
415 /* store w for final comparison */
416 ierr = VecDuplicate(pcis->vec1_B,&test_vec);CHKERRQ(ierr);
417 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
418 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
419
420 /* Jump operator P_D : results stored in pcis->vec1_B */
421 /* Action of B_delta */
422 ierr = MatMult(fetidpmat_ctx->B_delta,test_vec,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
423 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
424 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
425 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
426 /* Action of B_Ddelta^T */
427 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
428 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
429 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
430
431 /* Average operator E_D : results stored in pcis->vec2_B */
432 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,test_vec,pcis->vec1_global);CHKERRQ(ierr);
433 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
434 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
435
436 /* test E_D=I-P_D */
437 ierr = VecAXPY(pcis->vec1_B,1.0,pcis->vec2_B);CHKERRQ(ierr);
438 ierr = VecAXPY(pcis->vec1_B,-1.0,test_vec);CHKERRQ(ierr);
439 ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&val);CHKERRQ(ierr);
440 ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
441 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr);
442 ierr = PetscViewerASCIIPrintf(viewer,"D: CHECK infty norm of E_D + P_D - I: % 1.14e\n",PetscGlobalRank,val);CHKERRQ(ierr);
443
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 /******************************************************************/
448
449 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr);
450 /* set zero at vertices and essential dofs */
451 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr);
452 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0;
453 if (dirdofs) {
454 const PetscInt *idxs;
455 PetscInt ndir;
456
457 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr);
458 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr);
459 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0;
460 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr);
461 }
462 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr);
463
464 /* Jump operator P_D : results stored in pcis->vec1_B */
465
466 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
467 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
468 /* Action of B_delta */
469 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
470 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr);
471 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
472 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
473 /* Action of B_Ddelta^T */
474 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
475 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
476 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
477 /* scaling */
478 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr);
479 ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&val);CHKERRQ(ierr);
480 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of R^T_D P_D: % 1.14e\n",val);CHKERRQ(ierr);
481
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 ierr = VecDuplicate(fetidp_global,&test_vec);CHKERRQ(ierr);
488 ierr = VecSetRandom(fetidp_global,NULL);CHKERRQ(ierr);
489 if (fetidpmat_ctx->l2g_p) {
490 ierr = VecSet(fetidpmat_ctx->vP,0.);CHKERRQ(ierr);
491 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
492 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
493 }
494 /* Action of B_Ddelta^T */
495 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
496 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
497 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr);
498 /* Action of B_delta */
499 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr);
500 ierr = VecSet(test_vec,0.0);CHKERRQ(ierr);
501 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
502 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
503 ierr = VecAXPY(fetidp_global,-1.,test_vec);CHKERRQ(ierr);
504 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr);
505 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of P^T_D - I: % 1.14e\n",val);CHKERRQ(ierr);
506 ierr = VecDestroy(&test_vec);CHKERRQ(ierr);
507 }
508 ierr = PetscViewerASCIIPrintf(viewer,"-------------------------------------\n");CHKERRQ(ierr);
509 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
510 ierr = VecDestroy(&test_vec_p);CHKERRQ(ierr);
511 ierr = ISDestroy(&dirdofs);CHKERRQ(ierr);
512 ierr = VecDestroy(&fetidp_global);CHKERRQ(ierr);
513 ierr = ISRestoreIndices(isvert,&vertex_indices);CHKERRQ(ierr);
514 ierr = PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr);
515 PetscFunctionReturn(0);
516 }
517
KSPFETIDPSetUpOperators(KSP ksp)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;
535
536 PetscFunctionBegin;
537 pisz = PETSC_FALSE;
538 flip = PETSC_FALSE;
539 allp = PETSC_FALSE;
540 schp = PETSC_FALSE;
541 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");CHKERRQ(ierr);
542 ierr = PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);CHKERRQ(ierr);
543 ierr = PetscOptionsBool("-ksp_fetidp_pressure_all","Use the whole pressure set instead of just that at the interface",NULL,allp,&allp,NULL);CHKERRQ(ierr);
544 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);CHKERRQ(ierr);
545 ierr = PetscOptionsBool("-ksp_fetidp_pressure_schur","Use a BDDC solver for pressure",NULL,schp,&schp,NULL);CHKERRQ(ierr);
546 ierr = PetscOptionsEnd();CHKERRQ(ierr);
547
548 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)ksp),&size);CHKERRQ(ierr);
549 fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint);
550 if (size == 1) fetidp->saddlepoint = PETSC_FALSE;
551
552 ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr);
553 ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr);
554 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Amat should be of type MATIS");
555
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 ierr = PetscObjectStateGet((PetscObject)A,&matstate);CHKERRQ(ierr);
560 ierr = MatGetNonzeroState(A,&matnnzstate);CHKERRQ(ierr);
561 if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(0);
562 fetidp->matstate = matstate;
563 fetidp->matnnzstate = matnnzstate;
564 fetidp->statechanged = fetidp->saddlepoint;
565
566 /* see if we have some fields attached */
567 if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
568 DM dm;
569 PetscContainer c;
570
571 ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
572 ierr = PetscObjectQuery((PetscObject)A,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr);
573 if (dm) {
574 IS *fields;
575 PetscInt nf,i;
576
577 ierr = DMCreateFieldDecomposition(dm,&nf,NULL,&fields,NULL);CHKERRQ(ierr);
578 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,nf,fields);CHKERRQ(ierr);
579 for (i=0;i<nf;i++) {
580 ierr = ISDestroy(&fields[i]);CHKERRQ(ierr);
581 }
582 ierr = PetscFree(fields);CHKERRQ(ierr);
583 } else if (c) {
584 MatISLocalFields lf;
585
586 ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr);
587 ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr);
588 }
589 }
590
591 if (!fetidp->saddlepoint) {
592 ierr = PCSetOperators(fetidp->innerbddc,A,A);CHKERRQ(ierr);
593 } else {
594 Mat nA,lA,PPmat;
595 MatNullSpace nnsp;
596 IS pP;
597 PetscInt totP;
598
599 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr);
600 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr);
601
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;
613
614 ierr = MatGetLocalSize(A,&nl,NULL);CHKERRQ(ierr);
615 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr);
616 ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr);
617 ierr = MatGetLocalToGlobalMapping(A,&l2g,NULL);CHKERRQ(ierr);
618
619 if (!pcis->is_I_local) { /* need to compute interior dofs */
620 ierr = PetscCalloc1(n,&count);CHKERRQ(ierr);
621 ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
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 ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr);
627 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr);
628 } else {
629 ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr);
630 II = pcis->is_I_local;
631 }
632
633 /* interior dofs in layout */
634 ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
635 ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
636 ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr);
637 ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr);
638 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
639 ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr);
640 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
641 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
642 ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr);
643 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
644 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr);
645
646 /* pressure dofs */
647 Pall = NULL;
648 lPall = NULL;
649 ploc = PETSC_FALSE;
650 if (fid < 0) { /* zero pressure block */
651 PetscInt np;
652
653 ierr = MatFindZeroDiagonals(A,&Pall);CHKERRQ(ierr);
654 ierr = ISGetSize(Pall,&np);CHKERRQ(ierr);
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 ierr = ISDestroy(&Pall);CHKERRQ(ierr);
658 } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
659 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,1,&Pall);CHKERRQ(ierr);
660 }
661 }
662 if (!Pall) { /* look for registered fields */
663 if (pcbddc->n_ISForDofsLocal) {
664 PetscInt np;
665
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 ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr);
669 ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
670 ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
671 ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr);
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 ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr);
676 Pall = pcbddc->ISForDofs[fid];
677 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal");
678 }
679
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 ierr = ISDifference(II,lPall,&nII);CHKERRQ(ierr);
686 ierr = ISDestroy(&II);CHKERRQ(ierr);
687 II = nII;
688 } else {
689 IS nII;
690 ierr = ISDifference(pII,Pall,&nII);CHKERRQ(ierr);
691 ierr = ISDestroy(&pII);CHKERRQ(ierr);
692 pII = nII;
693 }
694 }
695 if (ploc) {
696 ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr);
697 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
698 } else {
699 ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr);
700 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
701 /* need all local pressure dofs */
702 ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
703 ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
704 ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr);
705 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
706 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
707 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
708 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
709 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
710 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
711 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr);
712 }
713
714 if (!Pall) {
715 ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
716 ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
717 ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr);
718 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr);
719 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1;
720 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr);
721 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
722 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
723 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
724 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr);
725 }
726 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject)Pall);CHKERRQ(ierr);
727
728 if (flip) {
729 PetscInt npl;
730 ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr);
731 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr);
732 ierr = MatCreateVecs(A,NULL,&fetidp->rhs_flip);CHKERRQ(ierr);
733 ierr = VecSet(fetidp->rhs_flip,1.);CHKERRQ(ierr);
734 ierr = VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
735 for (i=0;i<npl;i++) {
736 ierr = VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr);
737 }
738 ierr = VecAssemblyBegin(fetidp->rhs_flip);CHKERRQ(ierr);
739 ierr = VecAssemblyEnd(fetidp->rhs_flip);CHKERRQ(ierr);
740 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);CHKERRQ(ierr);
741 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr);
742 }
743 ierr = ISDestroy(&Pall);CHKERRQ(ierr);
744 ierr = ISDestroy(&pII);CHKERRQ(ierr);
745
746 /* local selected pressures in subdomain-wise and global ordering */
747 ierr = PetscArrayzero(matis->sf_leafdata,n);CHKERRQ(ierr);
748 ierr = PetscArrayzero(matis->sf_rootdata,nl);CHKERRQ(ierr);
749 if (!ploc) {
750 PetscInt *widxs2;
751
752 if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing parallel pressure IS");
753 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr);
754 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr);
755 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1;
756 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr);
757 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
758 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr);
759 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i;
760 ierr = PetscMalloc1(ni,&widxs2);CHKERRQ(ierr);
761 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs2);CHKERRQ(ierr);
762 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr);
763 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr);
764 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs2,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr);
765 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
766 ierr = ISDestroy(&is1);CHKERRQ(ierr);
767 } else {
768 if (!lP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing sequential pressure IS");
769 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr);
770 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
771 for (i=0;i<ni;i++)
772 if (idxs[i] >=0 && idxs[i] < n)
773 matis->sf_leafdata[idxs[i]] = 1;
774 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
775 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
776 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr);
777 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
778 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr);
779 ierr = ISDestroy(&is1);CHKERRQ(ierr);
780 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr);
781 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst;
782 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr);
783 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr);
784 }
785 ierr = PetscFree(widxs);CHKERRQ(ierr);
786
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;
794
795 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr);
796 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr);
797 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr);
798 ierr = ISDestroy(&is1);CHKERRQ(ierr);
799 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g_t);CHKERRQ(ierr);
800 ierr = ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr);
801 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr);
802 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr);
803 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iV",i,j);
804 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);CHKERRQ(ierr);
805 ierr = ISLocalToGlobalMappingDestroy(&l2g_t);CHKERRQ(ierr);
806 ierr = ISDestroy(&is1);CHKERRQ(ierr);
807 ierr = ISDestroy(&is2);CHKERRQ(ierr);
808 }
809 ierr = ISDestroy(&II);CHKERRQ(ierr);
810
811 /* exclude selected pressures from the inner BDDC */
812 if (pcbddc->DirichletBoundariesLocal) {
813 IS list[2],plP,isout;
814 PetscInt np;
815
816 /* need a parallel IS */
817 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
818 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
819 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr);
820 list[0] = plP;
821 list[1] = pcbddc->DirichletBoundariesLocal;
822 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
823 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr);
824 ierr = ISDestroy(&plP);CHKERRQ(ierr);
825 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
826 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr);
827 ierr = ISDestroy(&isout);CHKERRQ(ierr);
828 } else if (pcbddc->DirichletBoundaries) {
829 IS list[2],isout;
830
831 list[0] = pP;
832 list[1] = pcbddc->DirichletBoundaries;
833 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr);
834 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr);
835 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr);
836 ierr = ISDestroy(&isout);CHKERRQ(ierr);
837 } else {
838 IS plP;
839 PetscInt np;
840
841 /* need a parallel IS */
842 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
843 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr);
844 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr);
845 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr);
846 ierr = ISDestroy(&plP);CHKERRQ(ierr);
847 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr);
848 }
849
850 /* save CSR information for the pressure BDDC solver (if any) */
851 if (schp) {
852 PetscInt np,nt;
853
854 ierr = MatGetSize(matis->A,&nt,NULL);CHKERRQ(ierr);
855 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr);
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;
860
861 if (nv && nv == nt) {
862 ISLocalToGlobalMapping pmap;
863 PetscInt *schp_csr,*schp_xadj,*schp_adjn,p;
864 PetscContainer c;
865
866 ierr = ISLocalToGlobalMappingCreateIS(lPall,&pmap);CHKERRQ(ierr);
867 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr);
868 for (p = 0, nv = 0; p < np; p++) {
869 PetscInt x,n = idxs[p];
870
871 ierr = ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,NULL);CHKERRQ(ierr);
872 nv += x;
873 }
874 ierr = PetscMalloc1(np + 1 + nv,&schp_csr);CHKERRQ(ierr);
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];
879
880 ierr = ISGlobalToLocalMappingApply(pmap,IS_GTOLM_DROP,xadj[n+1]-xadj[n],adjn+xadj[n],&x,schp_adjn + schp_xadj[p]);CHKERRQ(ierr);
881 schp_xadj[p+1] = schp_xadj[p] + x;
882 }
883 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr);
884 ierr = ISLocalToGlobalMappingDestroy(&pmap);CHKERRQ(ierr);
885 ierr = PetscContainerCreate(PETSC_COMM_SELF,&c);CHKERRQ(ierr);
886 ierr = PetscContainerSetPointer(c,schp_csr);CHKERRQ(ierr);
887 ierr = PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault);CHKERRQ(ierr);
888 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pCSR",(PetscObject)c);CHKERRQ(ierr);
889 ierr = PetscContainerDestroy(&c);CHKERRQ(ierr);
890
891 }
892 }
893 }
894 ierr = ISDestroy(&lPall);CHKERRQ(ierr);
895 ierr = ISDestroy(&lP);CHKERRQ(ierr);
896 fetidp->pP = pP;
897 }
898
899 /* total number of selected pressure dofs */
900 ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr);
901
902 /* Set operator for inner BDDC */
903 if (totP || fetidp->rhs_flip) {
904 ierr = MatDuplicate(A,MAT_COPY_VALUES,&nA);CHKERRQ(ierr);
905 } else {
906 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
907 nA = A;
908 }
909 if (fetidp->rhs_flip) {
910 ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr);
911 if (totP) {
912 Mat lA2;
913
914 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr);
915 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr);
916 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr);
917 ierr = MatDestroy(&lA2);CHKERRQ(ierr);
918 }
919 }
920
921 if (totP) {
922 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
923 ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr);
924 } else {
925 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr);
926 }
927 ierr = MatGetNearNullSpace(Ap,&nnsp);CHKERRQ(ierr);
928 if (!nnsp) {
929 ierr = MatGetNullSpace(Ap,&nnsp);CHKERRQ(ierr);
930 }
931 if (!nnsp) {
932 ierr = MatGetNearNullSpace(A,&nnsp);CHKERRQ(ierr);
933 }
934 if (!nnsp) {
935 ierr = MatGetNullSpace(A,&nnsp);CHKERRQ(ierr);
936 }
937 ierr = MatSetNearNullSpace(nA,nnsp);CHKERRQ(ierr);
938 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr);
939 ierr = MatDestroy(&nA);CHKERRQ(ierr);
940
941 /* non-zero rhs on interior dofs when applying the preconditioner */
942 if (totP) pcbddc->switch_static = PETSC_TRUE;
943
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 }
949
950 /* Operators for pressure preconditioner */
951 if (totP) {
952 /* Extract pressure block if needed */
953 if (!pisz) {
954 Mat C;
955 IS nzrows = NULL;
956
957 ierr = MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
958 ierr = MatFindNonzeroRows(C,&nzrows);CHKERRQ(ierr);
959 if (nzrows) {
960 PetscInt i;
961
962 ierr = ISGetSize(nzrows,&i);CHKERRQ(ierr);
963 ierr = ISDestroy(&nzrows);CHKERRQ(ierr);
964 if (!i) pisz = PETSC_TRUE;
965 }
966 if (!pisz) {
967 ierr = MatScale(C,-1.);CHKERRQ(ierr); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized */
968 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr);
969 }
970 ierr = MatDestroy(&C);CHKERRQ(ierr);
971 }
972 /* Divergence mat */
973 if (!pcbddc->divudotp) {
974 Mat B;
975 IS P;
976 IS l2l = NULL;
977 PetscBool save;
978
979 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr);
980 if (!pisz) {
981 IS F,V;
982 PetscInt m,M;
983
984 ierr = MatGetOwnershipRange(A,&m,&M);CHKERRQ(ierr);
985 ierr = ISCreateStride(PetscObjectComm((PetscObject)A),M-m,m,1,&F);CHKERRQ(ierr);
986 ierr = ISComplement(P,m,M,&V);CHKERRQ(ierr);
987 ierr = MatCreateSubMatrix(A,P,V,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
988 {
989 Mat_IS *Bmatis = (Mat_IS*)B->data;
990 ierr = PetscObjectReference((PetscObject)Bmatis->getsub_cis);CHKERRQ(ierr);
991 l2l = Bmatis->getsub_cis;
992 }
993 ierr = ISDestroy(&V);CHKERRQ(ierr);
994 ierr = ISDestroy(&F);CHKERRQ(ierr);
995 } else {
996 ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
997 }
998 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */
999 ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,PETSC_FALSE,l2l);CHKERRQ(ierr);
1000 pcbddc->compute_nonetflux = save;
1001 ierr = MatDestroy(&B);CHKERRQ(ierr);
1002 ierr = ISDestroy(&l2l);CHKERRQ(ierr);
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 ierr = KSPFETIDPSetPressureOperator(ksp,Ap);CHKERRQ(ierr);
1007 }
1008 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr);
1009
1010 /* PPmat not present, use some default choice */
1011 if (!PPmat) {
1012 Mat C;
1013
1014 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr);
1015 if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */
1016 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1017 } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */
1018 IS P;
1019
1020 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr);
1021 ierr = MatCreateSubMatrix(A,P,P,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1022 ierr = MatScale(C,-1.);CHKERRQ(ierr);
1023 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1024 ierr = MatDestroy(&C);CHKERRQ(ierr);
1025 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */
1026 PetscInt nl;
1027
1028 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr);
1029 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&C);CHKERRQ(ierr);
1030 ierr = MatSetSizes(C,nl,nl,totP,totP);CHKERRQ(ierr);
1031 ierr = MatSetType(C,MATAIJ);CHKERRQ(ierr);
1032 ierr = MatMPIAIJSetPreallocation(C,1,NULL,0,NULL);CHKERRQ(ierr);
1033 ierr = MatSeqAIJSetPreallocation(C,1,NULL);CHKERRQ(ierr);
1034 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1035 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1036 ierr = MatShift(C,1.);CHKERRQ(ierr);
1037 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1038 ierr = MatDestroy(&C);CHKERRQ(ierr);
1039 }
1040 }
1041
1042 /* Preconditioned operator for the pressure block */
1043 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr);
1044 if (PPmat) {
1045 Mat C;
1046 IS Pall;
1047 PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg;
1048
1049 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr);
1050 ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr);
1051 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr);
1052 ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr);
1053 ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr);
1054 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr);
1055 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
1056 ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr);
1057 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr);
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 ierr = MatCreateSubMatrix(PPmat,Pall,Pall,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1065 } else {
1066 ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1067 }
1068 } else if (pAg == PAM) { /* global ordering for pressure only */
1069 if (!allp && !schp) { /* solving for interface pressure only */
1070 IS restr;
1071
1072 ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr);
1073 ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
1074 ierr = ISDestroy(&restr);CHKERRQ(ierr);
1075 } else {
1076 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr);
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 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr);
1082 C = PPmat;
1083 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix");
1084
1085 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr);
1086 ierr = MatDestroy(&C);CHKERRQ(ierr);
1087 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing Pmat for pressure block");
1088 } else { /* totP == 0 */
1089 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr);
1090 }
1091 }
1092 PetscFunctionReturn(0);
1093 }
1094
KSPSetUp_FETIDP(KSP ksp)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;
1100 PetscErrorCode ierr;
1101
1102 PetscFunctionBegin;
1103 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr);
1104 /* set up BDDC */
1105 ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr);
1106 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr);
1107 /* FETI-DP as it is implemented needs an exact coarse solver */
1108 if (pcbddc->coarse_ksp) {
1109 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
1110 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr);
1111 }
1112 /* FETI-DP as it is implemented needs exact local Neumann solvers */
1113 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr);
1114 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr);
1115
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 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr);
1125 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr);
1126 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr);
1127 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr);
1128 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr);
1129 ierr = PetscObjectIncrementTabLevel((PetscObject)D,(PetscObject)fetidp->innerksp,0);CHKERRQ(ierr);
1130 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr);
1131 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr);
1132 ierr = MatDestroy(&F);CHKERRQ(ierr);
1133 ierr = PCDestroy(&D);CHKERRQ(ierr);
1134 if (fetidp->check) {
1135 PetscViewer viewer;
1136
1137 if (!pcbddc->dbg_viewer) {
1138 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
1139 } else {
1140 viewer = pcbddc->dbg_viewer;
1141 }
1142 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr);
1143 }
1144 }
1145 fetidp->statechanged = PETSC_FALSE;
1146 pcbddc->new_primal_space = PETSC_FALSE;
1147
1148 /* propagate settings to the inner solve */
1149 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr);
1150 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr);
1151 if (ksp->res_hist) {
1152 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr);
1153 }
1154 ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr);
1155 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr);
1156 PetscFunctionReturn(0);
1157 }
1158
KSPSolve_FETIDP(KSP ksp)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;
1170
1171 PetscFunctionBegin;
1172 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
1173 if (fetidp->saddlepoint) {
1174 ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr);
1175 }
1176 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr);
1177 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr);
1178 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr);
1179 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr);
1180 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr);
1181 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr);
1182 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr);
1183 if (ksp->transpose_solve) {
1184 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
1185 } else {
1186 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr);
1187 }
1188 ierr = KSPGetConvergedReason(fetidp->innerksp,&reason);CHKERRQ(ierr);
1189 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr);
1190 ierr = PCGetFailedReason(pc,&pcreason);CHKERRQ(ierr);
1191 if ((reason < 0 && reason != KSP_DIVERGED_ITS) || pcreason) {
1192 PetscInt its;
1193 ierr = KSPGetIterationNumber(fetidp->innerksp,&its);CHKERRQ(ierr);
1194 ksp->reason = KSP_DIVERGED_PC_FAILED;
1195 ierr = VecSetInf(Xl);CHKERRQ(ierr);
1196 ierr = PetscInfo3(ksp,"Inner KSP solve failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);CHKERRQ(ierr);
1197 }
1198 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr);
1199 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr);
1200 if (nsp) {
1201 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr);
1202 }
1203 /* update ksp with stats from inner ksp */
1204 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr);
1205 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr);
1206 ksp->totalits += ksp->its;
1207 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1213 }
1214
KSPReset_FETIDP(KSP ksp)1215 static PetscErrorCode KSPReset_FETIDP(KSP ksp)
1216 {
1217 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
1218 PC_BDDC *pcbddc;
1219 PetscErrorCode ierr;
1220
1221 PetscFunctionBegin;
1222 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr);
1223 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr);
1224 /* avoid PCReset that does not take into account ref counting */
1225 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
1226 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr);
1227 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr);
1228 pcbddc = (PC_BDDC*)fetidp->innerbddc->data;
1229 pcbddc->symmetric_primal = PETSC_FALSE;
1230 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr);
1231 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr);
1232 fetidp->saddlepoint = PETSC_FALSE;
1233 fetidp->matstate = -1;
1234 fetidp->matnnzstate = -1;
1235 fetidp->statechanged = PETSC_TRUE;
1236 PetscFunctionReturn(0);
1237 }
1238
KSPDestroy_FETIDP(KSP ksp)1239 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
1240 {
1241 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
1242 PetscErrorCode ierr;
1243
1244 PetscFunctionBegin;
1245 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr);
1246 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr);
1247 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr);
1248 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr);
1249 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr);
1250 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr);
1251 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr);
1252 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr);
1253 ierr = PetscFree(ksp->data);CHKERRQ(ierr);
1254 PetscFunctionReturn(0);
1255 }
1256
KSPView_FETIDP(KSP ksp,PetscViewer viewer)1257 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer)
1258 {
1259 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
1260 PetscErrorCode ierr;
1261 PetscBool iascii;
1262
1263 PetscFunctionBegin;
1264 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1265 if (iascii) {
1266 ierr = PetscViewerASCIIPrintf(viewer," fully redundant: %d\n",fetidp->fully_redundant);CHKERRQ(ierr);
1267 ierr = PetscViewerASCIIPrintf(viewer," saddle point: %d\n",fetidp->saddlepoint);CHKERRQ(ierr);
1268 ierr = PetscViewerASCIIPrintf(viewer,"Inner KSP solver details\n");CHKERRQ(ierr);
1269 }
1270 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1271 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr);
1272 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1273 if (iascii) {
1274 ierr = PetscViewerASCIIPrintf(viewer,"Inner BDDC solver details\n");CHKERRQ(ierr);
1275 }
1276 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1277 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr);
1278 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1279 PetscFunctionReturn(0);
1280 }
1281
KSPSetFromOptions_FETIDP(PetscOptionItems * PetscOptionsObject,KSP ksp)1282 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp)
1283 {
1284 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data;
1285 PetscErrorCode ierr;
1286
1287 PetscFunctionBegin;
1288 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
1289 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
1290 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr);
1291 if (!fetidp->userbddc) {
1292 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
1293 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr);
1294 }
1295 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr);
1296 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr);
1297 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr);
1298 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr);
1299 ierr = PetscOptionsTail();CHKERRQ(ierr);
1300 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr);
1301 PetscFunctionReturn(0);
1302 }
1303
1304 /*MC
1305 KSPFETIDP - The FETI-DP method
1306
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.
1310
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.
1324
1325 Level: Advanced
1326
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.
1333
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
1346
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.
1348
1349 The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution.
1350
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.
1354
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
1360
1361 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP()
1362 M*/
KSPCreate_FETIDP(KSP ksp)1363 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
1364 {
1365 PetscErrorCode ierr;
1366 KSP_FETIDP *fetidp;
1367 KSP_FETIDPMon *monctx;
1368 PC_BDDC *pcbddc;
1369 PC pc;
1370
1371 PetscFunctionBegin;
1372 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);CHKERRQ(ierr);
1373 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);CHKERRQ(ierr);
1374 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
1375 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr);
1376 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
1377 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr);
1378 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr);
1379
1380 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr);
1381 fetidp->matstate = -1;
1382 fetidp->matnnzstate = -1;
1383 fetidp->statechanged = PETSC_TRUE;
1384
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 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1396 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1397 /* create the inner KSP for the Lagrange multipliers */
1398 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr);
1399 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr);
1400 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
1401 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr);
1402 /* monitor */
1403 ierr = PetscNew(&monctx);CHKERRQ(ierr);
1404 monctx->parentksp = ksp;
1405 fetidp->monctx = monctx;
1406 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr);
1407 /* create the inner BDDC */
1408 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr);
1409 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr);
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 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr);
1415 /* composed functions */
1416 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr);
1417 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr);
1418 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr);
1419 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);CHKERRQ(ierr);
1420 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
1421 ksp->setupnewmatrix = PETSC_TRUE;
1422 PetscFunctionReturn(0);
1423 }
1424