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 PetscBool save; 978 979 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 980 if (!pisz) { 981 IS F,V; 982 PetscInt m,M; 983 984 ierr = MatGetOwnershipRange(A,&m,&M);CHKERRQ(ierr); 985 ierr = ISCreateStride(PetscObjectComm((PetscObject)P),M-m,m,1,&F);CHKERRQ(ierr); 986 ierr = ISComplement(P,m,M,&V);CHKERRQ(ierr); 987 ierr = MatCreateSubMatrix(A,F,V,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 988 ierr = MatZeroRowsIS(B,P,0.0,NULL,NULL);CHKERRQ(ierr); 989 ierr = ISDestroy(&V);CHKERRQ(ierr); 990 ierr = ISDestroy(&F);CHKERRQ(ierr); 991 } else { 992 ierr = MatCreateSubMatrix(A,P,NULL,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 993 } 994 save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */ 995 ierr = PCBDDCSetDivergenceMat(fetidp->innerbddc,B,(PetscBool)!pisz,NULL);CHKERRQ(ierr); 996 pcbddc->compute_nonetflux = save; 997 ierr = MatDestroy(&B);CHKERRQ(ierr); 998 } 999 if (A != Ap) { /* user has provided a different Pmat, this always superseeds the setter (TODO: is it OK?) */ 1000 /* use monolithic operator, we restrict later */ 1001 ierr = KSPFETIDPSetPressureOperator(ksp,Ap);CHKERRQ(ierr); 1002 } 1003 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 1004 1005 /* PPmat not present, use some default choice */ 1006 if (!PPmat) { 1007 Mat C; 1008 1009 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject*)&C);CHKERRQ(ierr); 1010 if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */ 1011 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 1012 } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */ 1013 IS P; 1014 1015 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&P);CHKERRQ(ierr); 1016 ierr = MatCreateSubMatrix(A,P,P,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1017 ierr = MatScale(C,-1.);CHKERRQ(ierr); 1018 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 1019 ierr = MatDestroy(&C);CHKERRQ(ierr); 1020 } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */ 1021 PetscInt nl; 1022 1023 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr); 1024 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&C);CHKERRQ(ierr); 1025 ierr = MatSetSizes(C,nl,nl,totP,totP);CHKERRQ(ierr); 1026 ierr = MatSetType(C,MATAIJ);CHKERRQ(ierr); 1027 ierr = MatMPIAIJSetPreallocation(C,1,NULL,0,NULL);CHKERRQ(ierr); 1028 ierr = MatSeqAIJSetPreallocation(C,1,NULL);CHKERRQ(ierr); 1029 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1030 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1031 ierr = MatShift(C,1.);CHKERRQ(ierr); 1032 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 1033 ierr = MatDestroy(&C);CHKERRQ(ierr); 1034 } 1035 } 1036 1037 /* Preconditioned operator for the pressure block */ 1038 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 1039 if (PPmat) { 1040 Mat C; 1041 IS Pall; 1042 PetscInt AM,PAM,PAN,pam,pan,am,an,pl,pIl,pAg,pIg; 1043 1044 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_aP",(PetscObject*)&Pall);CHKERRQ(ierr); 1045 ierr = MatGetSize(A,&AM,NULL);CHKERRQ(ierr); 1046 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 1047 ierr = ISGetSize(Pall,&pAg);CHKERRQ(ierr); 1048 ierr = ISGetSize(fetidp->pP,&pIg);CHKERRQ(ierr); 1049 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 1050 ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr); 1051 ierr = ISGetLocalSize(Pall,&pIl);CHKERRQ(ierr); 1052 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 1053 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 1054 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 1055 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); 1056 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); 1057 if (PAM == AM) { /* monolithic ordering, restrict to pressure */ 1058 if (schp) { 1059 ierr = MatCreateSubMatrix(PPmat,Pall,Pall,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1060 } else { 1061 ierr = MatCreateSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1062 } 1063 } else if (pAg == PAM) { /* global ordering for pressure only */ 1064 if (!allp && !schp) { /* solving for interface pressure only */ 1065 IS restr; 1066 1067 ierr = ISRenumber(fetidp->pP,NULL,NULL,&restr);CHKERRQ(ierr); 1068 ierr = MatCreateSubMatrix(PPmat,restr,restr,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1069 ierr = ISDestroy(&restr);CHKERRQ(ierr); 1070 } else { 1071 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 1072 C = PPmat; 1073 } 1074 } else if (pIg == PAM) { /* global ordering for selected pressure only */ 1075 if (schp) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Need the entire matrix"); 1076 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 1077 C = PPmat; 1078 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Unable to use the pressure matrix"); 1079 1080 ierr = KSPFETIDPSetPressureOperator(ksp,C);CHKERRQ(ierr); 1081 ierr = MatDestroy(&C);CHKERRQ(ierr); 1082 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing Pmat for pressure block"); 1083 } else { /* totP == 0 */ 1084 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr); 1085 } 1086 } 1087 PetscFunctionReturn(0); 1088 } 1089 1090 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 1091 { 1092 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1093 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1094 PetscBool flg; 1095 PetscErrorCode ierr; 1096 1097 PetscFunctionBegin; 1098 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 1099 /* set up BDDC */ 1100 ierr = PCSetErrorIfFailure(fetidp->innerbddc,ksp->errorifnotconverged);CHKERRQ(ierr); 1101 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 1102 /* FETI-DP as it is implemented needs an exact coarse solver */ 1103 if (pcbddc->coarse_ksp) { 1104 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1105 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1106 } 1107 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 1108 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 1109 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 1110 1111 /* setup FETI-DP operators 1112 If fetidp->statechanged is true, we need to update the operators 1113 needed in the saddle-point case. This should be replaced 1114 by a better logic when the FETI-DP matrix and preconditioner will 1115 have their own classes */ 1116 if (pcbddc->new_primal_space || fetidp->statechanged) { 1117 Mat F; /* the FETI-DP matrix */ 1118 PC D; /* the FETI-DP preconditioner */ 1119 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 1120 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 1121 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 1122 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 1123 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 1124 ierr = PetscObjectIncrementTabLevel((PetscObject)D,(PetscObject)fetidp->innerksp,0);CHKERRQ(ierr); 1125 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 1126 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 1127 ierr = MatDestroy(&F);CHKERRQ(ierr); 1128 ierr = PCDestroy(&D);CHKERRQ(ierr); 1129 if (fetidp->check) { 1130 PetscViewer viewer; 1131 1132 if (!pcbddc->dbg_viewer) { 1133 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)); 1134 } else { 1135 viewer = pcbddc->dbg_viewer; 1136 } 1137 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr); 1138 } 1139 } 1140 fetidp->statechanged = PETSC_FALSE; 1141 pcbddc->new_primal_space = PETSC_FALSE; 1142 1143 /* propagate settings to the inner solve */ 1144 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 1145 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 1146 if (ksp->res_hist) { 1147 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 1148 } 1149 ierr = KSPSetErrorIfNotConverged(fetidp->innerksp,ksp->errorifnotconverged);CHKERRQ(ierr); 1150 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 1151 PetscFunctionReturn(0); 1152 } 1153 1154 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 1155 { 1156 PetscErrorCode ierr; 1157 Mat F,A; 1158 MatNullSpace nsp; 1159 Vec X,B,Xl,Bl; 1160 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1161 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1162 1163 PetscFunctionBegin; 1164 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1165 if (fetidp->saddlepoint) { 1166 ierr = PetscCitationsRegister(citation2,&cited2);CHKERRQ(ierr); 1167 } 1168 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 1169 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 1170 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 1171 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 1172 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 1173 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 1174 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 1175 if (ksp->transpose_solve) { 1176 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1177 } else { 1178 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1179 } 1180 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 1181 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr); 1182 if (nsp) { 1183 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr); 1184 } 1185 /* update ksp with stats from inner ksp */ 1186 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 1187 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 1188 ksp->totalits += ksp->its; 1189 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 1190 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 1191 pcbddc->temp_solution_used = PETSC_FALSE; 1192 pcbddc->rhs_change = PETSC_FALSE; 1193 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1194 PetscFunctionReturn(0); 1195 } 1196 1197 static PetscErrorCode KSPReset_FETIDP(KSP ksp) 1198 { 1199 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1200 PC_BDDC *pcbddc; 1201 PetscErrorCode ierr; 1202 1203 PetscFunctionBegin; 1204 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr); 1205 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr); 1206 /* avoid PCReset that does not take into account ref counting */ 1207 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1208 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1209 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1210 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1211 pcbddc->symmetric_primal = PETSC_FALSE; 1212 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1213 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1214 fetidp->saddlepoint = PETSC_FALSE; 1215 fetidp->matstate = -1; 1216 fetidp->matnnzstate = -1; 1217 fetidp->statechanged = PETSC_TRUE; 1218 PetscFunctionReturn(0); 1219 } 1220 1221 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 1222 { 1223 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1224 PetscErrorCode ierr; 1225 1226 PetscFunctionBegin; 1227 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr); 1228 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1229 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1230 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 1231 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 1232 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 1233 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 1234 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",NULL);CHKERRQ(ierr); 1235 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 1236 PetscFunctionReturn(0); 1237 } 1238 1239 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 1240 { 1241 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1242 PetscErrorCode ierr; 1243 PetscBool iascii; 1244 1245 PetscFunctionBegin; 1246 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1247 if (iascii) { 1248 ierr = PetscViewerASCIIPrintf(viewer," fully redundant: %d\n",fetidp->fully_redundant);CHKERRQ(ierr); 1249 ierr = PetscViewerASCIIPrintf(viewer," saddle point: %d\n",fetidp->saddlepoint);CHKERRQ(ierr); 1250 ierr = PetscViewerASCIIPrintf(viewer,"Inner KSP solver details\n");CHKERRQ(ierr); 1251 } 1252 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1253 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 1254 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1255 if (iascii) { 1256 ierr = PetscViewerASCIIPrintf(viewer,"Inner BDDC solver details\n");CHKERRQ(ierr); 1257 } 1258 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1259 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 1260 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 1265 { 1266 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1267 PetscErrorCode ierr; 1268 1269 PetscFunctionBegin; 1270 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 1271 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1272 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 1273 if (!fetidp->userbddc) { 1274 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1275 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 1276 } 1277 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 1278 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 1279 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 1280 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr); 1281 ierr = PetscOptionsTail();CHKERRQ(ierr); 1282 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 1286 /*MC 1287 KSPFETIDP - The FETI-DP method 1288 1289 This class implements the FETI-DP method [1]. 1290 The matrix for the KSP must be of type MATIS. 1291 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an internal KSP object. 1292 1293 Options Database Keys: 1294 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 1295 . -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 1296 . -ksp_fetidp_saddlepoint_flip <false> : usually, an incompressible Stokes problem is written as 1297 | A B^T | | v | = | f | 1298 | B 0 | | p | = | g | 1299 with B representing -\int_\Omega \nabla \cdot u q. 1300 If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as 1301 | A B^T | | v | = | f | 1302 |-B 0 | | p | = |-g | 1303 . -ksp_fetidp_pressure_field <-1> : activates support for saddle point problems, and identifies the pressure field id. 1304 If this information is not provided, the pressure field is detected by using MatFindZeroDiagonals(). 1305 - -ksp_fetidp_pressure_all <false> : if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field. 1306 1307 Level: Advanced 1308 1309 Notes: 1310 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., 1311 .vb 1312 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 1313 .ve 1314 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 1315 1316 For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with KSPFETIDPSetPressureOperator(). 1317 Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix). 1318 If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver. 1319 Options for the pressure solver can be prefixed with -fetidp_fielsplit_p_, E.g. 1320 .vb 1321 -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps 1322 .ve 1323 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 1324 non-redundant multipliers, i.e. -ksp_fetidp_fullyredundant false. Options for the scaling solver are prefixed by -fetidp_bddelta_, E.g. 1325 .vb 1326 -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu 1327 .ve 1328 1329 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. 1330 1331 The converged reason and number of iterations computed are passed from the inner KSP to this KSP at the end of the solution. 1332 1333 Developer Notes: 1334 Even though this method does not directly use any norms, the user is allowed to set the KSPNormType to any value. 1335 This is so users do not have to change KSPNormType options when they switch from other KSP methods to this one. 1336 1337 References: 1338 .vb 1339 . [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 1340 . [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 1341 .ve 1342 1343 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC(), KSPFETIDPGetInnerBDDC(), KSPFETIDPGetInnerKSP() 1344 M*/ 1345 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp) 1346 { 1347 PetscErrorCode ierr; 1348 KSP_FETIDP *fetidp; 1349 KSP_FETIDPMon *monctx; 1350 PC_BDDC *pcbddc; 1351 PC pc; 1352 1353 PetscFunctionBegin; 1354 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,3);CHKERRQ(ierr); 1355 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,2);CHKERRQ(ierr); 1356 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1357 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1358 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 1359 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);CHKERRQ(ierr); 1360 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr); 1361 1362 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr); 1363 fetidp->matstate = -1; 1364 fetidp->matnnzstate = -1; 1365 fetidp->statechanged = PETSC_TRUE; 1366 1367 ksp->data = (void*)fetidp; 1368 ksp->ops->setup = KSPSetUp_FETIDP; 1369 ksp->ops->solve = KSPSolve_FETIDP; 1370 ksp->ops->destroy = KSPDestroy_FETIDP; 1371 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP; 1372 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP; 1373 ksp->ops->view = KSPView_FETIDP; 1374 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP; 1375 ksp->ops->buildsolution = KSPBuildSolution_FETIDP; 1376 ksp->ops->buildresidual = KSPBuildResidualDefault; 1377 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 1378 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1379 /* create the inner KSP for the Lagrange multipliers */ 1380 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr); 1381 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr); 1382 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1383 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr); 1384 /* monitor */ 1385 ierr = PetscNew(&monctx);CHKERRQ(ierr); 1386 monctx->parentksp = ksp; 1387 fetidp->monctx = monctx; 1388 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr); 1389 /* create the inner BDDC */ 1390 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1391 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1392 /* make sure we always obtain a consistent FETI-DP matrix 1393 for symmetric problems, the user can always customize it through the command line */ 1394 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1395 pcbddc->symmetric_primal = PETSC_FALSE; 1396 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1397 /* composed functions */ 1398 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr); 1399 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr); 1400 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr); 1401 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperator_C",KSPFETIDPSetPressureOperator_FETIDP);CHKERRQ(ierr); 1402 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */ 1403 ksp->setupnewmatrix = PETSC_TRUE; 1404 PetscFunctionReturn(0); 1405 } 1406