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