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