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: The operator can be either passed in a) monolithic global ordering, b) pressure-only global ordering 90 or c) interface pressure ordering (if -ksp_fetidp_pressure_all false). 91 In cases b) and c), the pressure ordering of dofs needs to satisfy 92 pid_1 < pid_2 iff gid_1 < gid_2 93 where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding 94 id in the monolithic global ordering. 95 96 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators 97 @*/ 98 PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P) 99 { 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 104 if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,2); 105 ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperator_C",(KSP,Mat),(ksp,P));CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp) 110 { 111 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 112 113 PetscFunctionBegin; 114 *innerksp = fetidp->innerksp; 115 PetscFunctionReturn(0); 116 } 117 118 /*@ 119 KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers 120 121 Input Parameters: 122 + ksp - the FETI-DP KSP 123 - innerksp - the KSP for the multipliers 124 125 Level: advanced 126 127 Notes: 128 129 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC 130 @*/ 131 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp) 132 { 133 PetscErrorCode ierr; 134 135 PetscFunctionBegin; 136 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 137 PetscValidPointer(innerksp,2); 138 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc) 143 { 144 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 145 146 PetscFunctionBegin; 147 *pc = fetidp->innerbddc; 148 PetscFunctionReturn(0); 149 } 150 151 /*@ 152 KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 153 154 Input Parameters: 155 + ksp - the FETI-DP Krylov solver 156 - pc - the BDDC preconditioner 157 158 Level: advanced 159 160 Notes: 161 162 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP 163 @*/ 164 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc) 165 { 166 PetscErrorCode ierr; 167 168 PetscFunctionBegin; 169 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 170 PetscValidPointer(pc,2); 171 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc) 176 { 177 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr); 182 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 183 fetidp->innerbddc = pc; 184 fetidp->userbddc = PETSC_TRUE; 185 PetscFunctionReturn(0); 186 } 187 188 /*@ 189 KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 190 191 Collective on KSP 192 193 Input Parameters: 194 + ksp - the FETI-DP Krylov solver 195 - pc - the BDDC preconditioner 196 197 Level: advanced 198 199 Notes: 200 201 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP 202 @*/ 203 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc) 204 { 205 PetscBool isbddc; 206 PetscErrorCode ierr; 207 208 PetscFunctionBegin; 209 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 210 PetscValidHeaderSpecific(pc,PC_CLASSID,2); 211 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr); 212 if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner"); 213 ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V) 218 { 219 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 220 Mat F; 221 Vec Xl; 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 226 ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr); 227 if (v) { 228 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr); 229 *V = v; 230 } else { 231 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr); 232 } 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx) 237 { 238 KSP_FETIDPMon *monctx = (KSP_FETIDPMon*)ctx; 239 PetscErrorCode ierr; 240 241 PetscFunctionBegin; 242 ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig) 247 { 248 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin) 257 { 258 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr); 263 PetscFunctionReturn(0); 264 } 265 266 static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer) 267 { 268 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 269 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 270 PC_IS *pcis = (PC_IS*)fetidp->innerbddc->data; 271 Mat_IS *matis = (Mat_IS*)fetidp->innerbddc->pmat->data; 272 Mat F; 273 FETIDPMat_ctx fetidpmat_ctx; 274 Vec test_vec,test_vec_p = NULL,fetidp_global; 275 IS dirdofs,isvert; 276 MPI_Comm comm = PetscObjectComm((PetscObject)ksp); 277 PetscScalar sval,*array; 278 PetscReal val,rval; 279 const PetscInt *vertex_indices; 280 PetscInt i,n_vertices; 281 PetscBool isascii; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 PetscCheckSameComm(ksp,1,viewer,2); 286 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 287 if (!isascii) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported viewer"); 288 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP MAT --------------\n");CHKERRQ(ierr); 289 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 290 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 291 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 292 ierr = MatView(F,viewer);CHKERRQ(ierr); 293 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 294 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 295 ierr = MatShellGetContext(F,(void**)&fetidpmat_ctx);CHKERRQ(ierr); 296 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP TESTS--------------\n");CHKERRQ(ierr); 297 ierr = PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");CHKERRQ(ierr); 298 ierr = PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");CHKERRQ(ierr); 299 if (fetidp->fully_redundant) { 300 ierr = PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");CHKERRQ(ierr); 301 } else { 302 ierr = PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");CHKERRQ(ierr); 303 } 304 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 305 306 /* Get Vertices used to define the BDDC */ 307 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr); 308 ierr = ISGetLocalSize(isvert,&n_vertices);CHKERRQ(ierr); 309 ierr = ISGetIndices(isvert,&vertex_indices);CHKERRQ(ierr); 310 311 /******************************************************************/ 312 /* TEST A/B: Test numbering of global fetidp dofs */ 313 /******************************************************************/ 314 ierr = MatCreateVecs(F,&fetidp_global,NULL);CHKERRQ(ierr); 315 ierr = VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);CHKERRQ(ierr); 316 ierr = VecSet(fetidp_global,1.0);CHKERRQ(ierr); 317 ierr = VecSet(test_vec,1.);CHKERRQ(ierr); 318 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 319 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 320 if (fetidpmat_ctx->l2g_p) { 321 ierr = VecDuplicate(fetidpmat_ctx->vP,&test_vec_p);CHKERRQ(ierr); 322 ierr = VecSet(test_vec_p,1.);CHKERRQ(ierr); 323 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 324 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 325 } 326 ierr = VecAXPY(test_vec,-1.0,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 327 ierr = VecNorm(test_vec,NORM_INFINITY,&val);CHKERRQ(ierr); 328 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 329 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 330 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc: % 1.14e\n",rval);CHKERRQ(ierr); 331 332 if (fetidpmat_ctx->l2g_p) { 333 ierr = VecAXPY(test_vec_p,-1.0,fetidpmat_ctx->vP);CHKERRQ(ierr); 334 ierr = VecNorm(test_vec_p,NORM_INFINITY,&val);CHKERRQ(ierr); 335 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc (p): % 1.14e\n",rval);CHKERRQ(ierr); 337 } 338 339 if (fetidp->fully_redundant) { 340 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 341 ierr = VecSet(fetidpmat_ctx->lambda_local,0.5);CHKERRQ(ierr); 342 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 343 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 344 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr); 345 val = PetscRealPart(sval)-fetidpmat_ctx->n_lambda; 346 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 347 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob: % 1.14e\n",rval);CHKERRQ(ierr); 348 } 349 350 if (fetidpmat_ctx->l2g_p) { 351 ierr = VecSet(pcis->vec1_N,1.0);CHKERRQ(ierr); 352 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 353 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 354 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 355 356 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 357 ierr = VecSet(fetidpmat_ctx->vP,-1.0);CHKERRQ(ierr); 358 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 359 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 360 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 361 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 362 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 363 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 364 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr); 365 val = PetscRealPart(sval); 366 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 367 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob (p): % 1.14e\n",rval);CHKERRQ(ierr); 368 } 369 370 /******************************************************************/ 371 /* TEST C: It should hold B_delta*w=0, w\in\widehat{W} */ 372 /* This is the meaning of the B matrix */ 373 /******************************************************************/ 374 375 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 376 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 377 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 378 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 379 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 380 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 381 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 382 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 383 /* Action of B_delta */ 384 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 385 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 386 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 387 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 388 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr); 389 ierr = PetscViewerASCIIPrintf(viewer,"C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",val);CHKERRQ(ierr); 390 391 /******************************************************************/ 392 /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W} */ 393 /* E_D = R_D^TR */ 394 /* P_D = B_{D,delta}^T B_{delta} */ 395 /* eq.44 Mandel Tezaur and Dohrmann 2005 */ 396 /******************************************************************/ 397 398 /* compute a random vector in \widetilde{W} */ 399 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 400 /* set zero at vertices and essential dofs */ 401 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 402 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0; 403 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirdofs);CHKERRQ(ierr); 404 if (dirdofs) { 405 const PetscInt *idxs; 406 PetscInt ndir; 407 408 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr); 409 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr); 410 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0; 411 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr); 412 } 413 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 414 /* store w for final comparison */ 415 ierr = VecDuplicate(pcis->vec1_B,&test_vec);CHKERRQ(ierr); 416 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 417 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 418 419 /* Jump operator P_D : results stored in pcis->vec1_B */ 420 /* Action of B_delta */ 421 ierr = MatMult(fetidpmat_ctx->B_delta,test_vec,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 422 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 423 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 424 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 425 /* Action of B_Ddelta^T */ 426 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 427 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 428 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 429 430 /* Average operator E_D : results stored in pcis->vec2_B */ 431 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,test_vec,pcis->vec1_global);CHKERRQ(ierr); 432 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 433 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 434 435 /* test E_D=I-P_D */ 436 ierr = VecAXPY(pcis->vec1_B,1.0,pcis->vec2_B);CHKERRQ(ierr); 437 ierr = VecAXPY(pcis->vec1_B,-1.0,test_vec);CHKERRQ(ierr); 438 ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&val);CHKERRQ(ierr); 439 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 440 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 441 ierr = PetscViewerASCIIPrintf(viewer,"D: CHECK infty norm of E_D + P_D - I: % 1.14e\n",PetscGlobalRank,val);CHKERRQ(ierr); 442 443 /******************************************************************/ 444 /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W} */ 445 /* eq.48 Mandel Tezaur and Dohrmann 2005 */ 446 /******************************************************************/ 447 448 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 449 /* set zero at vertices and essential dofs */ 450 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 451 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0; 452 if (dirdofs) { 453 const PetscInt *idxs; 454 PetscInt ndir; 455 456 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr); 457 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr); 458 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0; 459 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr); 460 } 461 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 462 463 /* Jump operator P_D : results stored in pcis->vec1_B */ 464 465 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 466 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 467 /* Action of B_delta */ 468 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 469 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 470 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 471 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 472 /* Action of B_Ddelta^T */ 473 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 474 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 475 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 476 /* scaling */ 477 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr); 478 ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&val);CHKERRQ(ierr); 479 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of R^T_D P_D: % 1.14e\n",val);CHKERRQ(ierr); 480 481 if (!fetidp->fully_redundant) { 482 /******************************************************************/ 483 /* TEST F: It should holds B_{delta}B^T_{D,delta}=I */ 484 /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005 */ 485 /******************************************************************/ 486 ierr = VecDuplicate(fetidp_global,&test_vec);CHKERRQ(ierr); 487 ierr = VecSetRandom(fetidp_global,NULL);CHKERRQ(ierr); 488 if (fetidpmat_ctx->l2g_p) { 489 ierr = VecSet(fetidpmat_ctx->vP,0.);CHKERRQ(ierr); 490 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 491 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 492 } 493 /* Action of B_Ddelta^T */ 494 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 495 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 496 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 497 /* Action of B_delta */ 498 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 499 ierr = VecSet(test_vec,0.0);CHKERRQ(ierr); 500 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 501 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 502 ierr = VecAXPY(fetidp_global,-1.,test_vec);CHKERRQ(ierr); 503 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr); 504 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of P^T_D - I: % 1.14e\n",val);CHKERRQ(ierr); 505 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 506 } 507 ierr = PetscViewerASCIIPrintf(viewer,"-------------------------------------\n");CHKERRQ(ierr); 508 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 509 ierr = VecDestroy(&test_vec_p);CHKERRQ(ierr); 510 ierr = ISDestroy(&dirdofs);CHKERRQ(ierr); 511 ierr = VecDestroy(&fetidp_global);CHKERRQ(ierr); 512 ierr = ISRestoreIndices(isvert,&vertex_indices);CHKERRQ(ierr); 513 ierr = PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr); 514 PetscFunctionReturn(0); 515 } 516 517 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp) 518 { 519 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 520 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 521 Mat A,Ap; 522 PetscInt fid = -1; 523 PetscMPIInt size; 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 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)ksp),&size);CHKERRQ(ierr); 546 fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint); 547 if (size == 1) fetidp->saddlepoint = PETSC_FALSE; 548 549 ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr); 550 ierr = PetscObjectTypeCompare((PetscObject)A,MATIS,&ismatis);CHKERRQ(ierr); 551 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Amat should be of type MATIS"); 552 553 /* Quiet return if the matrix states are unchanged. 554 Needed only for the saddle point case since it uses MatZeroRows 555 on a matrix that may not have changed */ 556 ierr = PetscObjectStateGet((PetscObject)A,&matstate);CHKERRQ(ierr); 557 ierr = MatGetNonzeroState(A,&matnnzstate);CHKERRQ(ierr); 558 if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(0); 559 fetidp->matstate = matstate; 560 fetidp->matnnzstate = matnnzstate; 561 fetidp->statechanged = fetidp->saddlepoint; 562 563 /* see if we have some fields attached */ 564 if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) { 565 DM dm; 566 PetscContainer c; 567 568 ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr); 569 ierr = PetscObjectQuery((PetscObject)A,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr); 570 if (dm) { 571 IS *fields; 572 PetscInt nf,i; 573 574 ierr = DMCreateFieldDecomposition(dm,&nf,NULL,&fields,NULL);CHKERRQ(ierr); 575 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,nf,fields);CHKERRQ(ierr); 576 for (i=0;i<nf;i++) { 577 ierr = ISDestroy(&fields[i]);CHKERRQ(ierr); 578 } 579 ierr = PetscFree(fields);CHKERRQ(ierr); 580 } else if (c) { 581 MatISLocalFields lf; 582 ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr); 583 ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr); 584 } 585 } 586 587 if (!fetidp->saddlepoint) { 588 ierr = PCSetOperators(fetidp->innerbddc,A,A);CHKERRQ(ierr); 589 } else { 590 Mat nA,lA; 591 Mat PPmat; 592 IS pP; 593 PetscInt totP; 594 595 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 596 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr); 597 598 pP = fetidp->pP; 599 if (!pP) { /* first time, need to compute pressure dofs */ 600 PC_IS *pcis = (PC_IS*)fetidp->innerbddc->data; 601 Mat_IS *matis = (Mat_IS*)(A->data); 602 ISLocalToGlobalMapping l2g; 603 IS lP = NULL,II,pII,lPall,Pall,is1,is2; 604 const PetscInt *idxs; 605 PetscInt nl,ni,*widxs; 606 PetscInt i,j,n_neigh,*neigh,*n_shared,**shared,*count; 607 PetscInt rst,ren,n; 608 PetscBool ploc; 609 610 ierr = MatGetLocalSize(A,&nl,NULL);CHKERRQ(ierr); 611 ierr = MatGetOwnershipRange(A,&rst,&ren);CHKERRQ(ierr); 612 ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr); 613 ierr = MatGetLocalToGlobalMapping(A,&l2g,NULL);CHKERRQ(ierr); 614 615 if (!pcis->is_I_local) { /* need to compute interior dofs */ 616 ierr = PetscCalloc1(n,&count);CHKERRQ(ierr); 617 ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 618 for (i=1;i<n_neigh;i++) 619 for (j=0;j<n_shared[i];j++) 620 count[shared[i][j]] += 1; 621 for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i; 622 ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 623 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr); 624 } else { 625 ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr); 626 II = pcis->is_I_local; 627 } 628 629 /* interior dofs in layout */ 630 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 631 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 632 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 633 ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr); 634 ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr); 635 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 636 ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr); 637 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 638 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 639 ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr); 640 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 641 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr); 642 643 /* pressure dofs */ 644 Pall = NULL; 645 lPall = NULL; 646 ploc = PETSC_FALSE; 647 if (fid >= 0) { 648 if (pcbddc->n_ISForDofsLocal) { 649 PetscInt np; 650 651 if (fid >= pcbddc->n_ISForDofsLocal) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofsLocal); 652 /* need a sequential IS */ 653 ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr); 654 ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 655 ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 656 ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 657 ploc = PETSC_TRUE; 658 } else if (pcbddc->n_ISForDofs) { 659 if (fid >= pcbddc->n_ISForDofs) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofs); 660 ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr); 661 Pall = pcbddc->ISForDofs[fid]; 662 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Missing fields! Use PCBDDCSetDofsSplitting/Local"); 663 } else { /* fallback to zero pressure block */ 664 IS list[2]; 665 666 ierr = MatFindZeroDiagonals(A,&list[1]);CHKERRQ(ierr); 667 ierr = ISComplement(list[1],rst,ren,&list[0]);CHKERRQ(ierr); 668 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,2,list);CHKERRQ(ierr); 669 ierr = ISDestroy(&list[0]);CHKERRQ(ierr); 670 Pall = list[1]; 671 } 672 /* if the user requested the entire pressure, 673 remove the interior pressure dofs from II (or pII) */ 674 if (allp) { 675 if (ploc) { 676 IS nII; 677 ierr = ISDifference(II,lPall,&nII);CHKERRQ(ierr); 678 ierr = ISDestroy(&II);CHKERRQ(ierr); 679 II = nII; 680 } else { 681 IS nII; 682 ierr = ISDifference(pII,Pall,&nII);CHKERRQ(ierr); 683 ierr = ISDestroy(&pII);CHKERRQ(ierr); 684 pII = nII; 685 } 686 } 687 if (ploc) { 688 ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr); 689 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 690 } else { 691 ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr); 692 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 693 /* need all local pressure dofs */ 694 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 695 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 696 ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr); 697 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 698 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 699 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 700 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 701 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 702 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 703 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 704 } 705 706 if (!Pall) { 707 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 708 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 709 ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr); 710 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr); 711 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 712 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr); 713 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 714 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 715 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 716 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr); 717 } 718 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject)Pall);CHKERRQ(ierr); 719 720 if (flip) { 721 PetscInt npl; 722 ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr); 723 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 724 ierr = MatCreateVecs(A,NULL,&fetidp->rhs_flip);CHKERRQ(ierr); 725 ierr = VecSet(fetidp->rhs_flip,1.);CHKERRQ(ierr); 726 ierr = VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 727 for (i=0;i<npl;i++) { 728 ierr = VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr); 729 } 730 ierr = VecAssemblyBegin(fetidp->rhs_flip);CHKERRQ(ierr); 731 ierr = VecAssemblyEnd(fetidp->rhs_flip);CHKERRQ(ierr); 732 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);CHKERRQ(ierr); 733 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 734 } 735 ierr = ISDestroy(&Pall);CHKERRQ(ierr); 736 ierr = ISDestroy(&pII);CHKERRQ(ierr); 737 738 /* local selected pressures in subdomain-wise and global ordering */ 739 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 740 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 741 if (!ploc) { 742 PetscInt *widxs2; 743 744 if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing parallel pressure IS"); 745 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr); 746 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr); 747 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 748 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr); 749 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 750 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 751 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 752 ierr = PetscMalloc1(ni,&widxs2);CHKERRQ(ierr); 753 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs2);CHKERRQ(ierr); 754 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr); 755 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 756 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs2,PETSC_OWN_POINTER,&is1);CHKERRQ(ierr); 757 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 758 ierr = ISDestroy(&is1);CHKERRQ(ierr); 759 } else { 760 if (!lP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing sequential pressure IS"); 761 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr); 762 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 763 for (i=0;i<ni;i++) 764 if (idxs[i] >=0 && idxs[i] < n) 765 matis->sf_leafdata[idxs[i]] = 1; 766 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 767 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 768 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr); 769 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 770 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 771 ierr = ISDestroy(&is1);CHKERRQ(ierr); 772 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 773 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 774 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr); 775 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 776 } 777 ierr = PetscFree(widxs);CHKERRQ(ierr); 778 779 /* If there's any "interior pressure", 780 we may want to use a discrete harmonic solver instead 781 of a Stokes harmonic for the Dirichlet preconditioner 782 Need to extract the interior velocity dofs in interior dofs ordering (iV) 783 and interior pressure dofs in local ordering (iP) */ 784 if (!allp) { 785 ISLocalToGlobalMapping l2g_t; 786 787 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr); 788 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr); 789 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr); 790 ierr = ISDestroy(&is1);CHKERRQ(ierr); 791 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g_t);CHKERRQ(ierr); 792 ierr = ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr); 793 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr); 794 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr); 795 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iV",i,j); 796 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);CHKERRQ(ierr); 797 ierr = ISLocalToGlobalMappingDestroy(&l2g_t);CHKERRQ(ierr); 798 ierr = ISDestroy(&is1);CHKERRQ(ierr); 799 ierr = ISDestroy(&is2);CHKERRQ(ierr); 800 } 801 ierr = ISDestroy(&lPall);CHKERRQ(ierr); 802 ierr = ISDestroy(&II);CHKERRQ(ierr); 803 804 /* exclude selected pressures from the inner BDDC */ 805 if (pcbddc->DirichletBoundariesLocal) { 806 IS list[2],plP,isout; 807 PetscInt np; 808 809 /* need a parallel IS */ 810 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 811 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 812 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr); 813 list[0] = plP; 814 list[1] = pcbddc->DirichletBoundariesLocal; 815 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 816 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr); 817 ierr = ISDestroy(&plP);CHKERRQ(ierr); 818 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 819 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr); 820 ierr = ISDestroy(&isout);CHKERRQ(ierr); 821 } else if (pcbddc->DirichletBoundaries) { 822 IS list[2],isout; 823 824 list[0] = pP; 825 list[1] = pcbddc->DirichletBoundaries; 826 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 827 ierr = ISSortRemoveDups(isout);CHKERRQ(ierr); 828 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr); 829 ierr = ISDestroy(&isout);CHKERRQ(ierr); 830 } else { 831 IS plP; 832 PetscInt np; 833 834 /* need a parallel IS */ 835 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 836 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 837 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr); 838 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr); 839 ierr = ISDestroy(&plP);CHKERRQ(ierr); 840 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 841 } 842 ierr = ISDestroy(&lP);CHKERRQ(ierr); 843 fetidp->pP = pP; 844 } 845 846 /* total number of selected pressure dofs */ 847 ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr); 848 849 /* Set operator for inner BDDC */ 850 if (totP || fetidp->rhs_flip) { 851 ierr = MatDuplicate(A,MAT_COPY_VALUES,&nA);CHKERRQ(ierr); 852 } else { 853 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 854 nA = A; 855 } 856 if (fetidp->rhs_flip) { 857 ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr); 858 if (totP) { 859 Mat lA2; 860 861 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr); 862 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr); 863 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr); 864 ierr = MatDestroy(&lA2);CHKERRQ(ierr); 865 } 866 } 867 868 if (totP) { 869 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 870 ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr); 871 } else { 872 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr); 873 } 874 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr); 875 ierr = MatDestroy(&nA);CHKERRQ(ierr); 876 877 /* non-zero rhs on interior dofs when applying the preconditioner */ 878 if (totP) pcbddc->switch_static = PETSC_TRUE; 879 880 /* if there are no interface pressures, set inner bddc flag for benign saddle point */ 881 if (!totP) { 882 pcbddc->benign_saddle_point = PETSC_TRUE; 883 pcbddc->compute_nonetflux = PETSC_TRUE; 884 } 885 886 /* Divergence mat */ 887 if (totP) { 888 Mat B; 889 IS P; 890 PetscBool save; 891 892 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 893 ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 894 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */ 895 ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,PETSC_FALSE,NULL);CHKERRQ(ierr); 896 pcbddc->compute_nonetflux = save; 897 ierr = MatDestroy(&B);CHKERRQ(ierr); 898 } 899 900 /* Operators for pressure preconditioner */ 901 if (totP) { 902 /* Extract pressure block if needed */ 903 if (!pisz) { 904 Mat C; 905 IS nzrows = NULL; 906 907 ierr = MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 908 ierr = MatFindNonzeroRows(C,&nzrows);CHKERRQ(ierr); 909 if (nzrows) { 910 PetscInt i; 911 912 ierr = ISGetSize(nzrows,&i);CHKERRQ(ierr); 913 ierr = ISDestroy(&nzrows);CHKERRQ(ierr); 914 if (!i) pisz = PETSC_TRUE; 915 } 916 if (!pisz) { 917 ierr = MatScale(C,-1.);CHKERRQ(ierr); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized */ 918 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr); 919 } 920 ierr = MatDestroy(&C);CHKERRQ(ierr); 921 } 922 if (A != Ap) { /* user has provided a different Pmat, use it to extract the pressure preconditioner */ 923 Mat C; 924 925 ierr = MatCreateSubMatrix(Ap,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 926 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 927 ierr = MatDestroy(&C);CHKERRQ(ierr); 928 } 929 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 930 931 /* Preconditioned operator for the pressure block */ 932 if (PPmat) { 933 Mat C; 934 IS Pall; 935 PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg; 936 PetscBool ismatis; 937 938 ierr = PetscObjectTypeCompare((PetscObject)PPmat,MATIS,&ismatis);CHKERRQ(ierr); 939 if (ismatis) { 940 ierr = MatISGetMPIXAIJ(PPmat,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 941 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 942 ierr = MatDestroy(&C);CHKERRQ(ierr); 943 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 944 } 945 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr); 946 ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr); 947 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 948 ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr); 949 ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr); 950 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 951 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 952 ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr); 953 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 954 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 955 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 956 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); 957 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); 958 if (PAM == AM) { /* monolithic ordering, restrict to pressure */ 959 ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 960 } else if (pAg == PAM) { /* global ordering for pressure only */ 961 if (!allp) { /* solving for interface pressure only */ 962 IS restr; 963 964 ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr); 965 ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 966 ierr = ISDestroy(&restr);CHKERRQ(ierr); 967 } else { 968 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 969 C = PPmat; 970 } 971 } else if (pIg == PAM) { /* global ordering for selected pressure only */ 972 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 973 C = PPmat; 974 } else { 975 SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix"); 976 } 977 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 978 ierr = MatDestroy(&C);CHKERRQ(ierr); 979 } else { 980 Mat C; 981 982 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr); 983 if (C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */ 984 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 985 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */ 986 PetscInt nl; 987 988 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr); 989 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&PPmat);CHKERRQ(ierr); 990 ierr = MatSetSizes(PPmat,nl,nl,totP,totP);CHKERRQ(ierr); 991 ierr = MatSetType(PPmat,MATAIJ);CHKERRQ(ierr); 992 ierr = MatMPIAIJSetPreallocation(PPmat,1,NULL,0,NULL);CHKERRQ(ierr); 993 ierr = MatSeqAIJSetPreallocation(PPmat,1,NULL);CHKERRQ(ierr); 994 ierr = MatAssemblyBegin(PPmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 995 ierr = MatAssemblyEnd(PPmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 996 ierr = MatShift(PPmat,1.);CHKERRQ(ierr); 997 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 998 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 999 } 1000 } 1001 } else { /* totP == 0 */ 1002 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr); 1003 } 1004 } 1005 PetscFunctionReturn(0); 1006 } 1007 1008 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 1009 { 1010 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1011 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1012 PetscBool flg; 1013 PetscErrorCode ierr; 1014 1015 PetscFunctionBegin; 1016 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 1017 /* set up BDDC */ 1018 ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr); 1019 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 1020 /* FETI-DP as it is implemented needs an exact coarse solver */ 1021 if (pcbddc->coarse_ksp) { 1022 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1023 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1024 } 1025 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 1026 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1027 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1028 1029 /* setup FETI-DP operators 1030 If fetidp->statechanged is true, we need update the operators 1031 that are needed in the saddle-point case. This should be replaced 1032 by a better logic when the FETI-DP matrix and preconditioner will 1033 have their own classes */ 1034 if (pcbddc->new_primal_space || fetidp->statechanged) { 1035 Mat F; /* the FETI-DP matrix */ 1036 PC D; /* the FETI-DP preconditioner */ 1037 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 1038 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 1039 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 1040 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 1041 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 1042 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 1043 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 1044 ierr = MatDestroy(&F);CHKERRQ(ierr); 1045 ierr = PCDestroy(&D);CHKERRQ(ierr); 1046 if (fetidp->check) { 1047 PetscViewer viewer; 1048 1049 if (!pcbddc->dbg_viewer) { 1050 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)); 1051 } else { 1052 viewer = pcbddc->dbg_viewer; 1053 } 1054 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr); 1055 } 1056 } 1057 fetidp->statechanged = PETSC_FALSE; 1058 pcbddc->new_primal_space = PETSC_FALSE; 1059 1060 /* propagate settings to the inner solve */ 1061 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 1062 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 1063 if (ksp->res_hist) { 1064 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 1065 } 1066 ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr); 1067 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 1072 { 1073 PetscErrorCode ierr; 1074 Mat F,A; 1075 MatNullSpace nsp; 1076 Vec X,B,Xl,Bl; 1077 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1078 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1079 1080 PetscFunctionBegin; 1081 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1082 if (fetidp->saddlepoint) { 1083 ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr); 1084 } 1085 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 1086 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 1087 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 1088 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 1089 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 1090 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 1091 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 1092 if (ksp->transpose_solve) { 1093 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1094 } else { 1095 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1096 } 1097 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 1098 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr); 1099 if (nsp) { 1100 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr); 1101 } 1102 /* update ksp with stats from inner ksp */ 1103 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 1104 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 1105 ksp->totalits += ksp->its; 1106 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 1107 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 1108 pcbddc->temp_solution_used = PETSC_FALSE; 1109 pcbddc->rhs_change = PETSC_FALSE; 1110 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1111 PetscFunctionReturn(0); 1112 } 1113 1114 static PetscErrorCode KSPReset_FETIDP(KSP ksp) 1115 { 1116 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1117 PC_BDDC *pcbddc; 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr); 1122 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr); 1123 /* avoid PCReset that does not take into account ref counting */ 1124 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1125 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1126 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1127 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1128 pcbddc->symmetric_primal = PETSC_FALSE; 1129 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1130 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1131 fetidp->saddlepoint = PETSC_FALSE; 1132 fetidp->matstate = -1; 1133 fetidp->matnnzstate = -1; 1134 fetidp->statechanged = PETSC_TRUE; 1135 PetscFunctionReturn(0); 1136 } 1137 1138 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 1139 { 1140 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr); 1145 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1146 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1147 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 1148 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 1149 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 1150 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 1151 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr); 1152 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 1157 { 1158 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1159 PetscErrorCode ierr; 1160 PetscBool iascii; 1161 1162 PetscFunctionBegin; 1163 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1164 if (iascii) { 1165 ierr = PetscViewerASCIIPrintf(viewer," fully redundant: %d\n",fetidp->fully_redundant);CHKERRQ(ierr); 1166 ierr = PetscViewerASCIIPrintf(viewer," saddle point: %d\n",fetidp->saddlepoint);CHKERRQ(ierr); 1167 ierr = PetscViewerASCIIPrintf(viewer," inner solver details\n");CHKERRQ(ierr); 1168 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1169 } 1170 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 1171 if (iascii) { 1172 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1173 ierr = PetscViewerASCIIPrintf(viewer," BDDC solver details\n");CHKERRQ(ierr); 1174 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1175 } 1176 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 1177 if (iascii) { 1178 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1179 } 1180 PetscFunctionReturn(0); 1181 } 1182 1183 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 1184 { 1185 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1186 PetscErrorCode ierr; 1187 1188 PetscFunctionBegin; 1189 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 1190 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1191 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 1192 if (!fetidp->userbddc) { 1193 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1194 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 1195 } 1196 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 1197 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 1198 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 1199 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr); 1200 ierr = PetscOptionsTail();CHKERRQ(ierr); 1201 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 1202 PetscFunctionReturn(0); 1203 } 1204 1205 /*MC 1206 KSPFETIDP - The FETI-DP method 1207 1208 This class implements the FETI-DP method [1]. 1209 The matrix for the KSP must be of type MATIS. 1210 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object. 1211 1212 Options Database Keys: 1213 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 1214 . -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 1215 . -ksp_fetidp_saddlepoint_flip <false> : usually, an incompressible Stokes problem is written as 1216 | A B^T | | v | = | f | 1217 | B 0 | | p | = | g | 1218 with B representing -\int_\Omega \nabla \cdot u q. 1219 If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as 1220 | A B^T | | v | = | f | 1221 |-B 0 | | p | = |-g | 1222 . -ksp_fetidp_pressure_field <-1> : activates support for saddle point problems, and identifies the pressure field id. 1223 If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals(). 1224 - -ksp_fetidp_pressure_all <false> : if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field. 1225 1226 Level: Advanced 1227 1228 Notes: 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., 1229 .vb 1230 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 1231 .ve 1232 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 1233 1234 For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator(). 1235 Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix). 1236 If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver. 1237 Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g. 1238 .vb 1239 -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps 1240 .ve 1241 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 1242 non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g. 1243 .vb 1244 -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu 1245 .ve 1246 1247 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. 1248 1249 The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution. 1250 1251 Developer Notes: 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