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