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