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