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