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