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 = {http://www.worldscientific.com/doi/abs/10.1142/S0218202513500632},\n" 30 "eprint = {http://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 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 @*/ 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 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 @*/ 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 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 @*/ 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 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 @*/ 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 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 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 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 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 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 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 = MatISSetUpSF(A);CHKERRQ(ierr); 635 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 636 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 637 ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr); 638 ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr); 639 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 640 ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr); 641 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 642 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 643 ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr); 644 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 645 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr); 646 647 /* pressure dofs */ 648 Pall = NULL; 649 lPall = NULL; 650 ploc = PETSC_FALSE; 651 if (fid < 0) { /* zero pressure block */ 652 PetscInt np; 653 654 ierr = MatFindZeroDiagonals(A,&Pall);CHKERRQ(ierr); 655 ierr = ISGetSize(Pall,&np);CHKERRQ(ierr); 656 if (!np) { /* zero-block not found, defaults to last field (if set) */ 657 fid = pcbddc->n_ISForDofsLocal ? pcbddc->n_ISForDofsLocal - 1 : pcbddc->n_ISForDofs - 1; 658 ierr = ISDestroy(&Pall);CHKERRQ(ierr); 659 } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) { 660 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,1,&Pall);CHKERRQ(ierr); 661 } 662 } 663 if (!Pall) { /* look for registered fields */ 664 if (pcbddc->n_ISForDofsLocal) { 665 PetscInt np; 666 667 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); 668 /* need a sequential IS */ 669 ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr); 670 ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 671 ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 672 ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 673 ploc = PETSC_TRUE; 674 } else if (pcbddc->n_ISForDofs) { 675 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); 676 ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr); 677 Pall = pcbddc->ISForDofs[fid]; 678 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal"); 679 } 680 681 /* if the user requested the entire pressure, 682 remove the interior pressure dofs from II (or pII) */ 683 if (allp) { 684 if (ploc) { 685 IS nII; 686 ierr = ISDifference(II,lPall,&nII);CHKERRQ(ierr); 687 ierr = ISDestroy(&II);CHKERRQ(ierr); 688 II = nII; 689 } else { 690 IS nII; 691 ierr = ISDifference(pII,Pall,&nII);CHKERRQ(ierr); 692 ierr = ISDestroy(&pII);CHKERRQ(ierr); 693 pII = nII; 694 } 695 } 696 if (ploc) { 697 ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr); 698 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 699 } else { 700 ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr); 701 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 702 /* need all local pressure dofs */ 703 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 704 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 705 ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr); 706 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 707 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 708 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 709 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 710 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 711 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 712 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 713 } 714 715 if (!Pall) { 716 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 717 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 718 ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr); 719 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr); 720 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 721 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr); 722 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 723 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 724 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 725 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr); 726 } 727 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject)Pall);CHKERRQ(ierr); 728 729 if (flip) { 730 PetscInt npl; 731 ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr); 732 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 733 ierr = MatCreateVecs(A,NULL,&fetidp->rhs_flip);CHKERRQ(ierr); 734 ierr = VecSet(fetidp->rhs_flip,1.);CHKERRQ(ierr); 735 ierr = VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 736 for (i=0;i<npl;i++) { 737 ierr = VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr); 738 } 739 ierr = VecAssemblyBegin(fetidp->rhs_flip);CHKERRQ(ierr); 740 ierr = VecAssemblyEnd(fetidp->rhs_flip);CHKERRQ(ierr); 741 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);CHKERRQ(ierr); 742 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 743 } 744 ierr = ISDestroy(&Pall);CHKERRQ(ierr); 745 ierr = ISDestroy(&pII);CHKERRQ(ierr); 746 747 /* local selected pressures in subdomain-wise and global ordering */ 748 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 749 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 750 if (!ploc) { 751 PetscInt *widxs2; 752 753 if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing parallel pressure IS"); 754 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr); 755 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr); 756 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 757 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr); 758 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 759 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 760 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 761 ierr = PetscMalloc1(ni,&widxs2);CHKERRQ(ierr); 762 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs2);CHKERRQ(ierr); 763 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr); 764 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 765 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs2,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 766 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 767 ierr = ISDestroy(&is1);CHKERRQ(ierr); 768 } else { 769 if (!lP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing sequential pressure IS"); 770 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr); 771 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 772 for (i=0;i<ni;i++) 773 if (idxs[i] >=0 && idxs[i] < n) 774 matis->sf_leafdata[idxs[i]] = 1; 775 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 776 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 777 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr); 778 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 779 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 780 ierr = ISDestroy(&is1);CHKERRQ(ierr); 781 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 782 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 783 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr); 784 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 785 } 786 ierr = PetscFree(widxs);CHKERRQ(ierr); 787 788 /* If there's any "interior pressure", 789 we may want to use a discrete harmonic solver instead 790 of a Stokes harmonic for the Dirichlet preconditioner 791 Need to extract the interior velocity dofs in interior dofs ordering (iV) 792 and interior pressure dofs in local ordering (iP) */ 793 if (!allp) { 794 ISLocalToGlobalMapping l2g_t; 795 796 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr); 797 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr); 798 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr); 799 ierr = ISDestroy(&is1);CHKERRQ(ierr); 800 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g_t);CHKERRQ(ierr); 801 ierr = ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr); 802 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr); 803 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr); 804 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iV",i,j); 805 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);CHKERRQ(ierr); 806 ierr = ISLocalToGlobalMappingDestroy(&l2g_t);CHKERRQ(ierr); 807 ierr = ISDestroy(&is1);CHKERRQ(ierr); 808 ierr = ISDestroy(&is2);CHKERRQ(ierr); 809 } 810 ierr = ISDestroy(&lPall);CHKERRQ(ierr); 811 ierr = ISDestroy(&II);CHKERRQ(ierr); 812 813 /* exclude selected pressures from the inner BDDC */ 814 if (pcbddc->DirichletBoundariesLocal) { 815 IS list[2],plP,isout; 816 PetscInt np; 817 818 /* need a parallel IS */ 819 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 820 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 821 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr); 822 list[0] = plP; 823 list[1] = pcbddc->DirichletBoundariesLocal; 824 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 825 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr); 826 ierr = ISDestroy(&plP);CHKERRQ(ierr); 827 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 828 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr); 829 ierr = ISDestroy(&isout);CHKERRQ(ierr); 830 } else if (pcbddc->DirichletBoundaries) { 831 IS list[2],isout; 832 833 list[0] = pP; 834 list[1] = pcbddc->DirichletBoundaries; 835 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 836 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr); 837 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr); 838 ierr = ISDestroy(&isout);CHKERRQ(ierr); 839 } else { 840 IS plP; 841 PetscInt np; 842 843 /* need a parallel IS */ 844 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 845 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 846 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr); 847 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr); 848 ierr = ISDestroy(&plP);CHKERRQ(ierr); 849 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 850 } 851 ierr = ISDestroy(&lP);CHKERRQ(ierr); 852 fetidp->pP = pP; 853 } 854 855 /* total number of selected pressure dofs */ 856 ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr); 857 858 /* Set operator for inner BDDC */ 859 if (totP || fetidp->rhs_flip) { 860 ierr = MatDuplicate(A,MAT_COPY_VALUES,&nA);CHKERRQ(ierr); 861 } else { 862 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 863 nA = A; 864 } 865 if (fetidp->rhs_flip) { 866 ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr); 867 if (totP) { 868 Mat lA2; 869 870 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr); 871 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr); 872 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr); 873 ierr = MatDestroy(&lA2);CHKERRQ(ierr); 874 } 875 } 876 877 if (totP) { 878 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 879 ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr); 880 } else { 881 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr); 882 } 883 ierr = MatGetNearNullSpace(Ap,&nnsp);CHKERRQ(ierr); 884 if (!nnsp) { 885 ierr = MatGetNullSpace(Ap,&nnsp);CHKERRQ(ierr); 886 } 887 if (!nnsp) { 888 ierr = MatGetNearNullSpace(A,&nnsp);CHKERRQ(ierr); 889 } 890 if (!nnsp) { 891 ierr = MatGetNullSpace(A,&nnsp);CHKERRQ(ierr); 892 } 893 ierr = MatSetNearNullSpace(nA,nnsp);CHKERRQ(ierr); 894 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr); 895 ierr = MatDestroy(&nA);CHKERRQ(ierr); 896 897 /* non-zero rhs on interior dofs when applying the preconditioner */ 898 if (totP) pcbddc->switch_static = PETSC_TRUE; 899 900 /* if there are no interface pressures, set inner bddc flag for benign saddle point */ 901 if (!totP) { 902 pcbddc->benign_saddle_point = PETSC_TRUE; 903 pcbddc->compute_nonetflux = PETSC_TRUE; 904 } 905 906 /* Operators for pressure preconditioner */ 907 if (totP) { 908 /* Extract pressure block if needed */ 909 if (!pisz) { 910 Mat C; 911 IS nzrows = NULL; 912 913 ierr = MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 914 ierr = MatFindNonzeroRows(C,&nzrows);CHKERRQ(ierr); 915 if (nzrows) { 916 PetscInt i; 917 918 ierr = ISGetSize(nzrows,&i);CHKERRQ(ierr); 919 ierr = ISDestroy(&nzrows);CHKERRQ(ierr); 920 if (!i) pisz = PETSC_TRUE; 921 } 922 if (!pisz) { 923 ierr = MatScale(C,-1.);CHKERRQ(ierr); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized */ 924 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr); 925 } 926 ierr = MatDestroy(&C);CHKERRQ(ierr); 927 } 928 /* Divergence mat */ 929 if (!pcbddc->divudotp) { 930 Mat B; 931 IS P; 932 PetscBool save; 933 934 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 935 if (!pisz) { 936 IS F,V; 937 PetscInt m,M; 938 939 ierr = MatGetOwnershipRange(A,&m,&M);CHKERRQ(ierr); 940 ierr = ISCreateStride(PetscObjectComm((PetscObject)P),M-m,m,1,&F);CHKERRQ(ierr); 941 ierr = ISComplement(P,m,M,&V);CHKERRQ(ierr); 942 ierr = MatCreateSubMatrix(A,F,V,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 943 ierr = MatZeroRowsIS(B,P,0.0,NULL,NULL);CHKERRQ(ierr); 944 ierr = ISDestroy(&V);CHKERRQ(ierr); 945 ierr = ISDestroy(&F);CHKERRQ(ierr); 946 } else { 947 ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 948 } 949 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */ 950 ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,(PetscBool)!pisz,NULL);CHKERRQ(ierr); 951 pcbddc->compute_nonetflux = save; 952 ierr = MatDestroy(&B);CHKERRQ(ierr); 953 } 954 if (A != Ap) { /* user has provided a different Pmat, this always superseeds the setter (TODO: is it OK?) */ 955 /* use monolithic operator, we restrict later */ 956 ierr = KSPFETIDPSetPressureOperator(ksp,Ap);CHKERRQ(ierr); 957 } 958 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 959 960 /* PPmat not present, use some default choice */ 961 if (!PPmat) { 962 Mat C; 963 964 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr); 965 if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */ 966 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 967 } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */ 968 IS P; 969 970 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 971 ierr = MatCreateSubMatrix(A,P,P,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 972 ierr = MatScale(C,-1.);CHKERRQ(ierr); 973 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 974 ierr = MatDestroy(&C);CHKERRQ(ierr); 975 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */ 976 PetscInt nl; 977 978 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr); 979 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&C);CHKERRQ(ierr); 980 ierr = MatSetSizes(C,nl,nl,totP,totP);CHKERRQ(ierr); 981 ierr = MatSetType(C,MATAIJ);CHKERRQ(ierr); 982 ierr = MatMPIAIJSetPreallocation(C,1,NULL,0,NULL);CHKERRQ(ierr); 983 ierr = MatSeqAIJSetPreallocation(C,1,NULL);CHKERRQ(ierr); 984 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 985 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 986 ierr = MatShift(C,1.);CHKERRQ(ierr); 987 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 988 ierr = MatDestroy(&C);CHKERRQ(ierr); 989 } 990 } 991 992 /* Preconditioned operator for the pressure block */ 993 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 994 if (PPmat) { 995 Mat C; 996 IS Pall; 997 PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg; 998 999 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr); 1000 ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr); 1001 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 1002 ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr); 1003 ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr); 1004 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 1005 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 1006 ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr); 1007 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 1008 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 1009 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 1010 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); 1011 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); 1012 if (PAM == AM) { /* monolithic ordering, restrict to pressure */ 1013 if (schp) { 1014 ierr = MatCreateSubMatrix(PPmat,Pall,Pall,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1015 } else { 1016 ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1017 } 1018 } else if (pAg == PAM) { /* global ordering for pressure only */ 1019 if (!allp && !schp) { /* solving for interface pressure only */ 1020 IS restr; 1021 1022 ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr); 1023 ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1024 ierr = ISDestroy(&restr);CHKERRQ(ierr); 1025 } else { 1026 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 1027 C = PPmat; 1028 } 1029 } else if (pIg == PAM) { /* global ordering for selected pressure only */ 1030 if (schp) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Need the entire matrix"); 1031 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 1032 C = PPmat; 1033 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix"); 1034 1035 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 1036 ierr = MatDestroy(&C);CHKERRQ(ierr); 1037 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing Pmat for pressure block"); 1038 } else { /* totP == 0 */ 1039 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr); 1040 } 1041 } 1042 PetscFunctionReturn(0); 1043 } 1044 1045 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 1046 { 1047 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1048 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1049 PetscBool flg; 1050 PetscErrorCode ierr; 1051 1052 PetscFunctionBegin; 1053 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 1054 /* set up BDDC */ 1055 ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr); 1056 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 1057 /* FETI-DP as it is implemented needs an exact coarse solver */ 1058 if (pcbddc->coarse_ksp) { 1059 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1060 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1061 } 1062 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 1063 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1064 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1065 1066 /* setup FETI-DP operators 1067 If fetidp->statechanged is true, we need to update the operators 1068 needed in the saddle-point case. This should be replaced 1069 by a better logic when the FETI-DP matrix and preconditioner will 1070 have their own classes */ 1071 if (pcbddc->new_primal_space || fetidp->statechanged) { 1072 Mat F; /* the FETI-DP matrix */ 1073 PC D; /* the FETI-DP preconditioner */ 1074 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 1075 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 1076 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 1077 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 1078 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 1079 ierr = PetscObjectIncrementTabLevel((PetscObject)D,(PetscObject)fetidp->innerksp,0);CHKERRQ(ierr); 1080 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 1081 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 1082 ierr = MatDestroy(&F);CHKERRQ(ierr); 1083 ierr = PCDestroy(&D);CHKERRQ(ierr); 1084 if (fetidp->check) { 1085 PetscViewer viewer; 1086 1087 if (!pcbddc->dbg_viewer) { 1088 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)); 1089 } else { 1090 viewer = pcbddc->dbg_viewer; 1091 } 1092 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr); 1093 } 1094 } 1095 fetidp->statechanged = PETSC_FALSE; 1096 pcbddc->new_primal_space = PETSC_FALSE; 1097 1098 /* propagate settings to the inner solve */ 1099 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 1100 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 1101 if (ksp->res_hist) { 1102 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 1103 } 1104 ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr); 1105 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 1106 PetscFunctionReturn(0); 1107 } 1108 1109 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 1110 { 1111 PetscErrorCode ierr; 1112 Mat F,A; 1113 MatNullSpace nsp; 1114 Vec X,B,Xl,Bl; 1115 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1116 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1117 1118 PetscFunctionBegin; 1119 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1120 if (fetidp->saddlepoint) { 1121 ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr); 1122 } 1123 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 1124 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 1125 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 1126 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 1127 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 1128 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 1129 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 1130 if (ksp->transpose_solve) { 1131 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1132 } else { 1133 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1134 } 1135 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 1136 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr); 1137 if (nsp) { 1138 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr); 1139 } 1140 /* update ksp with stats from inner ksp */ 1141 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 1142 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 1143 ksp->totalits += ksp->its; 1144 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 1145 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 1146 pcbddc->temp_solution_used = PETSC_FALSE; 1147 pcbddc->rhs_change = PETSC_FALSE; 1148 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1149 PetscFunctionReturn(0); 1150 } 1151 1152 static PetscErrorCode KSPReset_FETIDP(KSP ksp) 1153 { 1154 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1155 PC_BDDC *pcbddc; 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr); 1160 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr); 1161 /* avoid PCReset that does not take into account ref counting */ 1162 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1163 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1164 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1165 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1166 pcbddc->symmetric_primal = PETSC_FALSE; 1167 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1168 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1169 fetidp->saddlepoint = PETSC_FALSE; 1170 fetidp->matstate = -1; 1171 fetidp->matnnzstate = -1; 1172 fetidp->statechanged = PETSC_TRUE; 1173 PetscFunctionReturn(0); 1174 } 1175 1176 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 1177 { 1178 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1179 PetscErrorCode ierr; 1180 1181 PetscFunctionBegin; 1182 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr); 1183 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1184 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1185 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 1186 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 1187 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 1188 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 1189 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr); 1190 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 1194 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 1195 { 1196 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1197 PetscErrorCode ierr; 1198 PetscBool iascii; 1199 1200 PetscFunctionBegin; 1201 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1202 if (iascii) { 1203 ierr = PetscViewerASCIIPrintf(viewer," fully redundant: %d\n",fetidp->fully_redundant);CHKERRQ(ierr); 1204 ierr = PetscViewerASCIIPrintf(viewer," saddle point: %d\n",fetidp->saddlepoint);CHKERRQ(ierr); 1205 ierr = PetscViewerASCIIPrintf(viewer,"Inner KSP solver details\n");CHKERRQ(ierr); 1206 } 1207 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1208 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 1209 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1210 if (iascii) { 1211 ierr = PetscViewerASCIIPrintf(viewer,"Inner BDDC solver details\n");CHKERRQ(ierr); 1212 } 1213 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1214 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 1215 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1216 PetscFunctionReturn(0); 1217 } 1218 1219 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 1220 { 1221 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1222 PetscErrorCode ierr; 1223 1224 PetscFunctionBegin; 1225 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 1226 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1227 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 1228 if (!fetidp->userbddc) { 1229 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1230 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 1231 } 1232 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 1233 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 1234 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 1235 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr); 1236 ierr = PetscOptionsTail();CHKERRQ(ierr); 1237 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 /*MC 1242 KSPFETIDP - The FETI-DP method 1243 1244 This class implements the FETI-DP method [1]. 1245 The matrix for the KSP must be of type MATIS. 1246 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object. 1247 1248 Options Database Keys: 1249 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 1250 . -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 1251 . -ksp_fetidp_saddlepoint_flip <false> : usually, an incompressible Stokes problem is written as 1252 | A B^T | | v | = | f | 1253 | B 0 | | p | = | g | 1254 with B representing -\int_\Omega \nabla \cdot u q. 1255 If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as 1256 | A B^T | | v | = | f | 1257 |-B 0 | | p | = |-g | 1258 . -ksp_fetidp_pressure_field <-1> : activates support for saddle point problems, and identifies the pressure field id. 1259 If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals(). 1260 - -ksp_fetidp_pressure_all <false> : if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field. 1261 1262 Level: Advanced 1263 1264 Notes: 1265 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., 1266 .vb 1267 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 1268 .ve 1269 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 1270 1271 For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator(). 1272 Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix). 1273 If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver. 1274 Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g. 1275 .vb 1276 -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps 1277 .ve 1278 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 1279 non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g. 1280 .vb 1281 -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu 1282 .ve 1283 1284 Some of the basic options such as the maximum number of iterations and tolerances are automatically passed from this KSP to the inner KSP that actually performs the iterations. 1285 1286 The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution. 1287 1288 Developer Notes: 1289 Even though this method does not directly use any norms, the user is allowed to set the KSPNormType to any value. 1290 This is so users do not have to change KSPNormType options when they switch from other KSP methods to this one. 1291 1292 References: 1293 .vb 1294 . [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 1295 . [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 1296 .ve 1297 1298 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP() 1299 M*/ 1300 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp) 1301 { 1302 PetscErrorCode ierr; 1303 KSP_FETIDP *fetidp; 1304 KSP_FETIDPMon *monctx; 1305 PC_BDDC *pcbddc; 1306 PC pc; 1307 1308 PetscFunctionBegin; 1309 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);CHKERRQ(ierr); 1310 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);CHKERRQ(ierr); 1311 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1312 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1313 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1314 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1315 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr); 1316 1317 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr); 1318 fetidp->matstate = -1; 1319 fetidp->matnnzstate = -1; 1320 fetidp->statechanged = PETSC_TRUE; 1321 1322 ksp->data = (void*)fetidp; 1323 ksp->ops->setup = KSPSetUp_FETIDP; 1324 ksp->ops->solve = KSPSolve_FETIDP; 1325 ksp->ops->destroy = KSPDestroy_FETIDP; 1326 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP; 1327 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP; 1328 ksp->ops->view = KSPView_FETIDP; 1329 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP; 1330 ksp->ops->buildsolution = KSPBuildSolution_FETIDP; 1331 ksp->ops->buildresidual = KSPBuildResidualDefault; 1332 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 1333 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1334 /* create the inner KSP for the Lagrange multipliers */ 1335 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr); 1336 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr); 1337 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1338 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr); 1339 /* monitor */ 1340 ierr = PetscNew(&monctx);CHKERRQ(ierr); 1341 monctx->parentksp = ksp; 1342 fetidp->monctx = monctx; 1343 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr); 1344 /* create the inner BDDC */ 1345 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1346 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1347 /* make sure we always obtain a consistent FETI-DP matrix 1348 for symmetric problems, the user can always customize it through the command line */ 1349 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1350 pcbddc->symmetric_primal = PETSC_FALSE; 1351 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1352 /* composed functions */ 1353 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr); 1354 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr); 1356 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);CHKERRQ(ierr); 1357 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */ 1358 ksp->setupnewmatrix = PETSC_TRUE; 1359 PetscFunctionReturn(0); 1360 } 1361