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