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]);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,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1090 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1091 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1092 PetscFunctionReturn(0); 1093 } 1094 1095 #undef __FUNCT__ 1096 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1097 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1098 { 1099 PetscErrorCode ierr; 1100 PetscInt bs; 1101 PetscBool flg,stokes = PETSC_FALSE; 1102 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1103 PCCompositeType ctype; 1104 1105 PetscFunctionBegin; 1106 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1107 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1108 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1109 if (flg) { 1110 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1111 } 1112 1113 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1114 if (stokes) { 1115 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1116 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1117 } 1118 1119 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1120 if (flg) { 1121 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1122 } 1123 /* Only setup fields once */ 1124 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1125 /* only allow user to set fields from command line if bs is already known. 1126 otherwise user can set them in PCFieldSplitSetDefaults() */ 1127 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1128 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1129 } 1130 if (jac->type == PC_COMPOSITE_SCHUR) { 1131 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1132 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1133 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); 1134 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1135 } 1136 ierr = PetscOptionsTail();CHKERRQ(ierr); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /*------------------------------------------------------------------------------------*/ 1141 1142 #undef __FUNCT__ 1143 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1144 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1145 { 1146 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1147 PetscErrorCode ierr; 1148 PC_FieldSplitLink ilink,next = jac->head; 1149 char prefix[128]; 1150 PetscInt i; 1151 1152 PetscFunctionBegin; 1153 if (jac->splitdefined) { 1154 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 for (i=0; i<n; i++) { 1158 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1159 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1160 } 1161 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1162 if (splitname) { 1163 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1164 } else { 1165 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1166 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1167 } 1168 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1169 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1170 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1171 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1172 1173 ilink->nfields = n; 1174 ilink->next = NULL; 1175 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1176 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1177 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1178 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1179 1180 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1181 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1182 1183 if (!next) { 1184 jac->head = ilink; 1185 ilink->previous = NULL; 1186 } else { 1187 while (next->next) { 1188 next = next->next; 1189 } 1190 next->next = ilink; 1191 ilink->previous = next; 1192 } 1193 jac->nsplits++; 1194 PetscFunctionReturn(0); 1195 } 1196 1197 #undef __FUNCT__ 1198 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1199 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1200 { 1201 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1202 PetscErrorCode ierr; 1203 1204 PetscFunctionBegin; 1205 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1206 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1207 1208 (*subksp)[1] = jac->kspschur; 1209 if (n) *n = jac->nsplits; 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1215 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1216 { 1217 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1218 PetscErrorCode ierr; 1219 PetscInt cnt = 0; 1220 PC_FieldSplitLink ilink = jac->head; 1221 1222 PetscFunctionBegin; 1223 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1224 while (ilink) { 1225 (*subksp)[cnt++] = ilink->ksp; 1226 ilink = ilink->next; 1227 } 1228 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); 1229 if (n) *n = jac->nsplits; 1230 PetscFunctionReturn(0); 1231 } 1232 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1235 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1236 { 1237 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1238 PetscErrorCode ierr; 1239 PC_FieldSplitLink ilink, next = jac->head; 1240 char prefix[128]; 1241 1242 PetscFunctionBegin; 1243 if (jac->splitdefined) { 1244 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1245 PetscFunctionReturn(0); 1246 } 1247 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1248 if (splitname) { 1249 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1250 } else { 1251 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1252 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1253 } 1254 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1255 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1256 ilink->is = is; 1257 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1258 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1259 ilink->is_col = is; 1260 ilink->next = NULL; 1261 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1262 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1263 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1264 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1265 1266 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1267 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1268 1269 if (!next) { 1270 jac->head = ilink; 1271 ilink->previous = NULL; 1272 } else { 1273 while (next->next) { 1274 next = next->next; 1275 } 1276 next->next = ilink; 1277 ilink->previous = next; 1278 } 1279 jac->nsplits++; 1280 PetscFunctionReturn(0); 1281 } 1282 1283 #undef __FUNCT__ 1284 #define __FUNCT__ "PCFieldSplitSetFields" 1285 /*@ 1286 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1287 1288 Logically Collective on PC 1289 1290 Input Parameters: 1291 + pc - the preconditioner context 1292 . splitname - name of this split, if NULL the number of the split is used 1293 . n - the number of fields in this split 1294 - fields - the fields in this split 1295 1296 Level: intermediate 1297 1298 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1299 1300 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1301 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 1302 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1303 where the numbered entries indicate what is in the field. 1304 1305 This function is called once per split (it creates a new split each time). Solve options 1306 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1307 1308 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1309 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1310 available when this routine is called. 1311 1312 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1313 1314 @*/ 1315 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1316 { 1317 PetscErrorCode ierr; 1318 1319 PetscFunctionBegin; 1320 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1321 PetscValidCharPointer(splitname,2); 1322 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1323 PetscValidIntPointer(fields,3); 1324 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 #undef __FUNCT__ 1329 #define __FUNCT__ "PCFieldSplitSetIS" 1330 /*@ 1331 PCFieldSplitSetIS - Sets the exact elements for field 1332 1333 Logically Collective on PC 1334 1335 Input Parameters: 1336 + pc - the preconditioner context 1337 . splitname - name of this split, if NULL the number of the split is used 1338 - is - the index set that defines the vector elements in this field 1339 1340 1341 Notes: 1342 Use PCFieldSplitSetFields(), for fields defined by strided types. 1343 1344 This function is called once per split (it creates a new split each time). Solve options 1345 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1346 1347 Level: intermediate 1348 1349 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1350 1351 @*/ 1352 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1353 { 1354 PetscErrorCode ierr; 1355 1356 PetscFunctionBegin; 1357 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1358 if (splitname) PetscValidCharPointer(splitname,2); 1359 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1360 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1361 PetscFunctionReturn(0); 1362 } 1363 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "PCFieldSplitGetIS" 1366 /*@ 1367 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1368 1369 Logically Collective on PC 1370 1371 Input Parameters: 1372 + pc - the preconditioner context 1373 - splitname - name of this split 1374 1375 Output Parameter: 1376 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1377 1378 Level: intermediate 1379 1380 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1381 1382 @*/ 1383 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1384 { 1385 PetscErrorCode ierr; 1386 1387 PetscFunctionBegin; 1388 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1389 PetscValidCharPointer(splitname,2); 1390 PetscValidPointer(is,3); 1391 { 1392 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1393 PC_FieldSplitLink ilink = jac->head; 1394 PetscBool found; 1395 1396 *is = NULL; 1397 while (ilink) { 1398 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1399 if (found) { 1400 *is = ilink->is; 1401 break; 1402 } 1403 ilink = ilink->next; 1404 } 1405 } 1406 PetscFunctionReturn(0); 1407 } 1408 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1411 /*@ 1412 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1413 fieldsplit preconditioner. If not set the matrix block size is used. 1414 1415 Logically Collective on PC 1416 1417 Input Parameters: 1418 + pc - the preconditioner context 1419 - bs - the block size 1420 1421 Level: intermediate 1422 1423 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1424 1425 @*/ 1426 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1427 { 1428 PetscErrorCode ierr; 1429 1430 PetscFunctionBegin; 1431 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1432 PetscValidLogicalCollectiveInt(pc,bs,2); 1433 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1434 PetscFunctionReturn(0); 1435 } 1436 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1439 /*@C 1440 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1441 1442 Collective on KSP 1443 1444 Input Parameter: 1445 . pc - the preconditioner context 1446 1447 Output Parameters: 1448 + n - the number of splits 1449 - pc - the array of KSP contexts 1450 1451 Note: 1452 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1453 (not the KSP just the array that contains them). 1454 1455 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1456 1457 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1458 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1459 KSP array must be. 1460 1461 1462 Level: advanced 1463 1464 .seealso: PCFIELDSPLIT 1465 @*/ 1466 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1467 { 1468 PetscErrorCode ierr; 1469 1470 PetscFunctionBegin; 1471 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1472 if (n) PetscValidIntPointer(n,2); 1473 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1474 PetscFunctionReturn(0); 1475 } 1476 1477 #undef __FUNCT__ 1478 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1479 /*@ 1480 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1481 A11 matrix. Otherwise no preconditioner is used. 1482 1483 Collective on PC 1484 1485 Input Parameters: 1486 + pc - the preconditioner context 1487 . 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 1488 - userpre - matrix to use for preconditioning, or NULL 1489 1490 Options Database: 1491 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11 1492 1493 Notes: 1494 $ If ptype is 1495 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1496 $ to this function). 1497 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1498 $ matrix associated with the Schur complement (i.e. A11) 1499 $ full then the preconditioner uses the exact Schur complement (this is expensive) 1500 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1501 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1502 $ preconditioner 1503 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1504 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. 1505 1506 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1507 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1508 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1509 1510 Level: intermediate 1511 1512 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1513 1514 @*/ 1515 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1516 { 1517 PetscErrorCode ierr; 1518 1519 PetscFunctionBegin; 1520 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1521 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1522 PetscFunctionReturn(0); 1523 } 1524 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1528 /*@ 1529 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1530 preconditioned. See PCFieldSplitSetSchurPre() for details. 1531 1532 Logically Collective on PC 1533 1534 Input Parameters: 1535 . pc - the preconditioner context 1536 1537 Output Parameters: 1538 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1539 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1540 1541 Level: intermediate 1542 1543 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1544 1545 @*/ 1546 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1547 { 1548 PetscErrorCode ierr; 1549 1550 PetscFunctionBegin; 1551 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1552 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1553 PetscFunctionReturn(0); 1554 } 1555 1556 #undef __FUNCT__ 1557 #define __FUNCT__ "PCFieldSplitSchurGetS" 1558 /*@ 1559 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1560 1561 Not collective 1562 1563 Input Parameter: 1564 . pc - the preconditioner context 1565 1566 Output Parameter: 1567 . S - the Schur complement matrix 1568 1569 Notes: 1570 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1571 1572 Level: advanced 1573 1574 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1575 1576 @*/ 1577 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1578 { 1579 PetscErrorCode ierr; 1580 const char* t; 1581 PetscBool isfs; 1582 PC_FieldSplit *jac; 1583 1584 PetscFunctionBegin; 1585 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1586 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1587 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1588 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1589 jac = (PC_FieldSplit*)pc->data; 1590 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1591 if (S) *S = jac->schur; 1592 PetscFunctionReturn(0); 1593 } 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1597 /*@ 1598 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1599 1600 Not collective 1601 1602 Input Parameters: 1603 + pc - the preconditioner context 1604 . S - the Schur complement matrix 1605 1606 Level: advanced 1607 1608 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1609 1610 @*/ 1611 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1612 { 1613 PetscErrorCode ierr; 1614 const char* t; 1615 PetscBool isfs; 1616 PC_FieldSplit *jac; 1617 1618 PetscFunctionBegin; 1619 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1620 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1621 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1622 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1623 jac = (PC_FieldSplit*)pc->data; 1624 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1625 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1626 PetscFunctionReturn(0); 1627 } 1628 1629 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 1632 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1633 { 1634 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1635 PetscErrorCode ierr; 1636 1637 PetscFunctionBegin; 1638 jac->schurpre = ptype; 1639 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1640 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1641 jac->schur_user = pre; 1642 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1643 } 1644 PetscFunctionReturn(0); 1645 } 1646 1647 #undef __FUNCT__ 1648 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 1649 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1650 { 1651 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1652 1653 PetscFunctionBegin; 1654 *ptype = jac->schurpre; 1655 *pre = jac->schur_user; 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1661 /*@ 1662 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1663 1664 Collective on PC 1665 1666 Input Parameters: 1667 + pc - the preconditioner context 1668 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1669 1670 Options Database: 1671 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1672 1673 1674 Level: intermediate 1675 1676 Notes: 1677 The FULL factorization is 1678 1679 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1680 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1681 1682 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, 1683 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). 1684 1685 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 1686 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 1687 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, 1688 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1689 1690 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 1691 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). 1692 1693 References: 1694 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1695 1696 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1697 1698 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1699 @*/ 1700 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1701 { 1702 PetscErrorCode ierr; 1703 1704 PetscFunctionBegin; 1705 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1706 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1712 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1713 { 1714 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1715 1716 PetscFunctionBegin; 1717 jac->schurfactorization = ftype; 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1723 /*@C 1724 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1725 1726 Collective on KSP 1727 1728 Input Parameter: 1729 . pc - the preconditioner context 1730 1731 Output Parameters: 1732 + A00 - the (0,0) block 1733 . A01 - the (0,1) block 1734 . A10 - the (1,0) block 1735 - A11 - the (1,1) block 1736 1737 Level: advanced 1738 1739 .seealso: PCFIELDSPLIT 1740 @*/ 1741 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1742 { 1743 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1744 1745 PetscFunctionBegin; 1746 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1747 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1748 if (A00) *A00 = jac->pmat[0]; 1749 if (A01) *A01 = jac->B; 1750 if (A10) *A10 = jac->C; 1751 if (A11) *A11 = jac->pmat[1]; 1752 PetscFunctionReturn(0); 1753 } 1754 1755 #undef __FUNCT__ 1756 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1757 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1758 { 1759 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1760 PetscErrorCode ierr; 1761 1762 PetscFunctionBegin; 1763 jac->type = type; 1764 if (type == PC_COMPOSITE_SCHUR) { 1765 pc->ops->apply = PCApply_FieldSplit_Schur; 1766 pc->ops->view = PCView_FieldSplit_Schur; 1767 1768 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1769 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 1770 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 1771 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1772 1773 } else { 1774 pc->ops->apply = PCApply_FieldSplit; 1775 pc->ops->view = PCView_FieldSplit; 1776 1777 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1778 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 1779 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1781 } 1782 PetscFunctionReturn(0); 1783 } 1784 1785 #undef __FUNCT__ 1786 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1787 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1788 { 1789 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1790 1791 PetscFunctionBegin; 1792 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1793 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); 1794 jac->bs = bs; 1795 PetscFunctionReturn(0); 1796 } 1797 1798 #undef __FUNCT__ 1799 #define __FUNCT__ "PCFieldSplitSetType" 1800 /*@ 1801 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1802 1803 Collective on PC 1804 1805 Input Parameter: 1806 . pc - the preconditioner context 1807 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1808 1809 Options Database Key: 1810 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1811 1812 Level: Intermediate 1813 1814 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1815 1816 .seealso: PCCompositeSetType() 1817 1818 @*/ 1819 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1820 { 1821 PetscErrorCode ierr; 1822 1823 PetscFunctionBegin; 1824 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1825 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "PCFieldSplitGetType" 1831 /*@ 1832 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1833 1834 Not collective 1835 1836 Input Parameter: 1837 . pc - the preconditioner context 1838 1839 Output Parameter: 1840 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1841 1842 Level: Intermediate 1843 1844 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1845 .seealso: PCCompositeSetType() 1846 @*/ 1847 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1848 { 1849 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1850 1851 PetscFunctionBegin; 1852 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1853 PetscValidIntPointer(type,2); 1854 *type = jac->type; 1855 PetscFunctionReturn(0); 1856 } 1857 1858 #undef __FUNCT__ 1859 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1860 /*@ 1861 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1862 1863 Logically Collective 1864 1865 Input Parameters: 1866 + pc - the preconditioner context 1867 - flg - boolean indicating whether to use field splits defined by the DM 1868 1869 Options Database Key: 1870 . -pc_fieldsplit_dm_splits 1871 1872 Level: Intermediate 1873 1874 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1875 1876 .seealso: PCFieldSplitGetDMSplits() 1877 1878 @*/ 1879 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1880 { 1881 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1882 PetscBool isfs; 1883 PetscErrorCode ierr; 1884 1885 PetscFunctionBegin; 1886 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1887 PetscValidLogicalCollectiveBool(pc,flg,2); 1888 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1889 if (isfs) { 1890 jac->dm_splits = flg; 1891 } 1892 PetscFunctionReturn(0); 1893 } 1894 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1898 /*@ 1899 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1900 1901 Logically Collective 1902 1903 Input Parameter: 1904 . pc - the preconditioner context 1905 1906 Output Parameter: 1907 . flg - boolean indicating whether to use field splits defined by the DM 1908 1909 Level: Intermediate 1910 1911 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1912 1913 .seealso: PCFieldSplitSetDMSplits() 1914 1915 @*/ 1916 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 1917 { 1918 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1919 PetscBool isfs; 1920 PetscErrorCode ierr; 1921 1922 PetscFunctionBegin; 1923 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1924 PetscValidPointer(flg,2); 1925 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1926 if (isfs) { 1927 if(flg) *flg = jac->dm_splits; 1928 } 1929 PetscFunctionReturn(0); 1930 } 1931 1932 /* -------------------------------------------------------------------------------------*/ 1933 /*MC 1934 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1935 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1936 1937 To set options on the solvers for each block append -fieldsplit_ to all the PC 1938 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1939 1940 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1941 and set the options directly on the resulting KSP object 1942 1943 Level: intermediate 1944 1945 Options Database Keys: 1946 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1947 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1948 been supplied explicitly by -pc_fieldsplit_%d_fields 1949 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1950 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1951 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11 1952 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1953 1954 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1955 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1956 1957 Notes: 1958 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1959 to define a field by an arbitrary collection of entries. 1960 1961 If no fields are set the default is used. The fields are defined by entries strided by bs, 1962 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1963 if this is not called the block size defaults to the blocksize of the second matrix passed 1964 to KSPSetOperators()/PCSetOperators(). 1965 1966 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1967 $ ( A10 A11 ) 1968 $ the preconditioner using full factorization is 1969 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 1970 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1971 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 1972 $ S = A11 - A10 ksp(A00) A01 1973 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 1974 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1975 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1976 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this 1977 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 1978 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 1979 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 1980 (e.g., -fieldsplit_1_pc_type asm). 1981 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1982 diag gives 1983 $ ( inv(A00) 0 ) 1984 $ ( 0 -ksp(S) ) 1985 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 1986 $ ( A00 0 ) 1987 $ ( A10 S ) 1988 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1989 $ ( A00 A01 ) 1990 $ ( 0 S ) 1991 where again the inverses of A00 and S are applied using KSPs. 1992 1993 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1994 is used automatically for a second block. 1995 1996 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1997 Generally it should be used with the AIJ format. 1998 1999 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2000 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2001 inside a smoother resulting in "Distributive Smoothers". 2002 2003 Concepts: physics based preconditioners, block preconditioners 2004 2005 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2006 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre() 2007 M*/ 2008 2009 #undef __FUNCT__ 2010 #define __FUNCT__ "PCCreate_FieldSplit" 2011 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2012 { 2013 PetscErrorCode ierr; 2014 PC_FieldSplit *jac; 2015 2016 PetscFunctionBegin; 2017 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2018 2019 jac->bs = -1; 2020 jac->nsplits = 0; 2021 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2022 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2023 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2024 jac->dm_splits = PETSC_TRUE; 2025 2026 pc->data = (void*)jac; 2027 2028 pc->ops->apply = PCApply_FieldSplit; 2029 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2030 pc->ops->setup = PCSetUp_FieldSplit; 2031 pc->ops->reset = PCReset_FieldSplit; 2032 pc->ops->destroy = PCDestroy_FieldSplit; 2033 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2034 pc->ops->view = PCView_FieldSplit; 2035 pc->ops->applyrichardson = 0; 2036 2037 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2038 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2039 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2040 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2041 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2042 PetscFunctionReturn(0); 2043 } 2044 2045 2046 2047