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 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_INITIAL_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 624 } else { 625 ierr = MatGetSubMatrix(pc->mat,ilink->next->is,ilink->is,MAT_REUSE_MATRIX,&jac->Afield[1]);CHKERRQ(ierr); 626 } 627 } else { 628 if (!jac->Afield) { 629 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 630 for (i=0; i<nsplit; i++) { 631 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 632 ilink = ilink->next; 633 } 634 } else { 635 for (i=0; i<nsplit; i++) { 636 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 637 ilink = ilink->next; 638 } 639 } 640 } 641 } 642 643 if (jac->type == PC_COMPOSITE_SCHUR) { 644 IS ccis; 645 PetscInt rstart,rend; 646 char lscname[256]; 647 PetscObject LSC_L; 648 649 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 650 651 /* When extracting off-diagonal submatrices, we take complements from this range */ 652 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 653 654 /* need to handle case when one is resetting up the preconditioner */ 655 if (jac->schur) { 656 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 657 658 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 659 ilink = jac->head; 660 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 661 if (jac->offdiag_use_amat) { 662 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 663 } else { 664 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 665 } 666 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 667 ilink = ilink->next; 668 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 669 if (jac->offdiag_use_amat) { 670 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 671 } else { 672 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 673 } 674 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 675 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 676 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 677 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 678 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 679 } 680 if (kspA != kspInner) { 681 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 682 } 683 if (kspUpper != kspA) { 684 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 685 } 686 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 687 } else { 688 const char *Dprefix; 689 char schurprefix[256], schurmatprefix[256]; 690 char schurtestoption[256]; 691 MatNullSpace sp; 692 PetscBool flg; 693 694 /* extract the A01 and A10 matrices */ 695 ilink = jac->head; 696 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 697 if (jac->offdiag_use_amat) { 698 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 699 } else { 700 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 701 } 702 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 703 ilink = ilink->next; 704 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 705 if (jac->offdiag_use_amat) { 706 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 707 } else { 708 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 709 } 710 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 711 712 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 713 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 714 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 715 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 716 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 717 /* 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? */ 718 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 719 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 720 ierr = MatGetNullSpace(jac->mat[1], &sp);CHKERRQ(ierr); 721 if (sp) { 722 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 723 } 724 725 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 726 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 727 if (flg) { 728 DM dmInner; 729 KSP kspInner; 730 731 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 732 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 733 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 734 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 735 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 736 737 /* Set DM for new solver */ 738 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 739 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 740 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 741 } else { 742 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 743 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 744 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 745 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 746 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 747 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 748 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 749 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 750 } 751 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 752 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 753 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 754 755 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 756 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options,((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 757 if (flg) { 758 DM dmInner; 759 760 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 761 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 762 ierr = KSPSetErrorIfNotConverged(jac->kspupper,pc->erroriffailure);CHKERRQ(ierr); 763 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 764 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 765 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 766 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 767 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 768 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0]);CHKERRQ(ierr); 769 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 770 } else { 771 jac->kspupper = jac->head->ksp; 772 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 773 } 774 775 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 776 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 777 } 778 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 779 ierr = KSPSetErrorIfNotConverged(jac->kspschur,pc->erroriffailure);CHKERRQ(ierr); 780 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 781 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 782 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 783 PC pcschur; 784 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 785 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 786 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 787 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 788 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 789 } 790 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac));CHKERRQ(ierr); 791 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 792 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 793 /* propogate DM */ 794 { 795 DM sdm; 796 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 797 if (sdm) { 798 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 799 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 800 } 801 } 802 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 803 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 804 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 805 } 806 807 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 808 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 809 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 810 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 811 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 812 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 813 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 814 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 815 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 816 } else { 817 /* set up the individual splits' PCs */ 818 i = 0; 819 ilink = jac->head; 820 while (ilink) { 821 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i]);CHKERRQ(ierr); 822 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 823 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 824 i++; 825 ilink = ilink->next; 826 } 827 } 828 829 jac->suboptionsset = PETSC_TRUE; 830 PetscFunctionReturn(0); 831 } 832 833 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 834 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 835 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 836 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 837 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 838 PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL) ||\ 839 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 840 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "PCApply_FieldSplit_Schur" 844 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 845 { 846 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 847 PetscErrorCode ierr; 848 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 849 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 850 851 PetscFunctionBegin; 852 switch (jac->schurfactorization) { 853 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 854 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 855 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 856 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 857 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 858 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 859 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 860 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 861 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 862 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 863 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 864 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 865 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 866 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 867 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 868 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 869 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 870 break; 871 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 872 /* [A00 0; A10 S], suitable for left preconditioning */ 873 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 874 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 875 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 876 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 877 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 878 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 879 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 880 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 881 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 882 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 883 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 884 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 885 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 886 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 887 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 888 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 889 break; 890 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 891 /* [A00 A01; 0 S], suitable for right preconditioning */ 892 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 893 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 894 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 895 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 896 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); 897 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 898 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 899 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 900 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 901 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 902 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 903 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 904 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 905 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 906 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 907 break; 908 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 909 /* [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 */ 910 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 911 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 912 ierr = PetscLogEventBegin(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 913 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 914 ierr = PetscLogEventEnd(KSP_Solve_FS_L,kspLower,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 915 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 916 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 917 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 918 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 919 920 ierr = PetscLogEventBegin(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 921 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 922 ierr = PetscLogEventEnd(KSP_Solve_FS_S,jac->kspschur,ilinkD->x,ilinkD->y,NULL);CHKERRQ(ierr); 923 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 924 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 925 926 if (kspUpper == kspA) { 927 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 928 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 929 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 930 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 931 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 932 } else { 933 ierr = PetscLogEventBegin(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 934 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 935 ierr = PetscLogEventEnd(ilinkA->event,kspA,ilinkA->x,ilinkA->y,NULL);CHKERRQ(ierr); 936 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 937 ierr = PetscLogEventBegin(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 938 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 939 ierr = PetscLogEventEnd(KSP_Solve_FS_U,kspUpper,ilinkA->x,ilinkA->z,NULL);CHKERRQ(ierr); 940 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 941 } 942 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 943 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 944 } 945 PetscFunctionReturn(0); 946 } 947 948 #undef __FUNCT__ 949 #define __FUNCT__ "PCApply_FieldSplit" 950 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 951 { 952 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 953 PetscErrorCode ierr; 954 PC_FieldSplitLink ilink = jac->head; 955 PetscInt cnt,bs; 956 KSPConvergedReason reason; 957 958 PetscFunctionBegin; 959 if (jac->type == PC_COMPOSITE_ADDITIVE) { 960 if (jac->defaultsplit) { 961 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 962 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); 963 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 964 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); 965 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 966 while (ilink) { 967 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 968 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 969 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 970 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 971 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 972 pc->failedreason = PC_SUBPC_ERROR; 973 } 974 ilink = ilink->next; 975 } 976 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 977 } else { 978 ierr = VecSet(y,0.0);CHKERRQ(ierr); 979 while (ilink) { 980 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 981 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 982 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 983 pc->failedreason = PC_SUBPC_ERROR; 984 } 985 ilink = ilink->next; 986 } 987 } 988 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 989 ierr = VecSet(y,0.0);CHKERRQ(ierr); 990 /* solve on first block for first block variables */ 991 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 992 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 993 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 994 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 995 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 996 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 997 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 998 pc->failedreason = PC_SUBPC_ERROR; 999 } 1000 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1001 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1002 1003 /* compute the residual only onto second block variables using first block variables */ 1004 ierr = MatMult(jac->Afield[1],ilink->y,ilink->next->x);CHKERRQ(ierr); 1005 ilink = ilink->next; 1006 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1007 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1008 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1009 1010 /* solve on second block variables */ 1011 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1012 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1013 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1014 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1015 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1016 pc->failedreason = PC_SUBPC_ERROR; 1017 } 1018 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1019 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1020 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1021 if (!jac->w1) { 1022 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1023 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1024 } 1025 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1026 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 1027 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1028 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1029 pc->failedreason = PC_SUBPC_ERROR; 1030 } 1031 cnt = 1; 1032 while (ilink->next) { 1033 ilink = ilink->next; 1034 /* compute the residual only over the part of the vector needed */ 1035 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 1036 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1037 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1038 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1039 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1040 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1041 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1042 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1043 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1044 pc->failedreason = PC_SUBPC_ERROR; 1045 } 1046 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1047 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1048 } 1049 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1050 cnt -= 2; 1051 while (ilink->previous) { 1052 ilink = ilink->previous; 1053 /* compute the residual only over the part of the vector needed */ 1054 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 1055 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 1056 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1057 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1058 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1059 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1060 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1061 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1062 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1063 pc->failedreason = PC_SUBPC_ERROR; 1064 } 1065 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1066 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1067 } 1068 } 1069 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 1074 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1075 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 1076 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1077 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 1078 PetscLogEventBegin(ilink->event,ilink->ksp,ilink->y,ilink->x,NULL) || \ 1079 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 1080 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 1081 1082 #undef __FUNCT__ 1083 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 1084 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 1085 { 1086 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1087 PetscErrorCode ierr; 1088 PC_FieldSplitLink ilink = jac->head; 1089 PetscInt bs; 1090 KSPConvergedReason reason; 1091 1092 PetscFunctionBegin; 1093 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1094 if (jac->defaultsplit) { 1095 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 1096 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); 1097 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 1098 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); 1099 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 1100 while (ilink) { 1101 ierr = PetscLogEventBegin(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1102 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 1103 ierr = PetscLogEventEnd(ilink->event,ilink->ksp,ilink->x,ilink->y,NULL);CHKERRQ(ierr); 1104 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1105 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1106 pc->failedreason = PC_SUBPC_ERROR; 1107 } 1108 ilink = ilink->next; 1109 } 1110 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 1111 } else { 1112 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1113 while (ilink) { 1114 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1115 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1116 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1117 pc->failedreason = PC_SUBPC_ERROR; 1118 } 1119 ilink = ilink->next; 1120 } 1121 } 1122 } else { 1123 if (!jac->w1) { 1124 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 1125 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 1126 } 1127 ierr = VecSet(y,0.0);CHKERRQ(ierr); 1128 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1129 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1130 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1131 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1132 pc->failedreason = PC_SUBPC_ERROR; 1133 } 1134 while (ilink->next) { 1135 ilink = ilink->next; 1136 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1137 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1138 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1139 } 1140 while (ilink->previous) { 1141 ilink = ilink->previous; 1142 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1143 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1144 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1145 } 1146 } else { 1147 while (ilink->next) { /* get to last entry in linked list */ 1148 ilink = ilink->next; 1149 } 1150 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1151 ierr = KSPGetConvergedReason(ilink->ksp,&reason);CHKERRQ(ierr); 1152 if (reason == KSP_DIVERGED_PCSETUP_FAILED) { 1153 pc->failedreason = PC_SUBPC_ERROR; 1154 } 1155 while (ilink->previous) { 1156 ilink = ilink->previous; 1157 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1158 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1159 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1160 } 1161 } 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "PCReset_FieldSplit" 1168 static PetscErrorCode PCReset_FieldSplit(PC pc) 1169 { 1170 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1171 PetscErrorCode ierr; 1172 PC_FieldSplitLink ilink = jac->head,next; 1173 1174 PetscFunctionBegin; 1175 while (ilink) { 1176 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1177 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1178 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1179 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1180 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1181 if (!ilink->is_orig) { /* save the original IS */ 1182 ierr = PetscObjectReference((PetscObject)ilink->is);CHKERRQ(ierr); 1183 ilink->is_orig = ilink->is; 1184 } 1185 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1186 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1187 next = ilink->next; 1188 ilink = next; 1189 } 1190 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1191 if (jac->mat && jac->mat != jac->pmat) { 1192 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1193 } else if (jac->mat) { 1194 jac->mat = NULL; 1195 } 1196 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1197 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1198 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1199 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1200 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1201 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1202 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1203 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1204 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1205 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1206 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1207 jac->reset = PETSC_TRUE; 1208 jac->isrestrict = PETSC_FALSE; 1209 PetscFunctionReturn(0); 1210 } 1211 1212 #undef __FUNCT__ 1213 #define __FUNCT__ "PCDestroy_FieldSplit" 1214 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1215 { 1216 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1217 PetscErrorCode ierr; 1218 PC_FieldSplitLink ilink = jac->head,next; 1219 1220 PetscFunctionBegin; 1221 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1222 while (ilink) { 1223 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1224 ierr = ISDestroy(&ilink->is_orig);CHKERRQ(ierr); 1225 next = ilink->next; 1226 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1227 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1228 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1229 ierr = PetscFree(ilink);CHKERRQ(ierr); 1230 ilink = next; 1231 } 1232 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1233 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1234 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1235 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1236 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1237 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1238 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1239 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1240 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1241 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1242 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",NULL);CHKERRQ(ierr); 1243 PetscFunctionReturn(0); 1244 } 1245 1246 #undef __FUNCT__ 1247 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1248 static PetscErrorCode PCSetFromOptions_FieldSplit(PetscOptionItems *PetscOptionsObject,PC pc) 1249 { 1250 PetscErrorCode ierr; 1251 PetscInt bs; 1252 PetscBool flg,stokes = PETSC_FALSE; 1253 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1254 PCCompositeType ctype; 1255 1256 PetscFunctionBegin; 1257 ierr = PetscOptionsHead(PetscOptionsObject,"FieldSplit options");CHKERRQ(ierr); 1258 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1259 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1260 if (flg) { 1261 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1262 } 1263 jac->diag_use_amat = pc->useAmat; 1264 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); 1265 jac->offdiag_use_amat = pc->useAmat; 1266 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); 1267 /* FIXME: No programmatic equivalent to the following. */ 1268 ierr = PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1269 if (stokes) { 1270 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1271 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1272 } 1273 1274 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1275 if (flg) { 1276 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1277 } 1278 /* Only setup fields once */ 1279 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1280 /* only allow user to set fields from command line if bs is already known. 1281 otherwise user can set them in PCFieldSplitSetDefaults() */ 1282 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1283 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1284 } 1285 if (jac->type == PC_COMPOSITE_SCHUR) { 1286 ierr = PetscOptionsGetEnum(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1287 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1288 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); 1289 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1290 } 1291 ierr = PetscOptionsTail();CHKERRQ(ierr); 1292 PetscFunctionReturn(0); 1293 } 1294 1295 /*------------------------------------------------------------------------------------*/ 1296 1297 #undef __FUNCT__ 1298 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1299 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1300 { 1301 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1302 PetscErrorCode ierr; 1303 PC_FieldSplitLink ilink,next = jac->head; 1304 char prefix[128]; 1305 PetscInt i; 1306 1307 PetscFunctionBegin; 1308 if (jac->splitdefined) { 1309 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1310 PetscFunctionReturn(0); 1311 } 1312 for (i=0; i<n; i++) { 1313 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1314 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1315 } 1316 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1317 if (splitname) { 1318 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1319 } else { 1320 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1321 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1322 } 1323 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 */ 1324 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1325 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1326 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1327 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1328 1329 ilink->nfields = n; 1330 ilink->next = NULL; 1331 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1332 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1333 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1334 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1335 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1336 1337 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1338 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1339 1340 if (!next) { 1341 jac->head = ilink; 1342 ilink->previous = NULL; 1343 } else { 1344 while (next->next) { 1345 next = next->next; 1346 } 1347 next->next = ilink; 1348 ilink->previous = next; 1349 } 1350 jac->nsplits++; 1351 PetscFunctionReturn(0); 1352 } 1353 1354 #undef __FUNCT__ 1355 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1356 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1357 { 1358 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1359 PetscErrorCode ierr; 1360 1361 PetscFunctionBegin; 1362 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1363 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1364 1365 (*subksp)[1] = jac->kspschur; 1366 if (n) *n = jac->nsplits; 1367 PetscFunctionReturn(0); 1368 } 1369 1370 #undef __FUNCT__ 1371 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1372 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1373 { 1374 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1375 PetscErrorCode ierr; 1376 PetscInt cnt = 0; 1377 PC_FieldSplitLink ilink = jac->head; 1378 1379 PetscFunctionBegin; 1380 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1381 while (ilink) { 1382 (*subksp)[cnt++] = ilink->ksp; 1383 ilink = ilink->next; 1384 } 1385 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); 1386 if (n) *n = jac->nsplits; 1387 PetscFunctionReturn(0); 1388 } 1389 1390 #undef __FUNCT__ 1391 #define __FUNCT__ "PCFieldSplitRestrictIS" 1392 /*@C 1393 PCFieldSplitRestrictIS - Restricts the fieldsplit ISs to be within a given IS. 1394 1395 Input Parameters: 1396 + pc - the preconditioner context 1397 + is - the index set that defines the indices to which the fieldsplit is to be restricted 1398 1399 Level: advanced 1400 1401 @*/ 1402 PetscErrorCode PCFieldSplitRestrictIS(PC pc,IS isy) 1403 { 1404 PetscErrorCode ierr; 1405 1406 PetscFunctionBegin; 1407 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1408 PetscValidHeaderSpecific(isy,IS_CLASSID,2); 1409 ierr = PetscTryMethod(pc,"PCFieldSplitRestrictIS_C",(PC,IS),(pc,isy));CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 1413 1414 #undef __FUNCT__ 1415 #define __FUNCT__ "PCFieldSplitRestrictIS_FieldSplit" 1416 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1417 { 1418 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1419 PetscErrorCode ierr; 1420 PC_FieldSplitLink ilink = jac->head, next; 1421 PetscInt localsize,size,sizez,i; 1422 const PetscInt *ind, *indz; 1423 PetscInt *indc, *indcz; 1424 PetscBool flg; 1425 1426 PetscFunctionBegin; 1427 ierr = ISGetLocalSize(isy,&localsize);CHKERRQ(ierr); 1428 ierr = MPI_Scan(&localsize,&size,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)isy));CHKERRQ(ierr); 1429 size -= localsize; 1430 while(ilink) { 1431 IS isrl,isr; 1432 PC subpc; 1433 if (jac->reset) { 1434 ierr = ISEmbed(ilink->is_orig, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1435 } else { 1436 ierr = ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl);CHKERRQ(ierr); 1437 } 1438 ierr = ISGetLocalSize(isrl,&localsize);CHKERRQ(ierr); 1439 ierr = PetscMalloc1(localsize,&indc);CHKERRQ(ierr); 1440 ierr = ISGetIndices(isrl,&ind);CHKERRQ(ierr); 1441 ierr = PetscMemcpy(indc,ind,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1442 ierr = ISRestoreIndices(isrl,&ind);CHKERRQ(ierr); 1443 ierr = ISDestroy(&isrl);CHKERRQ(ierr); 1444 for (i=0; i<localsize; i++) *(indc+i) += size; 1445 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)isy),localsize,indc,PETSC_OWN_POINTER,&isr);CHKERRQ(ierr); 1446 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1447 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1448 ilink->is = isr; 1449 ierr = PetscObjectReference((PetscObject)isr);CHKERRQ(ierr); 1450 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1451 ilink->is_col = isr; 1452 ierr = ISDestroy(&isr);CHKERRQ(ierr); 1453 ierr = KSPGetPC(ilink->ksp, &subpc);CHKERRQ(ierr); 1454 ierr = PetscObjectTypeCompare((PetscObject)subpc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1455 if(flg) { 1456 IS iszl,isz; 1457 MPI_Comm comm; 1458 if (jac->reset) { 1459 ierr = ISGetLocalSize(ilink->is_orig,&localsize);CHKERRQ(ierr); 1460 comm = PetscObjectComm((PetscObject)ilink->is_orig); 1461 ierr = ISEmbed(isy, ilink->is_orig, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1462 } else { 1463 ierr = ISGetLocalSize(ilink->is,&localsize);CHKERRQ(ierr); 1464 comm = PetscObjectComm((PetscObject)ilink->is); 1465 ierr = ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl);CHKERRQ(ierr); 1466 } 1467 ierr = MPI_Scan(&localsize,&sizez,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); 1468 sizez -= localsize; 1469 ierr = ISGetLocalSize(iszl,&localsize);CHKERRQ(ierr); 1470 ierr = PetscMalloc1(localsize,&indcz);CHKERRQ(ierr); 1471 ierr = ISGetIndices(iszl,&indz);CHKERRQ(ierr); 1472 ierr = PetscMemcpy(indcz,indz,localsize*sizeof(PetscInt));CHKERRQ(ierr); 1473 ierr = ISRestoreIndices(iszl,&indz);CHKERRQ(ierr); 1474 ierr = ISDestroy(&iszl);CHKERRQ(ierr); 1475 for (i=0; i<localsize; i++) *(indcz+i) += sizez; 1476 ierr = ISCreateGeneral(comm,localsize,indcz,PETSC_OWN_POINTER,&isz);CHKERRQ(ierr); 1477 ierr = PCFieldSplitRestrictIS(subpc,isz);CHKERRQ(ierr); 1478 ierr = ISDestroy(&isz);CHKERRQ(ierr); 1479 } 1480 next = ilink->next; 1481 ilink = next; 1482 } 1483 jac->isrestrict = PETSC_TRUE; 1484 PetscFunctionReturn(0); 1485 } 1486 1487 #undef __FUNCT__ 1488 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1489 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1490 { 1491 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1492 PetscErrorCode ierr; 1493 PC_FieldSplitLink ilink, next = jac->head; 1494 char prefix[128]; 1495 1496 PetscFunctionBegin; 1497 if (jac->splitdefined) { 1498 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1499 PetscFunctionReturn(0); 1500 } 1501 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1502 if (splitname) { 1503 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1504 } else { 1505 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1506 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1507 } 1508 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 */ 1509 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1510 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1511 ilink->is = is; 1512 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1513 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1514 ilink->is_col = is; 1515 ilink->next = NULL; 1516 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1517 ierr = KSPSetErrorIfNotConverged(ilink->ksp,pc->erroriffailure);CHKERRQ(ierr); 1518 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1519 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1520 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1521 1522 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1523 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1524 1525 if (!next) { 1526 jac->head = ilink; 1527 ilink->previous = NULL; 1528 } else { 1529 while (next->next) { 1530 next = next->next; 1531 } 1532 next->next = ilink; 1533 ilink->previous = next; 1534 } 1535 jac->nsplits++; 1536 PetscFunctionReturn(0); 1537 } 1538 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "PCFieldSplitSetFields" 1541 /*@ 1542 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1543 1544 Logically Collective on PC 1545 1546 Input Parameters: 1547 + pc - the preconditioner context 1548 . splitname - name of this split, if NULL the number of the split is used 1549 . n - the number of fields in this split 1550 - fields - the fields in this split 1551 1552 Level: intermediate 1553 1554 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1555 1556 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1557 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 1558 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1559 where the numbered entries indicate what is in the field. 1560 1561 This function is called once per split (it creates a new split each time). Solve options 1562 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1563 1564 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1565 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1566 available when this routine is called. 1567 1568 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1569 1570 @*/ 1571 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1572 { 1573 PetscErrorCode ierr; 1574 1575 PetscFunctionBegin; 1576 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1577 PetscValidCharPointer(splitname,2); 1578 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1579 PetscValidIntPointer(fields,3); 1580 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1581 PetscFunctionReturn(0); 1582 } 1583 1584 #undef __FUNCT__ 1585 #define __FUNCT__ "PCFieldSplitSetDiagUseAmat" 1586 /*@ 1587 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1588 1589 Logically Collective on PC 1590 1591 Input Parameters: 1592 + pc - the preconditioner object 1593 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1594 1595 Options Database: 1596 . -pc_fieldsplit_diag_use_amat 1597 1598 Level: intermediate 1599 1600 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1601 1602 @*/ 1603 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1604 { 1605 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1606 PetscBool isfs; 1607 PetscErrorCode ierr; 1608 1609 PetscFunctionBegin; 1610 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1611 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1612 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1613 jac->diag_use_amat = flg; 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1619 /*@ 1620 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1621 1622 Logically Collective on PC 1623 1624 Input Parameters: 1625 . pc - the preconditioner object 1626 1627 Output Parameters: 1628 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1629 1630 1631 Level: intermediate 1632 1633 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1634 1635 @*/ 1636 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1637 { 1638 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1639 PetscBool isfs; 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1644 PetscValidPointer(flg,2); 1645 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1646 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1647 *flg = jac->diag_use_amat; 1648 PetscFunctionReturn(0); 1649 } 1650 1651 #undef __FUNCT__ 1652 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1653 /*@ 1654 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1655 1656 Logically Collective on PC 1657 1658 Input Parameters: 1659 + pc - the preconditioner object 1660 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1661 1662 Options Database: 1663 . -pc_fieldsplit_off_diag_use_amat 1664 1665 Level: intermediate 1666 1667 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1668 1669 @*/ 1670 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1671 { 1672 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1673 PetscBool isfs; 1674 PetscErrorCode ierr; 1675 1676 PetscFunctionBegin; 1677 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1678 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1679 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1680 jac->offdiag_use_amat = flg; 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1686 /*@ 1687 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1688 1689 Logically Collective on PC 1690 1691 Input Parameters: 1692 . pc - the preconditioner object 1693 1694 Output Parameters: 1695 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1696 1697 1698 Level: intermediate 1699 1700 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1701 1702 @*/ 1703 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1704 { 1705 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1706 PetscBool isfs; 1707 PetscErrorCode ierr; 1708 1709 PetscFunctionBegin; 1710 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1711 PetscValidPointer(flg,2); 1712 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1713 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1714 *flg = jac->offdiag_use_amat; 1715 PetscFunctionReturn(0); 1716 } 1717 1718 1719 1720 #undef __FUNCT__ 1721 #define __FUNCT__ "PCFieldSplitSetIS" 1722 /*@C 1723 PCFieldSplitSetIS - Sets the exact elements for field 1724 1725 Logically Collective on PC 1726 1727 Input Parameters: 1728 + pc - the preconditioner context 1729 . splitname - name of this split, if NULL the number of the split is used 1730 - is - the index set that defines the vector elements in this field 1731 1732 1733 Notes: 1734 Use PCFieldSplitSetFields(), for fields defined by strided types. 1735 1736 This function is called once per split (it creates a new split each time). Solve options 1737 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1738 1739 Level: intermediate 1740 1741 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1742 1743 @*/ 1744 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1745 { 1746 PetscErrorCode ierr; 1747 1748 PetscFunctionBegin; 1749 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1750 if (splitname) PetscValidCharPointer(splitname,2); 1751 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1752 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1753 PetscFunctionReturn(0); 1754 } 1755 1756 #undef __FUNCT__ 1757 #define __FUNCT__ "PCFieldSplitGetIS" 1758 /*@ 1759 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1760 1761 Logically Collective on PC 1762 1763 Input Parameters: 1764 + pc - the preconditioner context 1765 - splitname - name of this split 1766 1767 Output Parameter: 1768 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1769 1770 Level: intermediate 1771 1772 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1773 1774 @*/ 1775 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1776 { 1777 PetscErrorCode ierr; 1778 1779 PetscFunctionBegin; 1780 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1781 PetscValidCharPointer(splitname,2); 1782 PetscValidPointer(is,3); 1783 { 1784 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1785 PC_FieldSplitLink ilink = jac->head; 1786 PetscBool found; 1787 1788 *is = NULL; 1789 while (ilink) { 1790 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1791 if (found) { 1792 *is = ilink->is; 1793 break; 1794 } 1795 ilink = ilink->next; 1796 } 1797 } 1798 PetscFunctionReturn(0); 1799 } 1800 1801 #undef __FUNCT__ 1802 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1803 /*@ 1804 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1805 fieldsplit preconditioner. If not set the matrix block size is used. 1806 1807 Logically Collective on PC 1808 1809 Input Parameters: 1810 + pc - the preconditioner context 1811 - bs - the block size 1812 1813 Level: intermediate 1814 1815 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1816 1817 @*/ 1818 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1819 { 1820 PetscErrorCode ierr; 1821 1822 PetscFunctionBegin; 1823 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1824 PetscValidLogicalCollectiveInt(pc,bs,2); 1825 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1831 /*@C 1832 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1833 1834 Collective on KSP 1835 1836 Input Parameter: 1837 . pc - the preconditioner context 1838 1839 Output Parameters: 1840 + n - the number of splits 1841 - pc - the array of KSP contexts 1842 1843 Note: 1844 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1845 (not the KSP just the array that contains them). 1846 1847 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1848 1849 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1850 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1851 KSP array must be. 1852 1853 1854 Level: advanced 1855 1856 .seealso: PCFIELDSPLIT 1857 @*/ 1858 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1859 { 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1864 if (n) PetscValidIntPointer(n,2); 1865 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 #undef __FUNCT__ 1870 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1871 /*@ 1872 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1873 A11 matrix. Otherwise no preconditioner is used. 1874 1875 Collective on PC 1876 1877 Input Parameters: 1878 + pc - the preconditioner context 1879 . 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 1880 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 1881 - userpre - matrix to use for preconditioning, or NULL 1882 1883 Options Database: 1884 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 1885 1886 Notes: 1887 $ If ptype is 1888 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1889 $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix 1890 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1891 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1892 $ preconditioner 1893 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1894 $ to this function). 1895 $ selfp then the preconditioning for the Schur complement is generated from an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1896 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 1897 $ lumped before extracting the diagonal using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 1898 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive) 1899 $ useful mostly as a test that the Schur complement approach can work for your problem 1900 1901 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1902 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1903 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1904 1905 Level: intermediate 1906 1907 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1908 MatSchurComplementSetAinvType(), PCLSC 1909 1910 @*/ 1911 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1912 { 1913 PetscErrorCode ierr; 1914 1915 PetscFunctionBegin; 1916 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1917 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1921 1922 #undef __FUNCT__ 1923 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1924 /*@ 1925 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1926 preconditioned. See PCFieldSplitSetSchurPre() for details. 1927 1928 Logically Collective on PC 1929 1930 Input Parameters: 1931 . pc - the preconditioner context 1932 1933 Output Parameters: 1934 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1935 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1936 1937 Level: intermediate 1938 1939 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1940 1941 @*/ 1942 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1943 { 1944 PetscErrorCode ierr; 1945 1946 PetscFunctionBegin; 1947 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1948 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 #undef __FUNCT__ 1953 #define __FUNCT__ "PCFieldSplitSchurGetS" 1954 /*@ 1955 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1956 1957 Not collective 1958 1959 Input Parameter: 1960 . pc - the preconditioner context 1961 1962 Output Parameter: 1963 . S - the Schur complement matrix 1964 1965 Notes: 1966 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1967 1968 Level: advanced 1969 1970 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1971 1972 @*/ 1973 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1974 { 1975 PetscErrorCode ierr; 1976 const char* t; 1977 PetscBool isfs; 1978 PC_FieldSplit *jac; 1979 1980 PetscFunctionBegin; 1981 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1982 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1983 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1984 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1985 jac = (PC_FieldSplit*)pc->data; 1986 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1987 if (S) *S = jac->schur; 1988 PetscFunctionReturn(0); 1989 } 1990 1991 #undef __FUNCT__ 1992 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1993 /*@ 1994 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1995 1996 Not collective 1997 1998 Input Parameters: 1999 + pc - the preconditioner context 2000 . S - the Schur complement matrix 2001 2002 Level: advanced 2003 2004 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 2005 2006 @*/ 2007 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2008 { 2009 PetscErrorCode ierr; 2010 const char* t; 2011 PetscBool isfs; 2012 PC_FieldSplit *jac; 2013 2014 PetscFunctionBegin; 2015 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2016 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 2017 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2018 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2019 jac = (PC_FieldSplit*)pc->data; 2020 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2021 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2022 PetscFunctionReturn(0); 2023 } 2024 2025 2026 #undef __FUNCT__ 2027 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 2028 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2029 { 2030 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2031 PetscErrorCode ierr; 2032 2033 PetscFunctionBegin; 2034 jac->schurpre = ptype; 2035 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2036 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 2037 jac->schur_user = pre; 2038 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 2039 } 2040 PetscFunctionReturn(0); 2041 } 2042 2043 #undef __FUNCT__ 2044 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 2045 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2046 { 2047 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2048 2049 PetscFunctionBegin; 2050 *ptype = jac->schurpre; 2051 *pre = jac->schur_user; 2052 PetscFunctionReturn(0); 2053 } 2054 2055 #undef __FUNCT__ 2056 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 2057 /*@ 2058 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 2059 2060 Collective on PC 2061 2062 Input Parameters: 2063 + pc - the preconditioner context 2064 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2065 2066 Options Database: 2067 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 2068 2069 2070 Level: intermediate 2071 2072 Notes: 2073 The FULL factorization is 2074 2075 $ (A B) = (1 0) (A 0) (1 Ainv*B) 2076 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 2077 2078 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, 2079 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). 2080 2081 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 2082 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 2083 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, 2084 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2085 2086 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 2087 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). 2088 2089 References: 2090 + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2091 - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2092 2093 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 2094 @*/ 2095 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2096 { 2097 PetscErrorCode ierr; 2098 2099 PetscFunctionBegin; 2100 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2101 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 2102 PetscFunctionReturn(0); 2103 } 2104 2105 #undef __FUNCT__ 2106 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 2107 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2108 { 2109 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2110 2111 PetscFunctionBegin; 2112 jac->schurfactorization = ftype; 2113 PetscFunctionReturn(0); 2114 } 2115 2116 #undef __FUNCT__ 2117 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 2118 /*@C 2119 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2120 2121 Collective on KSP 2122 2123 Input Parameter: 2124 . pc - the preconditioner context 2125 2126 Output Parameters: 2127 + A00 - the (0,0) block 2128 . A01 - the (0,1) block 2129 . A10 - the (1,0) block 2130 - A11 - the (1,1) block 2131 2132 Level: advanced 2133 2134 .seealso: PCFIELDSPLIT 2135 @*/ 2136 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2137 { 2138 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2139 2140 PetscFunctionBegin; 2141 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2142 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2143 if (A00) *A00 = jac->pmat[0]; 2144 if (A01) *A01 = jac->B; 2145 if (A10) *A10 = jac->C; 2146 if (A11) *A11 = jac->pmat[1]; 2147 PetscFunctionReturn(0); 2148 } 2149 2150 #undef __FUNCT__ 2151 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 2152 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2153 { 2154 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2155 PetscErrorCode ierr; 2156 2157 PetscFunctionBegin; 2158 jac->type = type; 2159 if (type == PC_COMPOSITE_SCHUR) { 2160 pc->ops->apply = PCApply_FieldSplit_Schur; 2161 pc->ops->view = PCView_FieldSplit_Schur; 2162 2163 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2164 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2165 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2166 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2167 2168 } else { 2169 pc->ops->apply = PCApply_FieldSplit; 2170 pc->ops->view = PCView_FieldSplit; 2171 2172 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2173 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2174 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2175 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2176 } 2177 PetscFunctionReturn(0); 2178 } 2179 2180 #undef __FUNCT__ 2181 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 2182 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2183 { 2184 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2185 2186 PetscFunctionBegin; 2187 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2188 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); 2189 jac->bs = bs; 2190 PetscFunctionReturn(0); 2191 } 2192 2193 #undef __FUNCT__ 2194 #define __FUNCT__ "PCFieldSplitSetType" 2195 /*@ 2196 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2197 2198 Collective on PC 2199 2200 Input Parameter: 2201 . pc - the preconditioner context 2202 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2203 2204 Options Database Key: 2205 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2206 2207 Level: Intermediate 2208 2209 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2210 2211 .seealso: PCCompositeSetType() 2212 2213 @*/ 2214 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2215 { 2216 PetscErrorCode ierr; 2217 2218 PetscFunctionBegin; 2219 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2220 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2221 PetscFunctionReturn(0); 2222 } 2223 2224 #undef __FUNCT__ 2225 #define __FUNCT__ "PCFieldSplitGetType" 2226 /*@ 2227 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2228 2229 Not collective 2230 2231 Input Parameter: 2232 . pc - the preconditioner context 2233 2234 Output Parameter: 2235 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2236 2237 Level: Intermediate 2238 2239 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2240 .seealso: PCCompositeSetType() 2241 @*/ 2242 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2243 { 2244 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2245 2246 PetscFunctionBegin; 2247 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2248 PetscValidIntPointer(type,2); 2249 *type = jac->type; 2250 PetscFunctionReturn(0); 2251 } 2252 2253 #undef __FUNCT__ 2254 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2255 /*@ 2256 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2257 2258 Logically Collective 2259 2260 Input Parameters: 2261 + pc - the preconditioner context 2262 - flg - boolean indicating whether to use field splits defined by the DM 2263 2264 Options Database Key: 2265 . -pc_fieldsplit_dm_splits 2266 2267 Level: Intermediate 2268 2269 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2270 2271 .seealso: PCFieldSplitGetDMSplits() 2272 2273 @*/ 2274 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2275 { 2276 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2277 PetscBool isfs; 2278 PetscErrorCode ierr; 2279 2280 PetscFunctionBegin; 2281 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2282 PetscValidLogicalCollectiveBool(pc,flg,2); 2283 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2284 if (isfs) { 2285 jac->dm_splits = flg; 2286 } 2287 PetscFunctionReturn(0); 2288 } 2289 2290 2291 #undef __FUNCT__ 2292 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2293 /*@ 2294 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2295 2296 Logically Collective 2297 2298 Input Parameter: 2299 . pc - the preconditioner context 2300 2301 Output Parameter: 2302 . flg - boolean indicating whether to use field splits defined by the DM 2303 2304 Level: Intermediate 2305 2306 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2307 2308 .seealso: PCFieldSplitSetDMSplits() 2309 2310 @*/ 2311 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2312 { 2313 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2314 PetscBool isfs; 2315 PetscErrorCode ierr; 2316 2317 PetscFunctionBegin; 2318 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2319 PetscValidPointer(flg,2); 2320 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2321 if (isfs) { 2322 if(flg) *flg = jac->dm_splits; 2323 } 2324 PetscFunctionReturn(0); 2325 } 2326 2327 /* -------------------------------------------------------------------------------------*/ 2328 /*MC 2329 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2330 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2331 2332 To set options on the solvers for each block append -fieldsplit_ to all the PC 2333 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2334 2335 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2336 and set the options directly on the resulting KSP object 2337 2338 Level: intermediate 2339 2340 Options Database Keys: 2341 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2342 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2343 been supplied explicitly by -pc_fieldsplit_%d_fields 2344 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2345 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2346 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2347 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2348 2349 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2350 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2351 2352 Notes: 2353 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2354 to define a field by an arbitrary collection of entries. 2355 2356 If no fields are set the default is used. The fields are defined by entries strided by bs, 2357 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2358 if this is not called the block size defaults to the blocksize of the second matrix passed 2359 to KSPSetOperators()/PCSetOperators(). 2360 2361 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2362 $ ( A10 A11 ) 2363 $ the preconditioner using full factorization is 2364 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2365 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2366 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2367 $ S = A11 - A10 ksp(A00) A01 2368 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 2369 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2370 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2371 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2372 2373 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2374 diag gives 2375 $ ( inv(A00) 0 ) 2376 $ ( 0 -ksp(S) ) 2377 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 2378 $ ( A00 0 ) 2379 $ ( A10 S ) 2380 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2381 $ ( A00 A01 ) 2382 $ ( 0 S ) 2383 where again the inverses of A00 and S are applied using KSPs. 2384 2385 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2386 is used automatically for a second block. 2387 2388 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2389 Generally it should be used with the AIJ format. 2390 2391 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2392 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2393 inside a smoother resulting in "Distributive Smoothers". 2394 2395 Concepts: physics based preconditioners, block preconditioners 2396 2397 There is a nice discussion of block preconditioners in 2398 2399 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2400 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2401 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2402 2403 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 2404 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 2405 of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables. 2406 2407 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2408 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2409 MatSchurComplementSetAinvType() 2410 M*/ 2411 2412 #undef __FUNCT__ 2413 #define __FUNCT__ "PCCreate_FieldSplit" 2414 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2415 { 2416 PetscErrorCode ierr; 2417 PC_FieldSplit *jac; 2418 2419 PetscFunctionBegin; 2420 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2421 2422 jac->bs = -1; 2423 jac->nsplits = 0; 2424 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2425 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2426 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2427 jac->dm_splits = PETSC_TRUE; 2428 2429 pc->data = (void*)jac; 2430 2431 pc->ops->apply = PCApply_FieldSplit; 2432 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2433 pc->ops->setup = PCSetUp_FieldSplit; 2434 pc->ops->reset = PCReset_FieldSplit; 2435 pc->ops->destroy = PCDestroy_FieldSplit; 2436 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2437 pc->ops->view = PCView_FieldSplit; 2438 pc->ops->applyrichardson = 0; 2439 2440 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2441 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2442 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2443 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2444 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2445 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr); 2446 PetscFunctionReturn(0); 2447 } 2448 2449 2450 2451