1 2 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 4 #include <petscdm.h> 5 6 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","SELFP","A11","USER","FULL","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 7 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactType","PC_FIELDSPLIT_SCHUR_FACT_",0}; 8 9 PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_S,KSP_Solve_FS_U,KSP_Solve_FS_L,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4; 10 11 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 12 struct _PC_FieldSplitLink { 13 KSP ksp; 14 Vec x,y,z; 15 char *splitname; 16 PetscInt nfields; 17 PetscInt *fields,*fields_col; 18 VecScatter sctx; 19 IS is,is_col,is_orig; 20 PC_FieldSplitLink next,previous; 21 PetscLogEvent event; 22 }; 23 24 typedef struct { 25 PCCompositeType type; 26 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 27 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 28 PetscInt bs; /* Block size for IS and Mat structures */ 29 PetscInt nsplits; /* Number of field divisions defined */ 30 Vec *x,*y,w1,w2; 31 Mat *mat; /* The diagonal block for each split */ 32 Mat *pmat; /* The preconditioning diagonal block for each split */ 33 Mat *Afield; /* The rows of the matrix associated with each split */ 34 PetscBool issetup; 35 36 /* Only used when Schur complement preconditioning is used */ 37 Mat B; /* The (0,1) block */ 38 Mat C; /* The (1,0) block */ 39 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 40 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 41 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 42 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 43 PCFieldSplitSchurFactType schurfactorization; 44 KSP kspschur; /* The solver for S */ 45 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 46 PC_FieldSplitLink head; 47 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 48 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 49 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 50 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 51 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 52 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 53 } PC_FieldSplit; 54 55 /* 56 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 57 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 58 PC you could change this. 59 */ 60 61 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 62 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 63 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 64 { 65 switch (jac->schurpre) { 66 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 67 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: return jac->schurp; 68 case PC_FIELDSPLIT_SCHUR_PRE_A11: return jac->pmat[1]; 69 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 70 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 71 default: 72 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 73 } 74 } 75 76 77 #include <petscdraw.h> 78 #undef __FUNCT__ 79 #define __FUNCT__ "PCView_FieldSplit" 80 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 81 { 82 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 83 PetscErrorCode ierr; 84 PetscBool iascii,isdraw; 85 PetscInt i,j; 86 PC_FieldSplitLink ilink = jac->head; 87 88 PetscFunctionBegin; 89 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 90 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 91 if (iascii) { 92 if (jac->bs > 0) { 93 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 94 } else { 95 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D\n",PCCompositeTypes[jac->type],jac->nsplits);CHKERRQ(ierr); 96 } 97 if (pc->useAmat) { 98 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 99 } 100 if (jac->diag_use_amat) { 101 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for diagonal blocks\n");CHKERRQ(ierr); 102 } 103 if (jac->offdiag_use_amat) { 104 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for off-diagonal blocks\n");CHKERRQ(ierr); 105 } 106 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 107 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 108 for (i=0; i<jac->nsplits; i++) { 109 if (ilink->fields) { 110 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 111 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 112 for (j=0; j<ilink->nfields; j++) { 113 if (j > 0) { 114 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 115 } 116 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 117 } 118 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 119 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 120 } else { 121 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 122 } 123 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 124 ilink = ilink->next; 125 } 126 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 127 } 128 129 if (isdraw) { 130 PetscDraw draw; 131 PetscReal x,y,w,wd; 132 133 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 134 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 135 w = 2*PetscMin(1.0 - x,x); 136 wd = w/(jac->nsplits + 1); 137 x = x - wd*(jac->nsplits-1)/2.0; 138 for (i=0; i<jac->nsplits; i++) { 139 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 140 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 141 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 142 x += wd; 143 ilink = ilink->next; 144 } 145 } 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "PCView_FieldSplit_Schur" 151 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 152 { 153 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 154 PetscErrorCode ierr; 155 PetscBool iascii,isdraw; 156 PetscInt i,j; 157 PC_FieldSplitLink ilink = jac->head; 158 159 PetscFunctionBegin; 160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 161 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 162 if (iascii) { 163 if (jac->bs > 0) { 164 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 165 } else { 166 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, factorization %s\n",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 167 } 168 if (pc->useAmat) { 169 ierr = PetscViewerASCIIPrintf(viewer," using Amat (not Pmat) as operator for blocks\n");CHKERRQ(ierr); 170 } 171 switch (jac->schurpre) { 172 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 173 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from S itself\n");CHKERRQ(ierr);break; 174 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 175 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses (lumped, if requested) A00's diagonal's inverse\n");CHKERRQ(ierr);break; 176 case PC_FIELDSPLIT_SCHUR_PRE_A11: 177 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr);break; 178 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 179 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from the exact Schur complement\n");CHKERRQ(ierr);break; 180 case PC_FIELDSPLIT_SCHUR_PRE_USER: 181 if (jac->schur_user) { 182 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from user provided matrix\n");CHKERRQ(ierr); 183 } else { 184 ierr = PetscViewerASCIIPrintf(viewer," Preconditioner for the Schur complement formed from A11\n");CHKERRQ(ierr); 185 } 186 break; 187 default: 188 SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 189 } 190 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 191 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 192 for (i=0; i<jac->nsplits; i++) { 193 if (ilink->fields) { 194 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 195 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 196 for (j=0; j<ilink->nfields; j++) { 197 if (j > 0) { 198 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 199 } 200 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 201 } 202 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 203 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 204 } else { 205 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 206 } 207 ilink = ilink->next; 208 } 209 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block\n");CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 211 if (jac->head) { 212 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 213 } else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 214 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 215 if (jac->head && jac->kspupper != jac->head->ksp) { 216 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for upper A00 in upper triangular factor \n");CHKERRQ(ierr); 217 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 218 if (jac->kspupper) {ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr);} 219 else {ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr);} 220 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 221 } 222 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 223 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 224 if (jac->kspschur) { 225 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 226 } else { 227 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 228 } 229 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 231 } else if (isdraw && jac->head) { 232 PetscDraw draw; 233 PetscReal x,y,w,wd,h; 234 PetscInt cnt = 2; 235 char str[32]; 236 237 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 238 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr); 239 if (jac->kspupper != jac->head->ksp) cnt++; 240 w = 2*PetscMin(1.0 - x,x); 241 wd = w/(cnt + 1); 242 243 ierr = PetscSNPrintf(str,32,"Schur fact. %s",PCFieldSplitSchurFactTypes[jac->schurfactorization]);CHKERRQ(ierr); 244 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 245 y -= h; 246 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 247 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]);CHKERRQ(ierr); 248 } else { 249 ierr = PetscSNPrintf(str,32,"Prec. for Schur from %s",PCFieldSplitSchurPreTypes[jac->schurpre]);CHKERRQ(ierr); 250 } 251 ierr = PetscDrawStringBoxed(draw,x+wd*(cnt-1)/2.0,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr); 252 y -= h; 253 x = x - wd*(cnt-1)/2.0; 254 255 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 256 ierr = KSPView(jac->head->ksp,viewer);CHKERRQ(ierr); 257 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 258 if (jac->kspupper != jac->head->ksp) { 259 x += wd; 260 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 261 ierr = KSPView(jac->kspupper,viewer);CHKERRQ(ierr); 262 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 263 } 264 x += wd; 265 ierr = PetscDrawPushCurrentPoint(draw,x,y);CHKERRQ(ierr); 266 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 267 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr); 268 } 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 274 /* Precondition: jac->bs is set to a meaningful value */ 275 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 276 { 277 PetscErrorCode ierr; 278 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 279 PetscInt i,nfields,*ifields,nfields_col,*ifields_col; 280 PetscBool flg,flg_col; 281 char optionname[128],splitname[8],optionname_col[128]; 282 283 PetscFunctionBegin; 284 ierr = PetscMalloc1(jac->bs,&ifields);CHKERRQ(ierr); 285 ierr = PetscMalloc1(jac->bs,&ifields_col);CHKERRQ(ierr); 286 for (i=0,flg=PETSC_TRUE;; i++) { 287 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 288 ierr = PetscSNPrintf(optionname,sizeof(optionname),"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 289 ierr = PetscSNPrintf(optionname_col,sizeof(optionname_col),"-pc_fieldsplit_%D_fields_col",i);CHKERRQ(ierr); 290 nfields = jac->bs; 291 nfields_col = jac->bs; 292 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 293 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,optionname_col,ifields_col,&nfields_col,&flg_col);CHKERRQ(ierr); 294 if (!flg) break; 295 else if (flg && !flg_col) { 296 if (!nfields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 297 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields);CHKERRQ(ierr); 298 } else { 299 if (!nfields || !nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot list zero fields"); 300 if (nfields != nfields_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Number of row and column fields must match"); 301 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields,ifields_col);CHKERRQ(ierr); 302 } 303 } 304 if (i > 0) { 305 /* Makes command-line setting of splits take precedence over setting them in code. 306 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 307 create new splits, which would probably not be what the user wanted. */ 308 jac->splitdefined = PETSC_TRUE; 309 } 310 ierr = PetscFree(ifields);CHKERRQ(ierr); 311 ierr = PetscFree(ifields_col);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "PCFieldSplitSetDefaults" 317 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 318 { 319 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 320 PetscErrorCode ierr; 321 PC_FieldSplitLink ilink = jac->head; 322 PetscBool fieldsplit_default = PETSC_FALSE,stokes = PETSC_FALSE,coupling = PETSC_FALSE; 323 PetscInt i; 324 325 PetscFunctionBegin; 326 /* 327 Kinda messy, but at least this now uses DMCreateFieldDecomposition() even with jac->reset. 328 Should probably be rewritten. 329 */ 330 if (!ilink || jac->reset) { 331 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 332 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_coupling",&coupling,NULL);CHKERRQ(ierr); 333 if (pc->dm && jac->dm_splits && !stokes && !coupling) { 334 PetscInt numFields, f, i, j; 335 char **fieldNames; 336 IS *fields; 337 DM *dms; 338 DM subdm[128]; 339 PetscBool flg; 340 341 ierr = DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms);CHKERRQ(ierr); 342 /* Allow the user to prescribe the splits */ 343 for (i = 0, flg = PETSC_TRUE;; i++) { 344 PetscInt ifields[128]; 345 IS compField; 346 char optionname[128], splitname[8]; 347 PetscInt nfields = numFields; 348 349 ierr = PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%D_fields", i);CHKERRQ(ierr); 350 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg);CHKERRQ(ierr); 351 if (!flg) break; 352 if (numFields > 128) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot currently support %d > 128 fields", numFields); 353 ierr = DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]);CHKERRQ(ierr); 354 if (nfields == 1) { 355 ierr = PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField);CHKERRQ(ierr); 356 } else { 357 ierr = PetscSNPrintf(splitname, sizeof(splitname), "%D", i);CHKERRQ(ierr); 358 ierr = PCFieldSplitSetIS(pc, splitname, compField);CHKERRQ(ierr); 359 } 360 ierr = ISDestroy(&compField);CHKERRQ(ierr); 361 for (j = 0; j < nfields; ++j) { 362 f = ifields[j]; 363 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 364 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 365 } 366 } 367 if (i == 0) { 368 for (f = 0; f < numFields; ++f) { 369 ierr = PCFieldSplitSetIS(pc, fieldNames[f], fields[f]);CHKERRQ(ierr); 370 ierr = PetscFree(fieldNames[f]);CHKERRQ(ierr); 371 ierr = ISDestroy(&fields[f]);CHKERRQ(ierr); 372 } 373 } else { 374 for (j=0; j<numFields; j++) { 375 ierr = DMDestroy(dms+j);CHKERRQ(ierr); 376 } 377 ierr = PetscFree(dms);CHKERRQ(ierr); 378 ierr = PetscMalloc1(i, &dms);CHKERRQ(ierr); 379 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 380 } 381 ierr = PetscFree(fieldNames);CHKERRQ(ierr); 382 ierr = PetscFree(fields);CHKERRQ(ierr); 383 if (dms) { 384 ierr = PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 385 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 386 const char *prefix; 387 ierr = PetscObjectGetOptionsPrefix((PetscObject)(ilink->ksp),&prefix);CHKERRQ(ierr); 388 ierr = PetscObjectSetOptionsPrefix((PetscObject)(dms[i]), prefix);CHKERRQ(ierr); 389 ierr = KSPSetDM(ilink->ksp, dms[i]);CHKERRQ(ierr); 390 ierr = KSPSetDMActive(ilink->ksp, PETSC_FALSE);CHKERRQ(ierr); 391 ierr = PetscObjectIncrementTabLevel((PetscObject)dms[i],(PetscObject)ilink->ksp,0);CHKERRQ(ierr); 392 ierr = DMDestroy(&dms[i]);CHKERRQ(ierr); 393 } 394 ierr = PetscFree(dms);CHKERRQ(ierr); 395 } 396 } else { 397 if (jac->bs <= 0) { 398 if (pc->pmat) { 399 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 400 } else jac->bs = 1; 401 } 402 403 if (stokes) { 404 IS zerodiags,rest; 405 PetscInt nmin,nmax; 406 407 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 408 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 409 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 410 if (jac->reset) { 411 if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup"); 412 jac->head->is = rest; 413 jac->head->next->is = zerodiags; 414 } else { 415 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 416 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 417 } 418 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 419 ierr = ISDestroy(&rest);CHKERRQ(ierr); 420 } else if (coupling) { 421 IS coupling,rest; 422 PetscInt nmin,nmax; 423 424 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 425 ierr = MatFindOffBlockDiagonalEntries(pc->mat,&coupling);CHKERRQ(ierr); 426 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc->mat),nmax-nmin,nmin,1,&rest);CHKERRQ(ierr); 427 ierr = ISSetIdentity(rest);CHKERRQ(ierr); 428 if (jac->reset) { 429 if (!jac->head) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"At reset jac must have head from previous setup"); 430 jac->head->is = rest; 431 jac->head->next->is = coupling; 432 } else { 433 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 434 ierr = PCFieldSplitSetIS(pc,"1",coupling);CHKERRQ(ierr); 435 } 436 ierr = ISDestroy(&coupling);CHKERRQ(ierr); 437 ierr = ISDestroy(&rest);CHKERRQ(ierr); 438 } else { 439 if (jac->reset && !jac->isrestrict) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 440 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&fieldsplit_default,NULL);CHKERRQ(ierr); 441 if (!fieldsplit_default) { 442 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 443 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 444 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 445 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 446 } 447 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 448 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 449 for (i=0; i<jac->bs; i++) { 450 char splitname[8]; 451 ierr = PetscSNPrintf(splitname,sizeof(splitname),"%D",i);CHKERRQ(ierr); 452 ierr = PCFieldSplitSetFields(pc,splitname,1,&i,&i);CHKERRQ(ierr); 453 } 454 jac->defaultsplit = PETSC_TRUE; 455 } 456 } 457 } 458 } else if (jac->nsplits == 1) { 459 if (ilink->is) { 460 IS is2; 461 PetscInt nmin,nmax; 462 463 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 464 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 465 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 466 ierr = ISDestroy(&is2);CHKERRQ(ierr); 467 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 468 } 469 470 if (jac->nsplits < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 471 PetscFunctionReturn(0); 472 } 473 474 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions,const char pre[], const char name[], char *value[], PetscBool *flg); 475 476 #undef __FUNCT__ 477 #define __FUNCT__ "PCSetUp_FieldSplit" 478 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 479 { 480 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 481 PetscErrorCode ierr; 482 PC_FieldSplitLink ilink; 483 PetscInt i,nsplit; 484 PetscBool sorted, sorted_col; 485 486 PetscFunctionBegin; 487 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 488 nsplit = jac->nsplits; 489 ilink = jac->head; 490 491 /* get the matrices for each split */ 492 if (!jac->issetup) { 493 PetscInt rstart,rend,nslots,bs; 494 495 jac->issetup = PETSC_TRUE; 496 497 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 498 if (jac->defaultsplit || !ilink->is) { 499 if (jac->bs <= 0) jac->bs = nsplit; 500 } 501 bs = jac->bs; 502 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 503 nslots = (rend - rstart)/bs; 504 for (i=0; i<nsplit; i++) { 505 if (jac->defaultsplit) { 506 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 507 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 508 } else if (!ilink->is) { 509 if (ilink->nfields > 1) { 510 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 511 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 512 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 513 for (j=0; j<nslots; j++) { 514 for (k=0; k<nfields; k++) { 515 ii[nfields*j + k] = rstart + bs*j + fields[k]; 516 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 517 } 518 } 519 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 520 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 521 ierr = ISSetBlockSize(ilink->is, nfields);CHKERRQ(ierr); 522 ierr = ISSetBlockSize(ilink->is_col, nfields);CHKERRQ(ierr); 523 } else { 524 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 525 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 526 } 527 } 528 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 529 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 530 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 531 ilink = ilink->next; 532 } 533 } 534 535 ilink = jac->head; 536 if (!jac->pmat) { 537 Vec xtmp; 538 539 ierr = MatCreateVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 540 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 541 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 542 for (i=0; i<nsplit; i++) { 543 MatNullSpace sp; 544 545 /* Check for preconditioning matrix attached to IS */ 546 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 547 if (jac->pmat[i]) { 548 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 549 if (jac->type == PC_COMPOSITE_SCHUR) { 550 jac->schur_user = jac->pmat[i]; 551 552 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 553 } 554 } else { 555 const char *prefix; 556 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 557 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 558 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 559 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 560 } 561 /* create work vectors for each split */ 562 ierr = MatCreateVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 563 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 564 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 565 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 566 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 567 if (sp) { 568 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 569 } 570 ilink = ilink->next; 571 } 572 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 573 } else { 574 for (i=0; i<nsplit; i++) { 575 Mat pmat; 576 577 /* Check for preconditioning matrix attached to IS */ 578 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 579 if (!pmat) { 580 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 581 } 582 ilink = ilink->next; 583 } 584 } 585 if (jac->diag_use_amat) { 586 ilink = jac->head; 587 if (!jac->mat) { 588 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 589 for (i=0; i<nsplit; i++) { 590 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 591 ilink = ilink->next; 592 } 593 } else { 594 for (i=0; i<nsplit; i++) { 595 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 596 ilink = ilink->next; 597 } 598 } 599 } else { 600 jac->mat = jac->pmat; 601 } 602 603 /* Check for null space attached to IS */ 604 ilink = jac->head; 605 for (i=0; i<nsplit; i++) { 606 MatNullSpace sp; 607 608 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 609 if (sp) { 610 ierr = MatSetNullSpace(jac->mat[i], sp);CHKERRQ(ierr); 611 } 612 ilink = ilink->next; 613 } 614 615 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 616 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 617 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 618 ilink = jac->head; 619 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 620 /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */ 621 if (!jac->Afield) { 622 ierr = PetscCalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 623 if (jac->offdiag_use_amat) { 624 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 625 } else { 626 ierr = MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 627 } 628 } else { 629 if (jac->offdiag_use_amat) { 630 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 631 } else { 632 ierr = MatGetSubMatrix(pc->pmat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 633 } 634 } 635 } else { 636 if (!jac->Afield) { 637 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 638 for (i=0; i<nsplit; i++) { 639 if (jac->offdiag_use_amat) { 640 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 641 } else { 642 ierr = MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 643 } 644 ilink = ilink->next; 645 } 646 } else { 647 for (i=0; i<nsplit; i++) { 648 if (jac->offdiag_use_amat) { 649 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 650 } else { 651 ierr = MatGetSubMatrix(pc->pmat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 652 } 653 ilink = ilink->next; 654 } 655 } 656 } 657 } 658 659 if (jac->type == PC_COMPOSITE_SCHUR) { 660 IS ccis; 661 PetscInt rstart,rend; 662 char lscname[256]; 663 PetscObject LSC_L; 664 665 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 666 667 /* When extracting off-diagonal submatrices, we take complements from this range */ 668 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 669 670 /* need to handle case when one is resetting up the preconditioner */ 671 if (jac->schur) { 672 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 673 674 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 675 ilink = jac->head; 676 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 677 if (jac->offdiag_use_amat) { 678 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 679 } else { 680 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 681 } 682 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 683 ilink = ilink->next; 684 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 685 if (jac->offdiag_use_amat) { 686 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 687 } else { 688 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 689 } 690 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 691 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 692 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 693 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 694 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 695 } 696 if (kspA != kspInner) { 697 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 698 } 699 if (kspUpper != kspA) { 700 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 701 } 702 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 703 } else { 704 const char *Dprefix; 705 char schurprefix[256], schurmatprefix[256]; 706 char schurtestoption[256]; 707 MatNullSpace sp; 708 PetscBool flg; 709 710 /* extract the A01 and A10 matrices */ 711 ilink = jac->head; 712 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 713 if (jac->offdiag_use_amat) { 714 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 715 } else { 716 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 717 } 718 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 719 ilink = ilink->next; 720 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 721 if (jac->offdiag_use_amat) { 722 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 723 } else { 724 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 725 } 726 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 727 728 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 729 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 730 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 731 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 732 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 733 /* Note that the inner KSP is NOT going to inherit this prefix, and if it did, it would be reset just below. Is that what we want? */ 734 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 735 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 736 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr); 737 if (sp) { 738 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 739 } 740 741 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 742 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 743 if (flg) { 744 DM dmInner; 745 KSP kspInner; 746 747 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 748 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 749 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 750 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 751 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 752 753 /* Set DM for new solver */ 754 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 755 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 756 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 757 } else { 758 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 759 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 760 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 761 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 762 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 763 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 764 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 765 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 766 } 767 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 768 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 769 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 770 771 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 772 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 773 if (flg) { 774 DM dmInner; 775 776 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 777 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 778 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr); 779 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 780 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 781 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 782 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 783 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 784 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 785 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 786 } else { 787 jac->kspupper = jac->head->ksp; 788 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 789 } 790 791 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 792 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 793 } 794 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 795 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr); 796 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 797 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 798 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 799 PC pcschur; 800 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 801 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 802 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 803 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 804 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 805 } 806 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 807 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 808 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 809 /* propogate DM */ 810 { 811 DM sdm; 812 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 813 if (sdm) { 814 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 815 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 816 } 817 } 818 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 819 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 820 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 821 } 822 823 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 824 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 825 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 826 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 827 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 828 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 829 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 830 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 831 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 832 } else { 833 /* set up the individual splits' PCs */ 834 i = 0; 835 ilink = jac->head; 836 while (ilink) { 837 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 838 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 839 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 840 i++; 841 ilink = ilink->next; 842 } 843 } 844 845 jac->suboptionsset = PETSC_TRUE; 846 PetscFunctionReturn(0); 847 } 848 849 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 850 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 851 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 852 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 853 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 854 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 855 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 856 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 857 858 #undef __FUNCT__ 859 #define __FUNCT__ "PCApply_FieldSplit_Schur" 860 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 861 { 862 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 863 PetscErrorCode ierr; 864 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 865 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 866 867 PetscFunctionBegin; 868 switch (jac->schurfactorization) { 869 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 870 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 871 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 873 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 875 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 876 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 877 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 878 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 879 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 880 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 881 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 882 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 883 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 884 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 885 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 886 break; 887 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 888 /* [A00 0; A10 S], suitable for left preconditioning */ 889 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 890 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 891 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 892 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 893 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 894 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 895 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 896 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 897 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 898 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 899 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 900 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 901 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 902 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 903 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 904 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 905 break; 906 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 907 /* [A00 A01; 0 S], suitable for right preconditioning */ 908 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 909 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 910 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 911 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 912 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 913 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 914 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 915 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 916 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 917 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 918 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 919 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 920 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 921 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 922 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 923 break; 924 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 925 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1], an exact solve if applied exactly, needs one extra solve with A */ 926 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 927 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 928 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 929 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 930 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 931 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 932 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 933 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 934 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 935 936 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 937 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 938 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 939 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 940 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 941 942 if (kspUpper == kspA) { 943 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 944 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 945 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 946 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 947 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 948 } else { 949 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 950 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 951 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 952 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 953 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 954 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 955 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 956 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 957 } 958 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 959 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 960 } 961 PetscFunctionReturn(0); 962 } 963 964 #undef __FUNCT__ 965 #define __FUNCT__ "PCApply_FieldSplit" 966 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 967 { 968 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 969 PetscErrorCode ierr; 970 PC_FieldSplitLink ilink = jac->head; 971 PetscInt cnt,bs; 972 KSPConvergedReason reason; 973 974 PetscFunctionBegin; 975 if (jac->type == PC_COMPOSITE_ADDITIVE) { 976 if (jac->defaultsplit) { 977 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 978 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 979 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 980 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 981 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 982 while (ilink) { 983 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 984 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 985 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 986 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 987 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 988 pc->failedreason = PC_SUBPC_ERROR; 989 } 990 ilink = ilink->next; 991 } 992 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 993 } else { 994 ierr = VecSet(y,0.0);CHKERRQ(ierr); 995 while (ilink) { 996 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 997 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 998 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 999 pc->failedreason = PC_SUBPC_ERROR; 1000 } 1001 ilink = ilink->next; 1002 } 1003 } 1004 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1005 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1006 /* solve on first block for first block variables */ 1007 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1008 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1009 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1010 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1011 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1012 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1013 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1014 pc->failedreason = PC_SUBPC_ERROR; 1015 } 1016 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1017 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1018 1019 /* compute the residual only onto second block variables using first block variables */ 1020 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 1021 ilink = ilink->next; 1022 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1023 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1024 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1025 1026 /* solve on second block variables */ 1027 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1028 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1029 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1030 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1031 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1032 pc->failedreason = PC_SUBPC_ERROR; 1033 } 1034 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1035 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1036 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1037 if (!jac->w1) { 1038 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1039 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1040 } 1041 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1042 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1043 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1044 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1045 pc->failedreason = PC_SUBPC_ERROR; 1046 } 1047 cnt = 1; 1048 while (ilink->next) { 1049 ilink = ilink->next; 1050 /* compute the residual only over the part of the vector needed */ 1051 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1052 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1053 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1054 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1055 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1056 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1057 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1058 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1059 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1060 pc->failedreason = PC_SUBPC_ERROR; 1061 } 1062 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1063 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1064 } 1065 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1066 cnt -= 2; 1067 while (ilink->previous) { 1068 ilink = ilink->previous; 1069 /* compute the residual only over the part of the vector needed */ 1070 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1071 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1072 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1073 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1074 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1075 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1076 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1077 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1078 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1079 pc->failedreason = PC_SUBPC_ERROR; 1080 } 1081 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1082 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1083 } 1084 } 1085 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1086 PetscFunctionReturn(0); 1087 } 1088 1089 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1090 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1091 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1092 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1093 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1094 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1095 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1096 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1097 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 1100 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1101 { 1102 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1103 PetscErrorCode ierr; 1104 PC_FieldSplitLink ilink = jac->head; 1105 PetscInt bs; 1106 KSPConvergedReason reason; 1107 1108 PetscFunctionBegin; 1109 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1110 if (jac->defaultsplit) { 1111 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1112 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of x vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1113 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1114 if (jac->bs > 0 && bs != jac->bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Blocksize of y vector %D does not match fieldsplit blocksize %D",bs,jac->bs); 1115 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1116 while (ilink) { 1117 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1118 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1119 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1120 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1121 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1122 pc->failedreason = PC_SUBPC_ERROR; 1123 } 1124 ilink = ilink->next; 1125 } 1126 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1127 } else { 1128 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1129 while (ilink) { 1130 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1131 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1132 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1133 pc->failedreason = PC_SUBPC_ERROR; 1134 } 1135 ilink = ilink->next; 1136 } 1137 } 1138 } else { 1139 if (!jac->w1) { 1140 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1141 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1142 } 1143 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1144 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1145 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1146 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1147 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1148 pc->failedreason = PC_SUBPC_ERROR; 1149 } 1150 while (ilink->next) { 1151 ilink = ilink->next; 1152 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1153 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1154 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1155 } 1156 while (ilink->previous) { 1157 ilink = ilink->previous; 1158 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1159 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1160 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1161 } 1162 } else { 1163 while (ilink->next) { /* get to last entry in linked list */ 1164 ilink = ilink->next; 1165 } 1166 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1167 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1168 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1169 pc->failedreason = PC_SUBPC_ERROR; 1170 } 1171 while (ilink->previous) { 1172 ilink = ilink->previous; 1173 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1174 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1175 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1176 } 1177 } 1178 } 1179 PetscFunctionReturn(0); 1180 } 1181 1182 #undef __FUNCT__ 1183 #define __FUNCT__ "PCReset_FieldSplit" 1184 static PetscErrorCode PCReset_FieldSplit(PC pc) 1185 { 1186 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1187 PetscErrorCode ierr; 1188 PC_FieldSplitLink ilink = jac->head,next; 1189 1190 PetscFunctionBegin; 1191 while (ilink) { 1192 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1193 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1194 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1195 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1196 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1197 if (!ilink->is_orig) { /* save the original IS */ 1198 ierr = PetscObjectReference((PetscObject)ilink->is);CHKERRQ(ierr); 1199 ilink->is_orig = ilink->is; 1200 } 1201 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1202 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1203 next = ilink->next; 1204 ilink = next; 1205 } 1206 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1207 if (jac->mat && jac->mat != jac->pmat) { 1208 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1209 } else if (jac->mat) { 1210 jac->mat = NULL; 1211 } 1212 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1213 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1214 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1215 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1216 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1217 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1218 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1219 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1220 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1221 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1222 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1223 jac->reset = PETSC_TRUE; 1224 jac->isrestrict = PETSC_FALSE; 1225 PetscFunctionReturn(0); 1226 } 1227 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "PCDestroy_FieldSplit" 1230 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1231 { 1232 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1233 PetscErrorCode ierr; 1234 PC_FieldSplitLink ilink = jac->head,next; 1235 1236 PetscFunctionBegin; 1237 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1238 while (ilink) { 1239 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1240 ierr = ISDestroy(&ilink->is_orig);CHKERRQ(ierr); 1241 next = ilink->next; 1242 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1243 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1244 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1245 ierr = PetscFree(ilink);CHKERRQ(ierr); 1246 ilink = next; 1247 } 1248 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1249 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1250 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1251 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1252 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1253 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1254 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1255 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1256 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1257 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1258 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr); 1259 PetscFunctionReturn(0); 1260 } 1261 1262 #undef __FUNCT__ 1263 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1264 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1265 { 1266 PetscErrorCode ierr; 1267 PetscInt bs; 1268 PetscBool flg,stokes = PETSC_FALSE; 1269 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1270 PCCompositeType ctype; 1271 1272 PetscFunctionBegin; 1273 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1274 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1275 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1276 if (flg) { 1277 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1278 } 1279 jac->diag_use_amat = pc->useAmat; 1280 ierr = PetscOptionsBool("-pc_fieldsplit_diag_use_amat","Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat",jac->diag_use_amat,&jac->diag_use_amat,NULL);CHKERRQ(ierr); 1281 jac->offdiag_use_amat = pc->useAmat; 1282 ierr = PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat","Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat",jac->offdiag_use_amat,&jac->offdiag_use_amat,NULL);CHKERRQ(ierr); 1283 /* FIXME: No programmatic equivalent to the following. */ 1284 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1285 if (stokes) { 1286 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1287 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1288 } 1289 1290 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1291 if (flg) { 1292 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1293 } 1294 /* Only setup fields once */ 1295 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1296 /* only allow user to set fields from command line if bs is already known. 1297 otherwise user can set them in PCFieldSplitSetDefaults() */ 1298 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1299 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1300 } 1301 if (jac->type == PC_COMPOSITE_SCHUR) { 1302 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1303 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1304 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_fact_type","Which off-diagonal parts of the block factorization to use","PCFieldSplitSetSchurFactType",PCFieldSplitSchurFactTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,NULL);CHKERRQ(ierr); 1305 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1306 } 1307 ierr = PetscOptionsTail();CHKERRQ(ierr); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 /*------------------------------------------------------------------------------------*/ 1312 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1315 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1316 { 1317 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1318 PetscErrorCode ierr; 1319 PC_FieldSplitLink ilink,next = jac->head; 1320 char prefix[128]; 1321 PetscInt i; 1322 1323 PetscFunctionBegin; 1324 if (jac->splitdefined) { 1325 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1326 PetscFunctionReturn(0); 1327 } 1328 for (i=0; i<n; i++) { 1329 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1330 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1331 } 1332 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1333 if (splitname) { 1334 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1335 } else { 1336 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1337 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1338 } 1339 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1340 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1341 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1342 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1343 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1344 1345 ilink->nfields = n; 1346 ilink->next = NULL; 1347 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1348 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1349 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1350 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1351 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1352 1353 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1354 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1355 1356 if (!next) { 1357 jac->head = ilink; 1358 ilink->previous = NULL; 1359 } else { 1360 while (next->next) { 1361 next = next->next; 1362 } 1363 next->next = ilink; 1364 ilink->previous = next; 1365 } 1366 jac->nsplits++; 1367 PetscFunctionReturn(0); 1368 } 1369 1370 #undef __FUNCT__ 1371 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1372 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1373 { 1374 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1375 PetscErrorCode ierr; 1376 1377 PetscFunctionBegin; 1378 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1379 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1380 1381 (*subksp)[1] = jac->kspschur; 1382 if (n) *n = jac->nsplits; 1383 PetscFunctionReturn(0); 1384 } 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1388 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1389 { 1390 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1391 PetscErrorCode ierr; 1392 PetscInt cnt = 0; 1393 PC_FieldSplitLink ilink = jac->head; 1394 1395 PetscFunctionBegin; 1396 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1397 while (ilink) { 1398 (*subksp)[cnt++] = ilink->ksp; 1399 ilink = ilink->next; 1400 } 1401 if (cnt != jac->nsplits) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number of splits in linked list %D does not match number in object %D",cnt,jac->nsplits); 1402 if (n) *n = jac->nsplits; 1403 PetscFunctionReturn(0); 1404 } 1405 1406 #undef __FUNCT__ 1407 #define __FUNCT__ "PCFieldSplitRestrictIS" 1408 /*@C 1409 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1410 1411 Input Parameters: 1412 + pc - the preconditioner context 1413 + is - the index set that defines the indices to which the fieldsplit is to be restricted 1414 1415 Level: advanced 1416 1417 @*/ 1418 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1419 { 1420 PetscErrorCode ierr; 1421 1422 PetscFunctionBegin; 1423 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1424 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1425 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr); 1426 PetscFunctionReturn(0); 1427 } 1428 1429 1430 #undef __FUNCT__ 1431 #define __FUNCT__ "PCFieldSplitRestrictIS_FieldSplit" 1432 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1433 { 1434 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1435 PetscErrorCode ierr; 1436 PC_FieldSplitLink ilink = jac->head, next; 1437 PetscInt localsize,size,sizez,i; 1438 const PetscInt *ind, *indz; 1439 PetscInt *indc, *indcz; 1440 PetscBool flg; 1441 1442 PetscFunctionBegin; 1443 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr); 1444 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr); 1445 size -= localsize; 1446 while(ilink) { 1447 IS isrl,isr; 1448 PC subpc; 1449 if (jac->reset) { 1450 ierr = ISEmbed(ilink->is_orig, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1451 } else { 1452 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1453 } 1454 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr); 1455 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr); 1456 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr); 1457 ierr = PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1458 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr); 1459 ierr = ISDestroy(&isrl);CHKERRQ(ierr); 1460 for (i=0; i<localsize; i++) *(indc+i) += size; 1461 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr); 1462 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1463 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1464 ilink->is = isr; 1465 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1466 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1467 ilink->is_col = isr; 1468 ierr = ISDestroy(&isr);CHKERRQ(ierr); 1469 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr); 1470 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1471 if(flg) { 1472 IS iszl,isz; 1473 MPI_Comm comm; 1474 if (jac->reset) { 1475 ierr = ISGetLocalSize(ilink->is_orig,&localsize);CHKERRQ(ierr); 1476 comm = PetscObjectComm((PetscObject)ilink->is_orig); 1477 ierr = ISEmbed(isy, ilink->is_orig, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1478 } else { 1479 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr); 1480 comm = PetscObjectComm((PetscObject)ilink->is); 1481 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1482 } 1483 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 1484 sizez -= localsize; 1485 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr); 1486 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr); 1487 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr); 1488 ierr = PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1489 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr); 1490 ierr = ISDestroy(&iszl);CHKERRQ(ierr); 1491 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1492 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr); 1493 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr); 1494 ierr = ISDestroy(&isz);CHKERRQ(ierr); 1495 } 1496 next = ilink->next; 1497 ilink = next; 1498 } 1499 jac->isrestrict = PETSC_TRUE; 1500 PetscFunctionReturn(0); 1501 } 1502 1503 #undef __FUNCT__ 1504 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1505 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1506 { 1507 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1508 PetscErrorCode ierr; 1509 PC_FieldSplitLink ilink, next = jac->head; 1510 char prefix[128]; 1511 1512 PetscFunctionBegin; 1513 if (jac->splitdefined) { 1514 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1515 PetscFunctionReturn(0); 1516 } 1517 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1518 if (splitname) { 1519 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1520 } else { 1521 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1522 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1523 } 1524 ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + jac->nsplits : KSP_Solve_FS_0 + 4; /* Any split great than 4 gets logged in the 4th split */ 1525 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1526 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1527 ilink->is = is; 1528 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1529 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1530 ilink->is_col = is; 1531 ilink->next = NULL; 1532 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1533 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1534 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1535 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1536 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1537 1538 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1539 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1540 1541 if (!next) { 1542 jac->head = ilink; 1543 ilink->previous = NULL; 1544 } else { 1545 while (next->next) { 1546 next = next->next; 1547 } 1548 next->next = ilink; 1549 ilink->previous = next; 1550 } 1551 jac->nsplits++; 1552 PetscFunctionReturn(0); 1553 } 1554 1555 #undef __FUNCT__ 1556 #define __FUNCT__ "PCFieldSplitSetFields" 1557 /*@ 1558 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1559 1560 Logically Collective on PC 1561 1562 Input Parameters: 1563 + pc - the preconditioner context 1564 . splitname - name of this split, if NULL the number of the split is used 1565 . n - the number of fields in this split 1566 - fields - the fields in this split 1567 1568 Level: intermediate 1569 1570 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1571 1572 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1573 size is three then one can define a field as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean 1574 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1575 where the numbered entries indicate what is in the field. 1576 1577 This function is called once per split (it creates a new split each time). Solve options 1578 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1579 1580 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1581 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1582 available when this routine is called. 1583 1584 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1585 1586 @*/ 1587 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1588 { 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1593 PetscValidCharPointer(splitname,2); 1594 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1595 PetscValidIntPointer(fields,3); 1596 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1597 PetscFunctionReturn(0); 1598 } 1599 1600 #undef __FUNCT__ 1601 #define __FUNCT__ "PCFieldSplitSetDiagUseAmat" 1602 /*@ 1603 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1604 1605 Logically Collective on PC 1606 1607 Input Parameters: 1608 + pc - the preconditioner object 1609 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1610 1611 Options Database: 1612 . -pc_fieldsplit_diag_use_amat 1613 1614 Level: intermediate 1615 1616 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1617 1618 @*/ 1619 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1620 { 1621 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1622 PetscBool isfs; 1623 PetscErrorCode ierr; 1624 1625 PetscFunctionBegin; 1626 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1627 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1628 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1629 jac->diag_use_amat = flg; 1630 PetscFunctionReturn(0); 1631 } 1632 1633 #undef __FUNCT__ 1634 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1635 /*@ 1636 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1637 1638 Logically Collective on PC 1639 1640 Input Parameters: 1641 . pc - the preconditioner object 1642 1643 Output Parameters: 1644 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1645 1646 1647 Level: intermediate 1648 1649 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1650 1651 @*/ 1652 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1653 { 1654 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1655 PetscBool isfs; 1656 PetscErrorCode ierr; 1657 1658 PetscFunctionBegin; 1659 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1660 PetscValidPointer(flg,2); 1661 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1662 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1663 *flg = jac->diag_use_amat; 1664 PetscFunctionReturn(0); 1665 } 1666 1667 #undef __FUNCT__ 1668 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1669 /*@ 1670 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1671 1672 Logically Collective on PC 1673 1674 Input Parameters: 1675 + pc - the preconditioner object 1676 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1677 1678 Options Database: 1679 . -pc_fieldsplit_off_diag_use_amat 1680 1681 Level: intermediate 1682 1683 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1684 1685 @*/ 1686 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1687 { 1688 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1689 PetscBool isfs; 1690 PetscErrorCode ierr; 1691 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1694 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1695 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1696 jac->offdiag_use_amat = flg; 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1702 /*@ 1703 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1704 1705 Logically Collective on PC 1706 1707 Input Parameters: 1708 . pc - the preconditioner object 1709 1710 Output Parameters: 1711 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1712 1713 1714 Level: intermediate 1715 1716 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1717 1718 @*/ 1719 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1720 { 1721 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1722 PetscBool isfs; 1723 PetscErrorCode ierr; 1724 1725 PetscFunctionBegin; 1726 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1727 PetscValidPointer(flg,2); 1728 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1729 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1730 *flg = jac->offdiag_use_amat; 1731 PetscFunctionReturn(0); 1732 } 1733 1734 1735 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "PCFieldSplitSetIS" 1738 /*@C 1739 PCFieldSplitSetIS - Sets the exact elements for field 1740 1741 Logically Collective on PC 1742 1743 Input Parameters: 1744 + pc - the preconditioner context 1745 . splitname - name of this split, if NULL the number of the split is used 1746 - is - the index set that defines the vector elements in this field 1747 1748 1749 Notes: 1750 Use PCFieldSplitSetFields(), for fields defined by strided types. 1751 1752 This function is called once per split (it creates a new split each time). Solve options 1753 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1754 1755 Level: intermediate 1756 1757 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1758 1759 @*/ 1760 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1761 { 1762 PetscErrorCode ierr; 1763 1764 PetscFunctionBegin; 1765 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1766 if (splitname) PetscValidCharPointer(splitname,2); 1767 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1768 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1769 PetscFunctionReturn(0); 1770 } 1771 1772 #undef __FUNCT__ 1773 #define __FUNCT__ "PCFieldSplitGetIS" 1774 /*@ 1775 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1776 1777 Logically Collective on PC 1778 1779 Input Parameters: 1780 + pc - the preconditioner context 1781 - splitname - name of this split 1782 1783 Output Parameter: 1784 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1785 1786 Level: intermediate 1787 1788 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1789 1790 @*/ 1791 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1792 { 1793 PetscErrorCode ierr; 1794 1795 PetscFunctionBegin; 1796 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1797 PetscValidCharPointer(splitname,2); 1798 PetscValidPointer(is,3); 1799 { 1800 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1801 PC_FieldSplitLink ilink = jac->head; 1802 PetscBool found; 1803 1804 *is = NULL; 1805 while (ilink) { 1806 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1807 if (found) { 1808 *is = ilink->is; 1809 break; 1810 } 1811 ilink = ilink->next; 1812 } 1813 } 1814 PetscFunctionReturn(0); 1815 } 1816 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1819 /*@ 1820 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1821 fieldsplit preconditioner. If not set the matrix block size is used. 1822 1823 Logically Collective on PC 1824 1825 Input Parameters: 1826 + pc - the preconditioner context 1827 - bs - the block size 1828 1829 Level: intermediate 1830 1831 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1832 1833 @*/ 1834 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1835 { 1836 PetscErrorCode ierr; 1837 1838 PetscFunctionBegin; 1839 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1840 PetscValidLogicalCollectiveInt(pc,bs,2); 1841 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 #undef __FUNCT__ 1846 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1847 /*@C 1848 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1849 1850 Collective on KSP 1851 1852 Input Parameter: 1853 . pc - the preconditioner context 1854 1855 Output Parameters: 1856 + n - the number of splits 1857 - pc - the array of KSP contexts 1858 1859 Note: 1860 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1861 (not the KSP just the array that contains them). 1862 1863 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1864 1865 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1866 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1867 KSP array must be. 1868 1869 1870 Level: advanced 1871 1872 .seealso: PCFIELDSPLIT 1873 @*/ 1874 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1875 { 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1880 if (n) PetscValidIntPointer(n,2); 1881 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1882 PetscFunctionReturn(0); 1883 } 1884 1885 #undef __FUNCT__ 1886 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1887 /*@ 1888 PCFieldSplitSetSchurPre - Indicates what operator is used to construct the preconditioner for the Schur complement. 1889 A11 matrix. Otherwise no preconditioner is used. 1890 1891 Collective on PC 1892 1893 Input Parameters: 1894 + pc - the preconditioner context 1895 . ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11 (default), PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_SCHUR_PRE_USER 1896 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 1897 - userpre - matrix to use for preconditioning, or NULL 1898 1899 Options Database: 1900 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 1901 1902 Notes: 1903 $ If ptype is 1904 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1905 $ matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 1906 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1907 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1908 $ preconditioner 1909 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1910 $ to this function). 1911 $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1912 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1913 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 1914 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive) 1915 $ useful mostly as a test that the Schur complement approach can work for your problem 1916 1917 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1918 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1919 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1920 1921 Level: intermediate 1922 1923 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1924 MatSchurComplementSetAinvType(), PCLSC 1925 1926 @*/ 1927 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1928 { 1929 PetscErrorCode ierr; 1930 1931 PetscFunctionBegin; 1932 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1933 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1934 PetscFunctionReturn(0); 1935 } 1936 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1937 1938 #undef __FUNCT__ 1939 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1940 /*@ 1941 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1942 preconditioned. See PCFieldSplitSetSchurPre() for details. 1943 1944 Logically Collective on PC 1945 1946 Input Parameters: 1947 . pc - the preconditioner context 1948 1949 Output Parameters: 1950 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1951 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1952 1953 Level: intermediate 1954 1955 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1956 1957 @*/ 1958 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1959 { 1960 PetscErrorCode ierr; 1961 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1964 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1965 PetscFunctionReturn(0); 1966 } 1967 1968 #undef __FUNCT__ 1969 #define __FUNCT__ "PCFieldSplitSchurGetS" 1970 /*@ 1971 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1972 1973 Not collective 1974 1975 Input Parameter: 1976 . pc - the preconditioner context 1977 1978 Output Parameter: 1979 . S - the Schur complement matrix 1980 1981 Notes: 1982 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1983 1984 Level: advanced 1985 1986 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1987 1988 @*/ 1989 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1990 { 1991 PetscErrorCode ierr; 1992 const char* t; 1993 PetscBool isfs; 1994 PC_FieldSplit *jac; 1995 1996 PetscFunctionBegin; 1997 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1998 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1999 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2000 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2001 jac = (PC_FieldSplit*)pc->data; 2002 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2003 if (S) *S = jac->schur; 2004 PetscFunctionReturn(0); 2005 } 2006 2007 #undef __FUNCT__ 2008 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 2009 /*@ 2010 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 2011 2012 Not collective 2013 2014 Input Parameters: 2015 + pc - the preconditioner context 2016 . S - the Schur complement matrix 2017 2018 Level: advanced 2019 2020 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 2021 2022 @*/ 2023 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2024 { 2025 PetscErrorCode ierr; 2026 const char* t; 2027 PetscBool isfs; 2028 PC_FieldSplit *jac; 2029 2030 PetscFunctionBegin; 2031 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2032 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 2033 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2034 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2035 jac = (PC_FieldSplit*)pc->data; 2036 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2037 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2038 PetscFunctionReturn(0); 2039 } 2040 2041 2042 #undef __FUNCT__ 2043 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 2044 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2045 { 2046 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2047 PetscErrorCode ierr; 2048 2049 PetscFunctionBegin; 2050 jac->schurpre = ptype; 2051 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2052 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 2053 jac->schur_user = pre; 2054 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 2055 } 2056 PetscFunctionReturn(0); 2057 } 2058 2059 #undef __FUNCT__ 2060 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 2061 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2062 { 2063 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2064 2065 PetscFunctionBegin; 2066 *ptype = jac->schurpre; 2067 *pre = jac->schur_user; 2068 PetscFunctionReturn(0); 2069 } 2070 2071 #undef __FUNCT__ 2072 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 2073 /*@ 2074 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 2075 2076 Collective on PC 2077 2078 Input Parameters: 2079 + pc - the preconditioner context 2080 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2081 2082 Options Database: 2083 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 2084 2085 2086 Level: intermediate 2087 2088 Notes: 2089 The FULL factorization is 2090 2091 $ (A B) = (1 0) (A 0) (1 Ainv*B) 2092 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 2093 2094 where S = D - C*Ainv*B. In practice, the full factorization is applied via block triangular solves with the grouping L*(D*U). UPPER uses D*U, LOWER uses L*D, 2095 and DIAG is the diagonal part with the sign of S flipped (because this makes the preconditioner positive definite for many formulations, thus allowing the use of KSPMINRES). 2096 2097 If applied exactly, FULL factorization is a direct solver. The preconditioned operator with LOWER or UPPER has all eigenvalues equal to 1 and minimal polynomial 2098 of degree 2, so KSPGMRES converges in 2 iterations. If the iteration count is very low, consider using KSPFGMRES or KSPGCR which can use one less preconditioner 2099 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. With DIAG, 2100 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2101 2102 For symmetric problems in which A is positive definite and S is negative definite, DIAG can be used with KSPMINRES. Note that a flexible method like KSPFGMRES 2103 or KSPGCR must be used if the fieldsplit preconditioner is nonlinear (e.g. a few iterations of a Krylov method is used inside a split). 2104 2105 References: 2106 + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2107 - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2108 2109 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 2110 @*/ 2111 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2112 { 2113 PetscErrorCode ierr; 2114 2115 PetscFunctionBegin; 2116 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2117 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 2118 PetscFunctionReturn(0); 2119 } 2120 2121 #undef __FUNCT__ 2122 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 2123 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2124 { 2125 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2126 2127 PetscFunctionBegin; 2128 jac->schurfactorization = ftype; 2129 PetscFunctionReturn(0); 2130 } 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 2134 /*@C 2135 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2136 2137 Collective on KSP 2138 2139 Input Parameter: 2140 . pc - the preconditioner context 2141 2142 Output Parameters: 2143 + A00 - the (0,0) block 2144 . A01 - the (0,1) block 2145 . A10 - the (1,0) block 2146 - A11 - the (1,1) block 2147 2148 Level: advanced 2149 2150 .seealso: PCFIELDSPLIT 2151 @*/ 2152 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2153 { 2154 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2155 2156 PetscFunctionBegin; 2157 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2158 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2159 if (A00) *A00 = jac->pmat[0]; 2160 if (A01) *A01 = jac->B; 2161 if (A10) *A10 = jac->C; 2162 if (A11) *A11 = jac->pmat[1]; 2163 PetscFunctionReturn(0); 2164 } 2165 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 2168 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2169 { 2170 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2171 PetscErrorCode ierr; 2172 2173 PetscFunctionBegin; 2174 jac->type = type; 2175 if (type == PC_COMPOSITE_SCHUR) { 2176 pc->ops->apply = PCApply_FieldSplit_Schur; 2177 pc->ops->view = PCView_FieldSplit_Schur; 2178 2179 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2180 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2181 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2182 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2183 2184 } else { 2185 pc->ops->apply = PCApply_FieldSplit; 2186 pc->ops->view = PCView_FieldSplit; 2187 2188 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2189 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2190 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2191 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2192 } 2193 PetscFunctionReturn(0); 2194 } 2195 2196 #undef __FUNCT__ 2197 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 2198 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2199 { 2200 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2201 2202 PetscFunctionBegin; 2203 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2204 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 2205 jac->bs = bs; 2206 PetscFunctionReturn(0); 2207 } 2208 2209 #undef __FUNCT__ 2210 #define __FUNCT__ "PCFieldSplitSetType" 2211 /*@ 2212 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2213 2214 Collective on PC 2215 2216 Input Parameter: 2217 . pc - the preconditioner context 2218 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2219 2220 Options Database Key: 2221 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2222 2223 Level: Intermediate 2224 2225 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2226 2227 .seealso: PCCompositeSetType() 2228 2229 @*/ 2230 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2231 { 2232 PetscErrorCode ierr; 2233 2234 PetscFunctionBegin; 2235 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2236 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2237 PetscFunctionReturn(0); 2238 } 2239 2240 #undef __FUNCT__ 2241 #define __FUNCT__ "PCFieldSplitGetType" 2242 /*@ 2243 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2244 2245 Not collective 2246 2247 Input Parameter: 2248 . pc - the preconditioner context 2249 2250 Output Parameter: 2251 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2252 2253 Level: Intermediate 2254 2255 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2256 .seealso: PCCompositeSetType() 2257 @*/ 2258 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2259 { 2260 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2261 2262 PetscFunctionBegin; 2263 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2264 PetscValidIntPointer(type,2); 2265 *type = jac->type; 2266 PetscFunctionReturn(0); 2267 } 2268 2269 #undef __FUNCT__ 2270 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2271 /*@ 2272 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2273 2274 Logically Collective 2275 2276 Input Parameters: 2277 + pc - the preconditioner context 2278 - flg - boolean indicating whether to use field splits defined by the DM 2279 2280 Options Database Key: 2281 . -pc_fieldsplit_dm_splits 2282 2283 Level: Intermediate 2284 2285 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2286 2287 .seealso: PCFieldSplitGetDMSplits() 2288 2289 @*/ 2290 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2291 { 2292 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2293 PetscBool isfs; 2294 PetscErrorCode ierr; 2295 2296 PetscFunctionBegin; 2297 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2298 PetscValidLogicalCollectiveBool(pc,flg,2); 2299 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2300 if (isfs) { 2301 jac->dm_splits = flg; 2302 } 2303 PetscFunctionReturn(0); 2304 } 2305 2306 2307 #undef __FUNCT__ 2308 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2309 /*@ 2310 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2311 2312 Logically Collective 2313 2314 Input Parameter: 2315 . pc - the preconditioner context 2316 2317 Output Parameter: 2318 . flg - boolean indicating whether to use field splits defined by the DM 2319 2320 Level: Intermediate 2321 2322 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2323 2324 .seealso: PCFieldSplitSetDMSplits() 2325 2326 @*/ 2327 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2328 { 2329 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2330 PetscBool isfs; 2331 PetscErrorCode ierr; 2332 2333 PetscFunctionBegin; 2334 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2335 PetscValidPointer(flg,2); 2336 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2337 if (isfs) { 2338 if(flg) *flg = jac->dm_splits; 2339 } 2340 PetscFunctionReturn(0); 2341 } 2342 2343 /* -------------------------------------------------------------------------------------*/ 2344 /*MC 2345 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2346 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2347 2348 To set options on the solvers for each block append -fieldsplit_ to all the PC 2349 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2350 2351 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2352 and set the options directly on the resulting KSP object 2353 2354 Level: intermediate 2355 2356 Options Database Keys: 2357 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2358 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2359 been supplied explicitly by -pc_fieldsplit_%d_fields 2360 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2361 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2362 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2363 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2364 2365 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2366 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2367 2368 Notes: 2369 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2370 to define a field by an arbitrary collection of entries. 2371 2372 If no fields are set the default is used. The fields are defined by entries strided by bs, 2373 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2374 if this is not called the block size defaults to the blocksize of the second matrix passed 2375 to KSPSetOperators()/PCSetOperators(). 2376 2377 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2378 $ ( A10 A11 ) 2379 $ the preconditioner using full factorization is 2380 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2381 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2382 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2383 $ S = A11 - A10 ksp(A00) A01 2384 which is usually dense and not stored explicitly. The action of ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 2385 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2386 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2387 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2388 2389 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2390 diag gives 2391 $ ( inv(A00) 0 ) 2392 $ ( 0 -ksp(S) ) 2393 note that slightly counter intuitively there is a negative in front of the ksp(S) so that the preconditioner is positive definite. The lower factorization is the inverse of 2394 $ ( A00 0 ) 2395 $ ( A10 S ) 2396 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2397 $ ( A00 A01 ) 2398 $ ( 0 S ) 2399 where again the inverses of A00 and S are applied using KSPs. 2400 2401 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2402 is used automatically for a second block. 2403 2404 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2405 Generally it should be used with the AIJ format. 2406 2407 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2408 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2409 inside a smoother resulting in "Distributive Smoothers". 2410 2411 Concepts: physics based preconditioners, block preconditioners 2412 2413 There is a nice discussion of block preconditioners in 2414 2415 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2416 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2417 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2418 2419 The Constrained Pressure Preconditioner (CPR) does not appear to be currently implementable directly with PCFIELDSPLIT. CPR solves first the Schur complemented pressure equation, updates the 2420 residual on all variables and then applies a simple ILU like preconditioner on all the variables. So it is very much like the full Schur complement with selfp representing the Schur complement but instead 2421 of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables. 2422 2423 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2424 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2425 MatSchurComplementSetAinvType() 2426 M*/ 2427 2428 #undef __FUNCT__ 2429 #define __FUNCT__ "PCCreate_FieldSplit" 2430 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2431 { 2432 PetscErrorCode ierr; 2433 PC_FieldSplit *jac; 2434 2435 PetscFunctionBegin; 2436 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2437 2438 jac->bs = -1; 2439 jac->nsplits = 0; 2440 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2441 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2442 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2443 jac->dm_splits = PETSC_TRUE; 2444 2445 pc->data = (void*)jac; 2446 2447 pc->ops->apply = PCApply_FieldSplit; 2448 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2449 pc->ops->setup = PCSetUp_FieldSplit; 2450 pc->ops->reset = PCReset_FieldSplit; 2451 pc->ops->destroy = PCDestroy_FieldSplit; 2452 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2453 pc->ops->view = PCView_FieldSplit; 2454 pc->ops->applyrichardson = 0; 2455 2456 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2457 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2458 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2459 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2460 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2461 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr); 2462 PetscFunctionReturn(0); 2463 } 2464 2465 2466 2467