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