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