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