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