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