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