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