1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/ 2 #include <../src/ksp/pc/impls/bddc/bddc.h> 3 4 /* 5 This file implements the FETI-DP method in PETSc as part of KSP. 6 */ 7 typedef struct { 8 KSP parentksp; 9 } KSP_FETIDPMon; 10 11 typedef struct { 12 KSP innerksp; /* the KSP for the Lagrange multipliers */ 13 PC innerbddc; /* the inner BDDC object */ 14 PetscBool fully_redundant; /* true for using a fully redundant set of multipliers */ 15 PetscBool userbddc; /* true if the user provided the PCBDDC object */ 16 PetscBool saddlepoint; /* support for saddle point problems */ 17 KSP_FETIDPMon *monctx; /* monitor context, used to pass user defined monitors 18 in the physical space */ 19 } KSP_FETIDP; 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "KSPFETIDPSetPressureOperators_FETIDP" 23 static PetscErrorCode KSPFETIDPSetPressureOperators_FETIDP(KSP ksp, Mat A, Mat P) 24 { 25 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 if (A || P) fetidp->saddlepoint = PETSC_TRUE; 30 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_AAmat",(PetscObject)A);CHKERRQ(ierr); 31 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)P);CHKERRQ(ierr); 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "KSPFETIDPSetPressureOperators" 37 /*@ 38 KSPFETIDPSetPressureOperators - Sets the operators used to setup the pressure preconditioner for saddle point FETI-DP. 39 40 Collective on KSP 41 42 Input Parameters: 43 + ksp - the FETI-DP Krylov solver 44 . A - the linear operator defining the pressure problem 45 - P - the linear operator to be preconditioned 46 47 Level: advanced 48 49 Notes: The operators can be either passed in monolithic global ordering or in interface pressure ordering. 50 In the latter case, the interface pressure ordering of dofs needs to satisfy 51 pid_1 < pid_2 iff gid_1 < gid_2 52 where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding 53 id in the global ordering. 54 55 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators 56 @*/ 57 PetscErrorCode KSPFETIDPSetPressureOperators(KSP ksp, Mat A, Mat P) 58 { 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 63 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 64 if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,3); 65 ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperators_C",(KSP,Mat,Mat),(ksp,A,P));CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "KSPFETIDPGetInnerKSP_FETIDP" 71 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp) 72 { 73 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 74 75 PetscFunctionBegin; 76 *innerksp = fetidp->innerksp; 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "KSPFETIDPGetInnerKSP" 82 /*@ 83 KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers 84 85 Input Parameters: 86 + ksp - the FETI-DP KSP 87 - innerksp - the KSP for the multipliers 88 89 Level: advanced 90 91 Notes: 92 93 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC 94 @*/ 95 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp) 96 { 97 PetscErrorCode ierr; 98 99 PetscFunctionBegin; 100 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 101 PetscValidPointer(innerksp,2); 102 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr); 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNCT__ 107 #define __FUNCT__ "KSPFETIDPGetInnerBDDC_FETIDP" 108 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc) 109 { 110 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 111 112 PetscFunctionBegin; 113 *pc = fetidp->innerbddc; 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "KSPFETIDPGetInnerBDDC" 119 /*@ 120 KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 121 122 Input Parameters: 123 + ksp - the FETI-DP Krylov solver 124 - pc - the BDDC preconditioner 125 126 Level: advanced 127 128 Notes: 129 130 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP 131 @*/ 132 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc) 133 { 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 138 PetscValidPointer(pc,2); 139 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "KSPFETIDPSetInnerBDDC_FETIDP" 145 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc) 146 { 147 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 148 PetscErrorCode ierr; 149 150 PetscFunctionBegin; 151 ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr); 152 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 153 fetidp->innerbddc = pc; 154 fetidp->userbddc = PETSC_TRUE; 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "KSPFETIDPSetInnerBDDC" 160 /*@ 161 KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 162 163 Collective on KSP 164 165 Input Parameters: 166 + ksp - the FETI-DP Krylov solver 167 - pc - the BDDC preconditioner 168 169 Level: advanced 170 171 Notes: 172 173 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP 174 @*/ 175 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc) 176 { 177 PetscBool isbddc; 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 182 PetscValidHeaderSpecific(pc,PC_CLASSID,2); 183 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr); 184 if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner"); 185 ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "KSPBuildSolution_FETIDP" 191 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V) 192 { 193 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 194 Mat F; 195 Vec Xl; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 200 ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr); 201 if (v) { 202 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr); 203 *V = v; 204 } else { 205 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "KSPMonitor_FETIDP" 212 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx) 213 { 214 KSP_FETIDPMon *monctx = (KSP_FETIDPMon*)ctx; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "KSPComputeEigenvalues_FETIDP" 224 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig) 225 { 226 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "KSPComputeExtremeSingularValues_FETIDP" 236 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin) 237 { 238 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 239 PetscErrorCode ierr; 240 241 PetscFunctionBegin; 242 ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNCT__ 247 #define __FUNCT__ "KSPFETIDPSetUpOperators" 248 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp) 249 { 250 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 251 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 252 Mat A,Ap; 253 PetscInt fid = -1; 254 PetscBool ismatis,pisz; 255 PetscBool flip; /* Usually, Stokes is written 256 | A B'| | v | = | f | 257 | B 0 | | p | = | g | 258 If saddlepoint_flip is true, the code assumes it is written as 259 | A B'| | v | = | f | 260 |-B 0 | | p | = |-g | 261 */ 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 pisz = PETSC_FALSE; 266 flip = PETSC_FALSE; 267 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");CHKERRQ(ierr); 268 ierr = PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);CHKERRQ(ierr); 269 ierr = PetscOptionsBool("-ksp_fetidp_pressure_iszero","Zero pressure block",NULL,pisz,&pisz,NULL);CHKERRQ(ierr); 270 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);CHKERRQ(ierr); 271 ierr = PetscOptionsEnd();CHKERRQ(ierr); 272 273 fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint); 274 275 ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr); 276 ierr = PetscObjectTypeCompare((PetscObject)Ap,MATIS,&ismatis);CHKERRQ(ierr); 277 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pmat should be of type MATIS"); 278 279 /* see if MATIS has same fields attached */ 280 if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) { 281 PetscContainer c; 282 283 ierr = PetscObjectQuery((PetscObject)Ap,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr); 284 if (c) { 285 MatISLocalFields lf; 286 ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr); 287 ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr); 288 } 289 } 290 291 if (!fetidp->saddlepoint) { 292 ierr = PCSetOperators(fetidp->innerbddc,A,Ap);CHKERRQ(ierr); 293 } else { 294 KSP kP; 295 Mat nA,lA; 296 Mat PAmat,PPmat; 297 Vec rhs_flip; 298 IS pP,lP; 299 300 ierr = MatISGetLocalMat(Ap,&lA);CHKERRQ(ierr); 301 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr); 302 303 /* saved index sets */ 304 pP = NULL; 305 lP = NULL; 306 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP" ,(PetscObject*)&pP);CHKERRQ(ierr); 307 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP" ,(PetscObject*)&lP);CHKERRQ(ierr); 308 309 /* vector for flipping rhs */ 310 rhs_flip = NULL; 311 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip" ,(PetscObject*)&rhs_flip);CHKERRQ(ierr); 312 if (!pP && !lP) { /* first time, need to compute boundary pressure dofs */ 313 PC_IS *pcis = (PC_IS*)fetidp->innerbddc->data; 314 Mat_IS *matis = (Mat_IS*)(Ap->data); 315 ISLocalToGlobalMapping l2g; 316 IS II,pII,lPall,Pall,is1,is2; 317 const PetscInt *idxs; 318 PetscInt nl,ni,*widxs; 319 PetscInt i,j,n_neigh,*neigh,*n_shared,**shared,*count; 320 PetscInt rst,ren,n; 321 PetscBool ploc; 322 323 ierr = MatGetLocalSize(Ap,&nl,NULL);CHKERRQ(ierr); 324 ierr = MatGetOwnershipRange(Ap,&rst,&ren);CHKERRQ(ierr); 325 ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr); 326 327 if (!pcis->is_I_local) { /* need to compute interior dofs */ 328 ierr = PetscCalloc1(n,&count);CHKERRQ(ierr); 329 ierr = MatGetLocalToGlobalMapping(Ap,&l2g,NULL);CHKERRQ(ierr); 330 ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 331 for (i=1;i<n_neigh;i++) 332 for (j=0;j<n_shared[i];j++) 333 count[shared[i][j]] += 1; 334 for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i; 335 ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 336 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr); 337 } else { 338 ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr); 339 II = pcis->is_I_local; 340 } 341 342 /* interior dofs in layout */ 343 ierr = MatISSetUpSF(Ap);CHKERRQ(ierr); 344 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 345 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 346 ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr); 347 ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr); 348 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 349 ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr); 350 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 351 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 352 ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr); 353 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 354 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr); 355 356 /* pressure space at the interface */ 357 Pall = NULL; 358 lPall = NULL; 359 ploc = PETSC_FALSE; 360 if (fid >= 0) { 361 if (pcbddc->n_ISForDofsLocal) { 362 PetscInt np; 363 364 if (fid >= pcbddc->n_ISForDofsLocal) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofsLocal); 365 /* need a sequential IS */ 366 ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr); 367 ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 368 ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 369 ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 370 ploc = PETSC_TRUE; 371 } else if (pcbddc->n_ISForDofs) { 372 if (fid >= pcbddc->n_ISForDofs) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofs); 373 ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr); 374 Pall = pcbddc->ISForDofs[fid]; 375 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Missing fields! Use PCBDDCSetDofsSplitting/Local"); 376 } else { /* fallback to zero pressure block */ 377 IS list[2]; 378 379 ierr = MatFindZeroDiagonals(Ap,&list[1]);CHKERRQ(ierr); 380 ierr = ISComplement(list[1],rst,ren,&list[0]);CHKERRQ(ierr); 381 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,2,list);CHKERRQ(ierr); 382 ierr = ISDestroy(&list[0]);CHKERRQ(ierr); 383 Pall = list[1]; 384 } 385 if (ploc) { 386 ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr); 387 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 388 } else { 389 ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr); 390 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 391 /* need all local pressure dofs */ 392 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 393 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 394 ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr); 395 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 396 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 397 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 398 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 399 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 400 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 401 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 402 } 403 404 if (flip) { 405 PetscInt npl; 406 if (!Pall) { 407 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 408 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 409 ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr); 410 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr); 411 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 412 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr); 413 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 414 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 415 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 416 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr); 417 } 418 ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr); 419 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 420 ierr = MatCreateVecs(Ap,NULL,&rhs_flip);CHKERRQ(ierr); 421 ierr = VecSet(rhs_flip,1.);CHKERRQ(ierr); 422 ierr = VecSetOption(rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 423 for (i=0;i<npl;i++) { 424 ierr = VecSetValue(rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr); 425 } 426 ierr = VecAssemblyBegin(rhs_flip);CHKERRQ(ierr); 427 ierr = VecAssemblyEnd(rhs_flip);CHKERRQ(ierr); 428 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)rhs_flip);CHKERRQ(ierr); 429 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 430 } 431 ierr = ISDestroy(&Pall);CHKERRQ(ierr); 432 ierr = ISDestroy(&pII);CHKERRQ(ierr); 433 434 /* local interface pressures in subdomain-wise and global ordering */ 435 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 436 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 437 if (pP) { 438 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr); 439 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr); 440 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 441 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr); 442 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 443 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 444 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 445 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs+ni);CHKERRQ(ierr); 446 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr); 447 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 448 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs+ni,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 449 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 450 ierr = ISDestroy(&is1);CHKERRQ(ierr); 451 } else { 452 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr); 453 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 454 for (i=0;i<ni;i++) 455 if (idxs[i] >=0 && idxs[i] < n) 456 matis->sf_leafdata[idxs[i]] = 1; 457 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 458 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 459 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr); 460 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 461 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 462 ierr = ISDestroy(&is1);CHKERRQ(ierr); 463 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 464 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 465 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr); 466 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 467 } 468 ierr = PetscFree(widxs);CHKERRQ(ierr); 469 470 /* We may want to use a discrete harmonic solver instead 471 of a Stokes harmonic for the Dirichler preconditioner 472 Need to extract the interior pressure dofs in interior dofs ordering */ 473 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr); 474 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr); 475 ierr = ISDestroy(&is1);CHKERRQ(ierr); 476 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g);CHKERRQ(ierr); 477 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr); 478 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr); 479 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr); 480 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iP",i,j); 481 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr); 482 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 483 ierr = ISDestroy(&is1);CHKERRQ(ierr); 484 ierr = ISDestroy(&is2);CHKERRQ(ierr); 485 ierr = ISDestroy(&lPall);CHKERRQ(ierr); 486 ierr = ISDestroy(&II);CHKERRQ(ierr); 487 488 /* exclude interface pressures from the inner BDDC */ 489 if (pcbddc->DirichletBoundariesLocal) { 490 IS list[2],plP,isout; 491 PetscInt np; 492 493 /* need a parallel IS */ 494 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 495 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 496 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr); 497 list[0] = plP; 498 list[1] = pcbddc->DirichletBoundariesLocal; 499 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 500 ierr = ISDestroy(&plP);CHKERRQ(ierr); 501 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 502 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr); 503 ierr = ISDestroy(&isout);CHKERRQ(ierr); 504 } else if (pcbddc->DirichletBoundaries) { 505 IS list[2],isout; 506 507 list[0] = pP; 508 list[1] = pcbddc->DirichletBoundaries; 509 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 510 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr); 511 ierr = ISDestroy(&isout);CHKERRQ(ierr); 512 } else { 513 IS plP; 514 PetscInt np; 515 516 /* need a parallel IS */ 517 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 518 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 519 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr); 520 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr); 521 ierr = ISDestroy(&plP);CHKERRQ(ierr); 522 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 523 } 524 } else { 525 if (!pP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing global interface pressure field"); 526 if (!lP) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Missing local interface pressure field"); 527 ierr = PetscObjectReference((PetscObject)pP);CHKERRQ(ierr); 528 ierr = PetscObjectReference((PetscObject)lP);CHKERRQ(ierr); 529 if (rhs_flip) { 530 ierr = PetscObjectReference((PetscObject)rhs_flip);CHKERRQ(ierr); 531 } 532 } 533 534 /* Set operator for inner BDDC */ 535 ierr = MatDuplicate(Ap,MAT_COPY_VALUES,&nA);CHKERRQ(ierr); 536 if (rhs_flip) { 537 Mat lA2; 538 539 ierr = MatDiagonalScale(nA,rhs_flip,NULL);CHKERRQ(ierr); 540 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr); 541 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr); 542 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr); 543 ierr = MatDestroy(&lA2);CHKERRQ(ierr); 544 } 545 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 546 ierr = MatZeroRowsColumnsIS(nA,pP,1.,NULL,NULL);CHKERRQ(ierr); 547 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr); 548 ierr = MatDestroy(&nA);CHKERRQ(ierr); 549 ierr = VecDestroy(&rhs_flip);CHKERRQ(ierr); 550 551 /* non-zero rhs on interior dofs when applying the preconditioner */ 552 pcbddc->switch_static = PETSC_TRUE; 553 554 /* Operators for pressure preconditioner */ 555 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject*)&PAmat);CHKERRQ(ierr); 556 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 557 if (PAmat) { 558 ierr = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr); 559 } 560 if (PPmat) { 561 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 562 } 563 564 /* Extract pressure block */ 565 if (!pisz) { 566 Mat C; 567 568 ierr = MatGetSubMatrix(Ap,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 569 ierr = MatScale(C,-1.);CHKERRQ(ierr); 570 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr); 571 /* default operators for the interface pressure solver */ 572 if (!PAmat) { 573 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 574 PAmat = C; 575 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 576 } 577 if (!PPmat) { 578 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 579 PPmat = C; 580 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 581 } 582 ierr = MatDestroy(&C);CHKERRQ(ierr); 583 } 584 585 if (!PAmat) { /* Just use the identity */ 586 PetscInt nl; 587 588 ierr = ISGetLocalSize(pP,&nl);CHKERRQ(ierr); 589 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&PAmat);CHKERRQ(ierr); 590 ierr = MatSetSizes(PAmat,nl,nl,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 591 ierr = MatSetType(PAmat,MATAIJ);CHKERRQ(ierr); 592 ierr = MatMPIAIJSetPreallocation(PAmat,1,NULL,0,NULL);CHKERRQ(ierr); 593 ierr = MatSeqAIJSetPreallocation(PAmat,1,NULL);CHKERRQ(ierr); 594 ierr = MatAssemblyBegin(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 595 ierr = MatAssemblyEnd(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 596 ierr = MatShift(PAmat,1.);CHKERRQ(ierr); 597 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 598 } else { /* check size of pressure operator and restrict if needed */ 599 PetscInt AM,PAM,PAN,pam,pan,am,an,pl; 600 601 ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr); 602 ierr = MatGetSize(PAmat,&PAM,&PAN);CHKERRQ(ierr); 603 ierr = MatGetLocalSize(PAmat,&pam,&pan);CHKERRQ(ierr); 604 ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr); 605 ierr = ISGetLocalSize(pP,&pl);CHKERRQ(ierr); 606 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 607 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 608 if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl); 609 if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl); 610 if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */ 611 Mat C; 612 613 ierr = MatGetSubMatrix(PAmat,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 614 ierr = MatDestroy(&PAmat);CHKERRQ(ierr); 615 PAmat = C; 616 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 617 } 618 } 619 if (PPmat) { 620 PetscInt AM,PAM,PAN,pam,pan,am,an,pl; 621 622 ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr); 623 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 624 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 625 ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr); 626 ierr = ISGetLocalSize(pP,&pl);CHKERRQ(ierr); 627 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 628 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 629 if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl); 630 if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl); 631 if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */ 632 Mat C; 633 634 ierr = MatGetSubMatrix(PPmat,pP,pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 635 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 636 PPmat = C; 637 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 638 } 639 } else { 640 ierr = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr); 641 PPmat = PAmat; 642 } 643 644 /* create pressure solver */ 645 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&kP);CHKERRQ(ierr); 646 ierr = KSPSetOperators(kP,PAmat,PPmat);CHKERRQ(ierr); 647 ierr = KSPSetOptionsPrefix(kP,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 648 ierr = KSPAppendOptionsPrefix(kP,"fetidp_p_");CHKERRQ(ierr); 649 ierr = KSPSetFromOptions(kP);CHKERRQ(ierr); 650 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PKSP",(PetscObject)kP);CHKERRQ(ierr); 651 ierr = MatDestroy(&PAmat);CHKERRQ(ierr); 652 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 653 ierr = KSPDestroy(&kP);CHKERRQ(ierr); 654 655 ierr = ISDestroy(&lP);CHKERRQ(ierr); 656 ierr = ISDestroy(&pP);CHKERRQ(ierr); 657 } 658 PetscFunctionReturn(0); 659 } 660 661 #undef __FUNCT__ 662 #define __FUNCT__ "KSPSetUp_FETIDP" 663 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 664 { 665 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 666 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 667 PetscBool flg; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 672 /* set up BDDC */ 673 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 674 675 /* FETI-DP as it is implemented needs an exact coarse solver */ 676 if (pcbddc->coarse_ksp) { 677 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 678 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 679 } 680 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 681 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 682 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 683 684 /* if the primal space is changed, setup F */ 685 if (pcbddc->new_primal_space || !pcbddc->coarse_size) { 686 Mat F; /* the FETI-DP matrix */ 687 PC D; /* the FETI-DP preconditioner */ 688 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 689 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 690 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 691 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 692 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 693 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 694 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 695 ierr = MatDestroy(&F);CHKERRQ(ierr); 696 ierr = PCDestroy(&D);CHKERRQ(ierr); 697 } 698 699 /* propagate settings to the inner solve */ 700 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 701 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 702 if (ksp->res_hist) { 703 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 704 } 705 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNCT__ 710 #define __FUNCT__ "KSPSolve_FETIDP" 711 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 712 { 713 PetscErrorCode ierr; 714 Mat F; 715 Vec X,B,Xl,Bl; 716 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 717 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 718 719 PetscFunctionBegin; 720 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 721 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 722 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 723 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 724 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 725 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 726 if (ksp->transpose_solve) { 727 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 728 } else { 729 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 730 } 731 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 732 /* update ksp with stats from inner ksp */ 733 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 734 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 735 ksp->totalits += ksp->its; 736 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 737 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 738 pcbddc->temp_solution_used = PETSC_FALSE; 739 pcbddc->rhs_change = PETSC_FALSE; 740 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 741 PetscFunctionReturn(0); 742 } 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "KSPDestroy_FETIDP" 746 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 747 { 748 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 749 PetscErrorCode ierr; 750 751 PetscFunctionBegin; 752 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 753 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 754 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 755 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 756 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 757 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 758 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",NULL);CHKERRQ(ierr); 759 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 760 PetscFunctionReturn(0); 761 } 762 763 #undef __FUNCT__ 764 #define __FUNCT__ "KSPView_FETIDP" 765 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 766 { 767 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 768 PetscErrorCode ierr; 769 PetscBool iascii; 770 771 PetscFunctionBegin; 772 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 773 if (iascii) { 774 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: fully redundant: %D\n",fetidp->fully_redundant);CHKERRQ(ierr); 775 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: saddle point: %D\n",fetidp->saddlepoint);CHKERRQ(ierr); 776 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: inner solver details\n");CHKERRQ(ierr); 777 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 778 } 779 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 780 if (iascii) { 781 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 782 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: BDDC solver details\n");CHKERRQ(ierr); 783 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 784 } 785 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 786 if (iascii) { 787 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 788 } 789 PetscFunctionReturn(0); 790 } 791 792 #undef __FUNCT__ 793 #define __FUNCT__ "KSPSetFromOptions_FETIDP" 794 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 795 { 796 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 797 PetscErrorCode ierr; 798 799 PetscFunctionBegin; 800 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 801 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 802 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 803 if (!fetidp->userbddc) { 804 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 805 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 806 } 807 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 808 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 809 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 810 ierr = PetscOptionsTail();CHKERRQ(ierr); 811 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 812 PetscFunctionReturn(0); 813 } 814 815 /*MC 816 KSPFETIDP - The FETI-DP method 817 818 This class implements the FETI-DP method [1]. 819 The preconditioning matrix for the KSP must be of type MATIS. 820 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an inner KSP object. 821 822 Options Database Keys: 823 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 824 - -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 825 826 Level: Advanced 827 828 Notes: Options for the inner KSP and for the customization of the PCBDDC object can be specified at command line by using the prefixes -fetidp_ and -fetidp_bddc_. E.g., 829 .vb 830 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 831 .ve 832 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 833 For saddle point problems with continuous pressures, the operators for the interface pressure solver can be specified with KSPFETIDPSetPressureOperators(). 834 835 References: 836 .vb 837 . [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 838 . [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 839 .ve 840 841 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP 842 M*/ 843 #undef __FUNCT__ 844 #define __FUNCT__ "KSPCreate_FETIDP" 845 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp) 846 { 847 PetscErrorCode ierr; 848 KSP_FETIDP *fetidp; 849 KSP_FETIDPMon *monctx; 850 PC_BDDC *pcbddc; 851 PC pc; 852 853 PetscFunctionBegin; 854 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr); 855 ksp->data = (void*)fetidp; 856 ksp->ops->setup = KSPSetUp_FETIDP; 857 ksp->ops->solve = KSPSolve_FETIDP; 858 ksp->ops->destroy = KSPDestroy_FETIDP; 859 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP; 860 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP; 861 ksp->ops->view = KSPView_FETIDP; 862 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP; 863 ksp->ops->buildsolution = KSPBuildSolution_FETIDP; 864 ksp->ops->buildresidual = KSPBuildResidualDefault; 865 /* create the inner KSP for the Lagrange multipliers */ 866 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr); 867 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr); 868 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 869 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr); 870 /* monitor */ 871 ierr = PetscNew(&monctx);CHKERRQ(ierr); 872 monctx->parentksp = ksp; 873 fetidp->monctx = monctx; 874 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr); 875 /* create the inner BDDC */ 876 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 877 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 878 /* make sure we always obtain a consistent FETI-DP matrix 879 for symmetric problems, the user can always customize it through the command line */ 880 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 881 pcbddc->symmetric_primal = PETSC_FALSE; 882 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 883 /* composed functions */ 884 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr); 885 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr); 886 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr); 887 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",KSPFETIDPSetPressureOperators_FETIDP);CHKERRQ(ierr); 888 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */ 889 ksp->setupnewmatrix = PETSC_TRUE; 890 PetscFunctionReturn(0); 891 } 892