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