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