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 if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing parallel pressure IS"); 730 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr); 731 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr); 732 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 733 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr); 734 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 735 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 736 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 737 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs+ni);CHKERRQ(ierr); 738 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr); 739 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 740 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs+ni,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 741 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 742 ierr = ISDestroy(&is1);CHKERRQ(ierr); 743 } else { 744 if (!lP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing sequential pressure IS"); 745 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr); 746 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 747 for (i=0;i<ni;i++) 748 if (idxs[i] >=0 && idxs[i] < n) 749 matis->sf_leafdata[idxs[i]] = 1; 750 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 751 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 752 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr); 753 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 754 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 755 ierr = ISDestroy(&is1);CHKERRQ(ierr); 756 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 757 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 758 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr); 759 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 760 } 761 ierr = PetscFree(widxs);CHKERRQ(ierr); 762 763 /* If there's any "interior pressure", 764 we may want to use a discrete harmonic solver instead 765 of a Stokes harmonic for the Dirichlet preconditioner 766 Need to extract the interior velocity dofs in interior dofs ordering (iV) 767 and interior pressure dofs in local ordering (iP) */ 768 if (!allp) { 769 ISLocalToGlobalMapping l2g_t; 770 771 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr); 772 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr); 773 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr); 774 ierr = ISDestroy(&is1);CHKERRQ(ierr); 775 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g_t);CHKERRQ(ierr); 776 ierr = ISGlobalToLocalMappingApplyIS(l2g_t,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr); 777 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr); 778 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr); 779 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iV",i,j); 780 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iV",(PetscObject)is1);CHKERRQ(ierr); 781 ierr = ISLocalToGlobalMappingDestroy(&l2g_t);CHKERRQ(ierr); 782 ierr = ISDestroy(&is1);CHKERRQ(ierr); 783 ierr = ISDestroy(&is2);CHKERRQ(ierr); 784 } 785 ierr = ISDestroy(&lPall);CHKERRQ(ierr); 786 ierr = ISDestroy(&II);CHKERRQ(ierr); 787 788 /* exclude selected pressures from the inner BDDC */ 789 if (pcbddc->DirichletBoundariesLocal) { 790 IS list[2],plP,isout; 791 PetscInt np; 792 793 /* need a parallel IS */ 794 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 795 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 796 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr); 797 list[0] = plP; 798 list[1] = pcbddc->DirichletBoundariesLocal; 799 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 800 ierr = ISDestroy(&plP);CHKERRQ(ierr); 801 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 802 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr); 803 ierr = ISDestroy(&isout);CHKERRQ(ierr); 804 } else if (pcbddc->DirichletBoundaries) { 805 IS list[2],isout; 806 807 list[0] = pP; 808 list[1] = pcbddc->DirichletBoundaries; 809 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 810 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr); 811 ierr = ISDestroy(&isout);CHKERRQ(ierr); 812 } else { 813 IS plP; 814 PetscInt np; 815 816 /* need a parallel IS */ 817 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 818 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 819 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr); 820 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr); 821 ierr = ISDestroy(&plP);CHKERRQ(ierr); 822 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 823 } 824 ierr = ISDestroy(&lP);CHKERRQ(ierr); 825 fetidp->pP = pP; 826 } 827 828 /* total number of selected pressure dofs */ 829 ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr); 830 831 /* Set operator for inner BDDC */ 832 if (totP || fetidp->rhs_flip) { 833 ierr = MatDuplicate(A,MAT_COPY_VALUES,&nA);CHKERRQ(ierr); 834 } else { 835 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 836 nA = A; 837 } 838 if (fetidp->rhs_flip) { 839 ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr); 840 if (totP) { 841 Mat lA2; 842 843 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr); 844 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr); 845 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr); 846 ierr = MatDestroy(&lA2);CHKERRQ(ierr); 847 } 848 } 849 850 if (totP) { 851 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 852 ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr); 853 } else { 854 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr); 855 } 856 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr); 857 ierr = MatDestroy(&nA);CHKERRQ(ierr); 858 859 /* non-zero rhs on interior dofs when applying the preconditioner */ 860 if (totP) pcbddc->switch_static = PETSC_TRUE; 861 862 /* if there are no pressures, set inner bddc flag for benign saddle point */ 863 if (!totP) { 864 pcbddc->benign_saddle_point = PETSC_TRUE; 865 pcbddc->compute_nonetflux = PETSC_TRUE; 866 } 867 868 /* Divergence mat */ 869 if (totP) { 870 Mat B; 871 IS P; 872 PetscBool save; 873 874 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 875 ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 876 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */ 877 ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,PETSC_FALSE,NULL);CHKERRQ(ierr); 878 pcbddc->compute_nonetflux = save; 879 ierr = MatDestroy(&B);CHKERRQ(ierr); 880 } 881 882 /* Operators for pressure preconditioner */ 883 if (totP) { 884 885 /* Extract pressure block */ 886 if (!pisz) { 887 Mat C; 888 889 ierr = MatCreateSubMatrix(A,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 890 ierr = MatScale(C,-1.);CHKERRQ(ierr); 891 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr); 892 ierr = MatDestroy(&C);CHKERRQ(ierr); 893 } else if (A != Ap) { /* user has provided a different Pmat, use it to extract the pressure preconditioner */ 894 Mat C; 895 896 ierr = MatCreateSubMatrix(Ap,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 897 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 898 ierr = MatDestroy(&C);CHKERRQ(ierr); 899 } 900 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 901 902 /* Preconditioned operator for the pressure block */ 903 if (PPmat) { 904 Mat C; 905 IS Pall; 906 PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg; 907 PetscBool ismatis; 908 909 ierr = PetscObjectTypeCompare((PetscObject)PPmat,MATIS,&ismatis);CHKERRQ(ierr); 910 if (ismatis) { 911 ierr = MatISGetMPIXAIJ(PPmat,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 912 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 913 ierr = MatDestroy(&C);CHKERRQ(ierr); 914 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 915 } 916 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr); 917 ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr); 918 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 919 ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr); 920 ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr); 921 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 922 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 923 ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr); 924 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 925 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 926 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 927 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); 928 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); 929 if (PAM == AM) { /* monolithic ordering, restrict to pressure */ 930 ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 931 } else if (pAg == PAM) { /* global ordering for pressure only */ 932 if (!allp) { /* solving for interface pressure only */ 933 IS restr; 934 935 ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr); 936 ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 937 ierr = ISDestroy(&restr);CHKERRQ(ierr); 938 } else { 939 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 940 C = PPmat; 941 } 942 } else if (pIg == PAM) { /* global ordering for selected pressure only */ 943 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 944 C = PPmat; 945 } else { 946 SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix"); 947 } 948 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 949 ierr = MatDestroy(&C);CHKERRQ(ierr); 950 } else { 951 Mat C; 952 953 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr); 954 if (C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */ 955 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)C);CHKERRQ(ierr); 956 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */ 957 PetscInt nl; 958 959 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr); 960 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&PPmat);CHKERRQ(ierr); 961 ierr = MatSetSizes(PPmat,nl,nl,totP,totP);CHKERRQ(ierr); 962 ierr = MatSetType(PPmat,MATAIJ);CHKERRQ(ierr); 963 ierr = MatMPIAIJSetPreallocation(PPmat,1,NULL,0,NULL);CHKERRQ(ierr); 964 ierr = MatSeqAIJSetPreallocation(PPmat,1,NULL);CHKERRQ(ierr); 965 ierr = MatAssemblyBegin(PPmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 966 ierr = MatAssemblyEnd(PPmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 967 ierr = MatShift(PPmat,1.);CHKERRQ(ierr); 968 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 969 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 970 } 971 } 972 } else { /* totP == 0 */ 973 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr); 974 } 975 } 976 PetscFunctionReturn(0); 977 } 978 979 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 980 { 981 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 982 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 983 PetscBool flg; 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 988 /* set up BDDC */ 989 ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr); 990 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 991 /* FETI-DP as it is implemented needs an exact coarse solver */ 992 if (pcbddc->coarse_ksp) { 993 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 994 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 995 } 996 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 997 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 998 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 999 1000 /* setup FETI-DP operators 1001 If fetidp->statechanged is true, we need update the operators 1002 that are needed in the saddle-point case. This should be replaced 1003 by a better logic when the FETI-DP matrix and preconditioner will 1004 have their own classes */ 1005 if (pcbddc->new_primal_space || fetidp->statechanged) { 1006 Mat F; /* the FETI-DP matrix */ 1007 PC D; /* the FETI-DP preconditioner */ 1008 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 1009 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 1010 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 1011 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 1012 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 1013 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 1014 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 1015 ierr = MatDestroy(&F);CHKERRQ(ierr); 1016 ierr = PCDestroy(&D);CHKERRQ(ierr); 1017 if (fetidp->check) { 1018 PetscViewer viewer; 1019 1020 if (!pcbddc->dbg_viewer) { 1021 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)); 1022 } else { 1023 viewer = pcbddc->dbg_viewer; 1024 } 1025 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr); 1026 } 1027 } 1028 fetidp->statechanged = PETSC_FALSE; 1029 pcbddc->new_primal_space = PETSC_FALSE; 1030 1031 /* propagate settings to the inner solve */ 1032 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 1033 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 1034 if (ksp->res_hist) { 1035 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 1036 } 1037 ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr); 1038 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 1043 { 1044 PetscErrorCode ierr; 1045 Mat F,A; 1046 MatNullSpace nsp; 1047 Vec X,B,Xl,Bl; 1048 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1049 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1050 1051 PetscFunctionBegin; 1052 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1053 if (fetidp->saddlepoint) { 1054 ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr); 1055 } 1056 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 1057 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 1058 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 1059 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 1060 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 1061 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 1062 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 1063 if (ksp->transpose_solve) { 1064 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1065 } else { 1066 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1067 } 1068 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 1069 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr); 1070 if (nsp) { 1071 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr); 1072 } 1073 /* update ksp with stats from inner ksp */ 1074 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 1075 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 1076 ksp->totalits += ksp->its; 1077 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 1078 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 1079 pcbddc->temp_solution_used = PETSC_FALSE; 1080 pcbddc->rhs_change = PETSC_FALSE; 1081 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1082 PetscFunctionReturn(0); 1083 } 1084 1085 static PetscErrorCode KSPReset_FETIDP(KSP ksp) 1086 { 1087 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1088 PC_BDDC *pcbddc; 1089 PetscErrorCode ierr; 1090 1091 PetscFunctionBegin; 1092 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr); 1093 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr); 1094 /* avoid PCReset that does not take into account ref counting */ 1095 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1096 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1097 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1098 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1099 pcbddc->symmetric_primal = PETSC_FALSE; 1100 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1101 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1102 fetidp->saddlepoint = PETSC_FALSE; 1103 fetidp->matstate = -1; 1104 fetidp->matnnzstate = -1; 1105 fetidp->statechanged = PETSC_TRUE; 1106 PetscFunctionReturn(0); 1107 } 1108 1109 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 1110 { 1111 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1112 PetscErrorCode ierr; 1113 1114 PetscFunctionBegin; 1115 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr); 1116 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1117 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1118 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 1119 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 1120 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 1121 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 1122 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr); 1123 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 1128 { 1129 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1130 PetscErrorCode ierr; 1131 PetscBool iascii; 1132 1133 PetscFunctionBegin; 1134 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1135 if (iascii) { 1136 ierr = PetscViewerASCIIPrintf(viewer," fully redundant: %D\n",fetidp->fully_redundant);CHKERRQ(ierr); 1137 ierr = PetscViewerASCIIPrintf(viewer," saddle point: %D\n",fetidp->saddlepoint);CHKERRQ(ierr); 1138 ierr = PetscViewerASCIIPrintf(viewer," inner solver details\n");CHKERRQ(ierr); 1139 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1140 } 1141 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 1142 if (iascii) { 1143 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1144 ierr = PetscViewerASCIIPrintf(viewer," BDDC solver details\n");CHKERRQ(ierr); 1145 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1146 } 1147 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 1148 if (iascii) { 1149 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1150 } 1151 PetscFunctionReturn(0); 1152 } 1153 1154 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 1155 { 1156 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 1161 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1162 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 1163 if (!fetidp->userbddc) { 1164 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1165 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 1166 } 1167 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 1168 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 1169 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 1170 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr); 1171 ierr = PetscOptionsTail();CHKERRQ(ierr); 1172 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 /*MC 1177 KSPFETIDP - The FETI-DP method 1178 1179 This class implements the FETI-DP method [1]. 1180 The matrix for the KSP must be of type MATIS. 1181 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object. 1182 1183 Options Database Keys: 1184 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 1185 . -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 1186 . -ksp_fetidp_saddlepoint_flip <false> : usually, an incompressible Stokes problem is written as 1187 | A B^T | | v | = | f | 1188 | B 0 | | p | = | g | 1189 with B representing -\int_\Omega \nabla \cdot u q. 1190 If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as 1191 | A B^T | | v | = | f | 1192 |-B 0 | | p | = |-g | 1193 . -ksp_fetidp_pressure_field <-1> : activates support for saddle point problems, and identifies the pressure field id. 1194 If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals(). 1195 . -ksp_fetidp_pressure_iszero <true> : if false, extracts the pressure block from the matrix (i.e. for Almost Incompressible Elasticity) 1196 - -ksp_fetidp_pressure_all <false> : if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field. 1197 1198 Level: Advanced 1199 1200 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., 1201 .vb 1202 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 1203 .ve 1204 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 1205 1206 For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator(). 1207 If it's not provided, an identity matrix will be created; the user then needs to scale it through a Richardson solver. 1208 Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g. 1209 .vb 1210 -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_package mumps 1211 .ve 1212 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 1213 non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g. 1214 .vb 1215 -fetidp_bddelta_pc_factor_mat_solver_package mumps -my_fetidp_bddelta_pc_type lu 1216 .ve 1217 1218 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. 1219 1220 The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution. 1221 1222 Developer Notes: Even though this method does not directly use any norms, the user is allowed to set the KSPNormType to any value. 1223 This is so users do not have to change KSPNormType options when they switch from other KSP methods to this one. 1224 1225 References: 1226 .vb 1227 . [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 1228 . [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 1229 .ve 1230 1231 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP() 1232 M*/ 1233 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp) 1234 { 1235 PetscErrorCode ierr; 1236 KSP_FETIDP *fetidp; 1237 KSP_FETIDPMon *monctx; 1238 PC_BDDC *pcbddc; 1239 PC pc; 1240 1241 PetscFunctionBegin; 1242 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);CHKERRQ(ierr); 1243 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);CHKERRQ(ierr); 1244 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1245 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1246 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1247 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1248 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr); 1249 1250 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr); 1251 fetidp->matstate = -1; 1252 fetidp->matnnzstate = -1; 1253 fetidp->statechanged = PETSC_TRUE; 1254 1255 ksp->data = (void*)fetidp; 1256 ksp->ops->setup = KSPSetUp_FETIDP; 1257 ksp->ops->solve = KSPSolve_FETIDP; 1258 ksp->ops->destroy = KSPDestroy_FETIDP; 1259 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP; 1260 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP; 1261 ksp->ops->view = KSPView_FETIDP; 1262 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP; 1263 ksp->ops->buildsolution = KSPBuildSolution_FETIDP; 1264 ksp->ops->buildresidual = KSPBuildResidualDefault; 1265 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 1266 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1267 /* create the inner KSP for the Lagrange multipliers */ 1268 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr); 1269 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr); 1270 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1271 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr); 1272 /* monitor */ 1273 ierr = PetscNew(&monctx);CHKERRQ(ierr); 1274 monctx->parentksp = ksp; 1275 fetidp->monctx = monctx; 1276 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr); 1277 /* create the inner BDDC */ 1278 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1279 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1280 /* make sure we always obtain a consistent FETI-DP matrix 1281 for symmetric problems, the user can always customize it through the command line */ 1282 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1283 pcbddc->symmetric_primal = PETSC_FALSE; 1284 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1285 /* composed functions */ 1286 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr); 1287 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr); 1288 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr); 1289 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);CHKERRQ(ierr); 1290 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */ 1291 ksp->setupnewmatrix = PETSC_TRUE; 1292 PetscFunctionReturn(0); 1293 } 1294