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