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 MatStructure flag = pc->flag; 482 PetscBool sorted, sorted_col; 483 484 PetscFunctionBegin; 485 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 486 nsplit = jac->nsplits; 487 ilink = jac->head; 488 489 /* get the matrices for each split */ 490 if (!jac->issetup) { 491 PetscInt rstart,rend,nslots,bs; 492 493 jac->issetup = PETSC_TRUE; 494 495 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 496 if (jac->defaultsplit || !ilink->is) { 497 if (jac->bs <= 0) jac->bs = nsplit; 498 } 499 bs = jac->bs; 500 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 501 nslots = (rend - rstart)/bs; 502 for (i=0; i<nsplit; i++) { 503 if (jac->defaultsplit) { 504 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 505 ierr = ISDuplicate(ilink->is,&ilink->is_col);CHKERRQ(ierr); 506 } else if (!ilink->is) { 507 if (ilink->nfields > 1) { 508 PetscInt *ii,*jj,j,k,nfields = ilink->nfields,*fields = ilink->fields,*fields_col = ilink->fields_col; 509 ierr = PetscMalloc1(ilink->nfields*nslots,&ii);CHKERRQ(ierr); 510 ierr = PetscMalloc1(ilink->nfields*nslots,&jj);CHKERRQ(ierr); 511 for (j=0; j<nslots; j++) { 512 for (k=0; k<nfields; k++) { 513 ii[nfields*j + k] = rstart + bs*j + fields[k]; 514 jj[nfields*j + k] = rstart + bs*j + fields_col[k]; 515 } 516 } 517 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 518 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)pc),nslots*nfields,jj,PETSC_OWN_POINTER,&ilink->is_col);CHKERRQ(ierr); 519 } else { 520 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 521 ierr = ISCreateStride(PetscObjectComm((PetscObject)pc),nslots,rstart+ilink->fields_col[0],bs,&ilink->is_col);CHKERRQ(ierr); 522 } 523 } 524 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 525 if (ilink->is_col) { ierr = ISSorted(ilink->is_col,&sorted_col);CHKERRQ(ierr); } 526 if (!sorted || !sorted_col) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 527 ilink = ilink->next; 528 } 529 } 530 531 ilink = jac->head; 532 if (!jac->pmat) { 533 Vec xtmp; 534 535 ierr = MatGetVecs(pc->pmat,&xtmp,NULL);CHKERRQ(ierr); 536 ierr = PetscMalloc1(nsplit,&jac->pmat);CHKERRQ(ierr); 537 ierr = PetscMalloc2(nsplit,&jac->x,nsplit,&jac->y);CHKERRQ(ierr); 538 for (i=0; i<nsplit; i++) { 539 MatNullSpace sp; 540 541 /* Check for preconditioning matrix attached to IS */ 542 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &jac->pmat[i]);CHKERRQ(ierr); 543 if (jac->pmat[i]) { 544 ierr = PetscObjectReference((PetscObject) jac->pmat[i]);CHKERRQ(ierr); 545 if (jac->type == PC_COMPOSITE_SCHUR) { 546 jac->schur_user = jac->pmat[i]; 547 548 ierr = PetscObjectReference((PetscObject) jac->schur_user);CHKERRQ(ierr); 549 } 550 } else { 551 const char *prefix; 552 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 553 ierr = KSPGetOptionsPrefix(ilink->ksp,&prefix);CHKERRQ(ierr); 554 ierr = MatSetOptionsPrefix(jac->pmat[i],prefix);CHKERRQ(ierr); 555 ierr = MatViewFromOptions(jac->pmat[i],NULL,"-mat_view");CHKERRQ(ierr); 556 } 557 /* create work vectors for each split */ 558 ierr = MatGetVecs(jac->pmat[i],&jac->x[i],&jac->y[i]);CHKERRQ(ierr); 559 ilink->x = jac->x[i]; ilink->y = jac->y[i]; ilink->z = NULL; 560 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 561 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],NULL,&ilink->sctx);CHKERRQ(ierr); 562 /* Check for null space attached to IS */ 563 ierr = PetscObjectQuery((PetscObject) ilink->is, "nullspace", (PetscObject*) &sp);CHKERRQ(ierr); 564 if (sp) { 565 ierr = MatSetNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 566 } 567 ierr = PetscObjectQuery((PetscObject) ilink->is, "nearnullspace", (PetscObject*) &sp);CHKERRQ(ierr); 568 if (sp) { 569 ierr = MatSetNearNullSpace(jac->pmat[i], sp);CHKERRQ(ierr); 570 } 571 ilink = ilink->next; 572 } 573 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 574 } else { 575 for (i=0; i<nsplit; i++) { 576 Mat pmat; 577 578 /* Check for preconditioning matrix attached to IS */ 579 ierr = PetscObjectQuery((PetscObject) ilink->is, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr); 580 if (!pmat) { 581 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 582 } 583 ilink = ilink->next; 584 } 585 } 586 if (pc->useAmat) { 587 ilink = jac->head; 588 if (!jac->mat) { 589 ierr = PetscMalloc1(nsplit,&jac->mat);CHKERRQ(ierr); 590 for (i=0; i<nsplit; i++) { 591 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 592 ilink = ilink->next; 593 } 594 } else { 595 for (i=0; i<nsplit; i++) { 596 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is_col,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 597 ilink = ilink->next; 598 } 599 } 600 } else { 601 jac->mat = jac->pmat; 602 } 603 604 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 605 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 606 ilink = jac->head; 607 if (!jac->Afield) { 608 ierr = PetscMalloc1(nsplit,&jac->Afield);CHKERRQ(ierr); 609 for (i=0; i<nsplit; i++) { 610 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 611 ilink = ilink->next; 612 } 613 } else { 614 for (i=0; i<nsplit; i++) { 615 ierr = MatGetSubMatrix(pc->mat,ilink->is,NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 616 ilink = ilink->next; 617 } 618 } 619 } 620 621 if (jac->type == PC_COMPOSITE_SCHUR) { 622 IS ccis; 623 PetscInt rstart,rend; 624 char lscname[256]; 625 PetscObject LSC_L; 626 627 if (nsplit != 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 628 629 /* When extracting off-diagonal submatrices, we take complements from this range */ 630 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 631 632 /* need to handle case when one is resetting up the preconditioner */ 633 if (jac->schur) { 634 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 635 636 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 637 ilink = jac->head; 638 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 639 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 640 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 641 ilink = ilink->next; 642 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 643 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 644 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 645 ierr = MatSchurComplementUpdateSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],pc->flag);CHKERRQ(ierr); 646 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 647 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 648 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 649 } 650 if (kspA != kspInner) { 651 ierr = KSPSetOperators(kspA,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 652 } 653 if (kspUpper != kspA) { 654 ierr = KSPSetOperators(kspUpper,jac->mat[0],jac->pmat[0],pc->flag);CHKERRQ(ierr); 655 } 656 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 657 } else { 658 const char *Dprefix; 659 char schurprefix[256], schurmatprefix[256]; 660 char schurtestoption[256]; 661 MatNullSpace sp; 662 PetscBool flg; 663 664 /* extract the A01 and A10 matrices */ 665 ilink = jac->head; 666 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 667 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 668 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 669 ilink = ilink->next; 670 ierr = ISComplement(ilink->is_col,rstart,rend,&ccis);CHKERRQ(ierr); 671 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 672 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 673 674 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 675 ierr = MatCreate(((PetscObject)jac->mat[0])->comm,&jac->schur);CHKERRQ(ierr); 676 ierr = MatSetType(jac->schur,MATSCHURCOMPLEMENT);CHKERRQ(ierr); 677 ierr = MatSchurComplementSetSubMatrices(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1]);CHKERRQ(ierr); 678 ierr = PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 679 /* 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? */ 680 ierr = MatSetOptionsPrefix(jac->schur,schurmatprefix);CHKERRQ(ierr); 681 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 682 ierr = MatGetNullSpace(jac->pmat[1], &sp);CHKERRQ(ierr); 683 if (sp) { 684 ierr = MatSetNullSpace(jac->schur, sp);CHKERRQ(ierr); 685 } 686 687 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname);CHKERRQ(ierr); 688 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 689 if (flg) { 690 DM dmInner; 691 KSP kspInner; 692 693 ierr = MatSchurComplementGetKSP(jac->schur, &kspInner);CHKERRQ(ierr); 694 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 695 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 696 ierr = PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject) pc, 2);CHKERRQ(ierr); 697 ierr = KSPSetOptionsPrefix(kspInner, schurprefix);CHKERRQ(ierr); 698 699 /* Set DM for new solver */ 700 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 701 ierr = KSPSetDM(kspInner, dmInner);CHKERRQ(ierr); 702 ierr = KSPSetDMActive(kspInner, PETSC_FALSE);CHKERRQ(ierr); 703 } else { 704 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 705 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 706 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 707 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 708 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 709 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 710 ierr = KSPSetType(jac->head->ksp,KSPGMRES);CHKERRQ(ierr); 711 ierr = MatSchurComplementSetKSP(jac->schur,jac->head->ksp);CHKERRQ(ierr); 712 } 713 ierr = KSPSetOperators(jac->head->ksp,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 714 ierr = KSPSetFromOptions(jac->head->ksp);CHKERRQ(ierr); 715 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 716 717 ierr = PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname);CHKERRQ(ierr); 718 ierr = PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->prefix, schurtestoption, NULL, &flg);CHKERRQ(ierr); 719 if (flg) { 720 DM dmInner; 721 722 ierr = PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname);CHKERRQ(ierr); 723 ierr = KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper);CHKERRQ(ierr); 724 ierr = KSPSetOptionsPrefix(jac->kspupper, schurprefix);CHKERRQ(ierr); 725 ierr = KSPGetDM(jac->head->ksp, &dmInner);CHKERRQ(ierr); 726 ierr = KSPSetDM(jac->kspupper, dmInner);CHKERRQ(ierr); 727 ierr = KSPSetDMActive(jac->kspupper, PETSC_FALSE);CHKERRQ(ierr); 728 ierr = KSPSetFromOptions(jac->kspupper);CHKERRQ(ierr); 729 ierr = KSPSetOperators(jac->kspupper,jac->mat[0],jac->pmat[0],flag);CHKERRQ(ierr); 730 ierr = VecDuplicate(jac->head->x, &jac->head->z);CHKERRQ(ierr); 731 } else { 732 jac->kspupper = jac->head->ksp; 733 ierr = PetscObjectReference((PetscObject) jac->head->ksp);CHKERRQ(ierr); 734 } 735 736 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 737 ierr = MatSchurComplementGetPmat(jac->schur,MAT_INITIAL_MATRIX,&jac->schurp);CHKERRQ(ierr); 738 } 739 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&jac->kspschur);CHKERRQ(ierr); 740 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 741 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 742 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 743 PC pcschur; 744 ierr = KSPGetPC(jac->kspschur,&pcschur);CHKERRQ(ierr); 745 ierr = PCSetType(pcschur,PCNONE);CHKERRQ(ierr); 746 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 747 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 748 ierr = MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user);CHKERRQ(ierr); 749 } 750 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 751 ierr = KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix);CHKERRQ(ierr); 752 ierr = KSPSetOptionsPrefix(jac->kspschur, Dprefix);CHKERRQ(ierr); 753 /* propogate DM */ 754 { 755 DM sdm; 756 ierr = KSPGetDM(jac->head->next->ksp, &sdm);CHKERRQ(ierr); 757 if (sdm) { 758 ierr = KSPSetDM(jac->kspschur, sdm);CHKERRQ(ierr); 759 ierr = KSPSetDMActive(jac->kspschur, PETSC_FALSE);CHKERRQ(ierr); 760 } 761 } 762 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 763 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 764 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 765 } 766 767 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 768 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_L",ilink->splitname);CHKERRQ(ierr); 769 ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 770 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 771 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_L",(PetscObject)LSC_L);CHKERRQ(ierr);} 772 ierr = PetscSNPrintf(lscname,sizeof(lscname),"%s_LSC_Lp",ilink->splitname);CHKERRQ(ierr); 773 ierr = PetscObjectQuery((PetscObject)pc->pmat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr); 774 if (!LSC_L) {ierr = PetscObjectQuery((PetscObject)pc->mat,lscname,(PetscObject*)&LSC_L);CHKERRQ(ierr);} 775 if (LSC_L) {ierr = PetscObjectCompose((PetscObject)jac->schur,"LSC_Lp",(PetscObject)LSC_L);CHKERRQ(ierr);} 776 } else { 777 /* set up the individual splits' PCs */ 778 i = 0; 779 ilink = jac->head; 780 while (ilink) { 781 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 782 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 783 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 784 i++; 785 ilink = ilink->next; 786 } 787 } 788 789 jac->suboptionsset = PETSC_TRUE; 790 PetscFunctionReturn(0); 791 } 792 793 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 794 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 795 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 796 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 797 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 798 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 799 800 #undef __FUNCT__ 801 #define __FUNCT__ "PCApply_FieldSplit_Schur" 802 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 803 { 804 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 805 PetscErrorCode ierr; 806 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 807 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 808 809 PetscFunctionBegin; 810 switch (jac->schurfactorization) { 811 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 812 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 813 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 814 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 815 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 816 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 817 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 818 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 819 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 820 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 821 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 822 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 823 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 824 break; 825 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 826 /* [A00 0; A10 S], suitable for left preconditioning */ 827 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 828 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 829 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 830 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 831 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 832 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 833 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 834 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 835 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 836 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 837 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 838 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 839 break; 840 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 841 /* [A00 A01; 0 S], suitable for right preconditioning */ 842 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 843 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 844 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 845 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 846 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 847 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 848 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 849 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 850 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 851 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 852 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 853 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 854 break; 855 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 856 /* [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 */ 857 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 858 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 859 ierr = KSPSolve(kspLower,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 860 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 861 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 862 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 863 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 864 865 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 866 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 867 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 868 869 if (kspUpper == kspA) { 870 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 871 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 872 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 873 } else { 874 ierr = KSPSolve(kspA,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 875 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 876 ierr = KSPSolve(kspUpper,ilinkA->x,ilinkA->z);CHKERRQ(ierr); 877 ierr = VecAXPY(ilinkA->y,-1.0,ilinkA->z);CHKERRQ(ierr); 878 } 879 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 880 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 881 } 882 PetscFunctionReturn(0); 883 } 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "PCApply_FieldSplit" 887 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 888 { 889 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 890 PetscErrorCode ierr; 891 PC_FieldSplitLink ilink = jac->head; 892 PetscInt cnt,bs; 893 894 PetscFunctionBegin; 895 if (jac->type == PC_COMPOSITE_ADDITIVE) { 896 if (jac->defaultsplit) { 897 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 898 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); 899 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 900 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); 901 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 902 while (ilink) { 903 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 904 ilink = ilink->next; 905 } 906 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 907 } else { 908 ierr = VecSet(y,0.0);CHKERRQ(ierr); 909 while (ilink) { 910 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 911 ilink = ilink->next; 912 } 913 } 914 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 915 if (!jac->w1) { 916 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 917 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 918 } 919 ierr = VecSet(y,0.0);CHKERRQ(ierr); 920 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 921 cnt = 1; 922 while (ilink->next) { 923 ilink = ilink->next; 924 /* compute the residual only over the part of the vector needed */ 925 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 926 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 927 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 928 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 929 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 930 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 931 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 932 } 933 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 934 cnt -= 2; 935 while (ilink->previous) { 936 ilink = ilink->previous; 937 /* compute the residual only over the part of the vector needed */ 938 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 939 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 940 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 941 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 942 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 943 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 944 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 945 } 946 } 947 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 948 PetscFunctionReturn(0); 949 } 950 951 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 952 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 953 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 954 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 955 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 956 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 957 958 #undef __FUNCT__ 959 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 960 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 961 { 962 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 963 PetscErrorCode ierr; 964 PC_FieldSplitLink ilink = jac->head; 965 PetscInt bs; 966 967 PetscFunctionBegin; 968 if (jac->type == PC_COMPOSITE_ADDITIVE) { 969 if (jac->defaultsplit) { 970 ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr); 971 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); 972 ierr = VecGetBlockSize(y,&bs);CHKERRQ(ierr); 973 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); 974 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 975 while (ilink) { 976 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 977 ilink = ilink->next; 978 } 979 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 980 } else { 981 ierr = VecSet(y,0.0);CHKERRQ(ierr); 982 while (ilink) { 983 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 984 ilink = ilink->next; 985 } 986 } 987 } else { 988 if (!jac->w1) { 989 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 990 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 991 } 992 ierr = VecSet(y,0.0);CHKERRQ(ierr); 993 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 994 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 995 while (ilink->next) { 996 ilink = ilink->next; 997 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 998 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 999 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1000 } 1001 while (ilink->previous) { 1002 ilink = ilink->previous; 1003 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1004 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1005 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1006 } 1007 } else { 1008 while (ilink->next) { /* get to last entry in linked list */ 1009 ilink = ilink->next; 1010 } 1011 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 1012 while (ilink->previous) { 1013 ilink = ilink->previous; 1014 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 1015 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 1016 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 1017 } 1018 } 1019 } 1020 PetscFunctionReturn(0); 1021 } 1022 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "PCReset_FieldSplit" 1025 static PetscErrorCode PCReset_FieldSplit(PC pc) 1026 { 1027 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1028 PetscErrorCode ierr; 1029 PC_FieldSplitLink ilink = jac->head,next; 1030 1031 PetscFunctionBegin; 1032 while (ilink) { 1033 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 1034 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 1035 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 1036 ierr = VecDestroy(&ilink->z);CHKERRQ(ierr); 1037 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 1038 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1039 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1040 next = ilink->next; 1041 ilink = next; 1042 } 1043 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1044 if (jac->mat && jac->mat != jac->pmat) { 1045 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 1046 } else if (jac->mat) { 1047 jac->mat = NULL; 1048 } 1049 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 1050 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 1051 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 1052 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 1053 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 1054 ierr = MatDestroy(&jac->schurp);CHKERRQ(ierr); 1055 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1056 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 1057 ierr = KSPDestroy(&jac->kspupper);CHKERRQ(ierr); 1058 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 1059 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1060 jac->reset = PETSC_TRUE; 1061 PetscFunctionReturn(0); 1062 } 1063 1064 #undef __FUNCT__ 1065 #define __FUNCT__ "PCDestroy_FieldSplit" 1066 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1067 { 1068 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1069 PetscErrorCode ierr; 1070 PC_FieldSplitLink ilink = jac->head,next; 1071 1072 PetscFunctionBegin; 1073 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 1074 while (ilink) { 1075 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 1076 next = ilink->next; 1077 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 1078 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 1079 ierr = PetscFree(ilink->fields_col);CHKERRQ(ierr); 1080 ierr = PetscFree(ilink);CHKERRQ(ierr); 1081 ilink = next; 1082 } 1083 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 1084 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1085 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",NULL);CHKERRQ(ierr); 1086 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",NULL);CHKERRQ(ierr); 1087 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",NULL);CHKERRQ(ierr); 1088 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",NULL);CHKERRQ(ierr); 1089 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",NULL);CHKERRQ(ierr); 1090 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",NULL);CHKERRQ(ierr); 1091 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",NULL);CHKERRQ(ierr); 1092 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",NULL);CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 #undef __FUNCT__ 1097 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 1098 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 1099 { 1100 PetscErrorCode ierr; 1101 PetscInt bs; 1102 PetscBool flg,stokes = PETSC_FALSE; 1103 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1104 PCCompositeType ctype; 1105 1106 PetscFunctionBegin; 1107 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 1108 ierr = PetscOptionsBool("-pc_fieldsplit_dm_splits","Whether to use DMCreateFieldDecomposition() for splits","PCFieldSplitSetDMSplits",jac->dm_splits,&jac->dm_splits,NULL);CHKERRQ(ierr); 1109 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 1110 if (flg) { 1111 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 1112 } 1113 1114 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,NULL);CHKERRQ(ierr); 1115 if (stokes) { 1116 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 1117 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 1118 } 1119 1120 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 1121 if (flg) { 1122 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 1123 } 1124 /* Only setup fields once */ 1125 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1126 /* only allow user to set fields from command line if bs is already known. 1127 otherwise user can set them in PCFieldSplitSetDefaults() */ 1128 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 1129 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 1130 } 1131 if (jac->type == PC_COMPOSITE_SCHUR) { 1132 ierr = PetscOptionsGetEnum(((PetscObject)pc)->prefix,"-pc_fieldsplit_schur_factorization_type",PCFieldSplitSchurFactTypes,(PetscEnum*)&jac->schurfactorization,&flg);CHKERRQ(ierr); 1133 if (flg) {ierr = PetscInfo(pc,"Deprecated use of -pc_fieldsplit_schur_factorization_type\n");CHKERRQ(ierr);} 1134 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); 1135 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSetSchurPre",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,NULL);CHKERRQ(ierr); 1136 } 1137 ierr = PetscOptionsTail();CHKERRQ(ierr); 1138 PetscFunctionReturn(0); 1139 } 1140 1141 /*------------------------------------------------------------------------------------*/ 1142 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 1145 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1146 { 1147 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1148 PetscErrorCode ierr; 1149 PC_FieldSplitLink ilink,next = jac->head; 1150 char prefix[128]; 1151 PetscInt i; 1152 1153 PetscFunctionBegin; 1154 if (jac->splitdefined) { 1155 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1156 PetscFunctionReturn(0); 1157 } 1158 for (i=0; i<n; i++) { 1159 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 1160 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 1161 } 1162 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1163 if (splitname) { 1164 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1165 } else { 1166 ierr = PetscMalloc1(3,&ilink->splitname);CHKERRQ(ierr); 1167 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 1168 } 1169 ierr = PetscMalloc1(n,&ilink->fields);CHKERRQ(ierr); 1170 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 1171 ierr = PetscMalloc1(n,&ilink->fields_col);CHKERRQ(ierr); 1172 ierr = PetscMemcpy(ilink->fields_col,fields_col,n*sizeof(PetscInt));CHKERRQ(ierr); 1173 1174 ilink->nfields = n; 1175 ilink->next = NULL; 1176 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1177 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1178 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1179 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1180 1181 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1182 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1183 1184 if (!next) { 1185 jac->head = ilink; 1186 ilink->previous = NULL; 1187 } else { 1188 while (next->next) { 1189 next = next->next; 1190 } 1191 next->next = ilink; 1192 ilink->previous = next; 1193 } 1194 jac->nsplits++; 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1200 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1201 { 1202 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1203 PetscErrorCode ierr; 1204 1205 PetscFunctionBegin; 1206 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1207 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1208 1209 (*subksp)[1] = jac->kspschur; 1210 if (n) *n = jac->nsplits; 1211 PetscFunctionReturn(0); 1212 } 1213 1214 #undef __FUNCT__ 1215 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1216 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1217 { 1218 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1219 PetscErrorCode ierr; 1220 PetscInt cnt = 0; 1221 PC_FieldSplitLink ilink = jac->head; 1222 1223 PetscFunctionBegin; 1224 ierr = PetscMalloc1(jac->nsplits,subksp);CHKERRQ(ierr); 1225 while (ilink) { 1226 (*subksp)[cnt++] = ilink->ksp; 1227 ilink = ilink->next; 1228 } 1229 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); 1230 if (n) *n = jac->nsplits; 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNCT__ 1235 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1236 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1237 { 1238 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1239 PetscErrorCode ierr; 1240 PC_FieldSplitLink ilink, next = jac->head; 1241 char prefix[128]; 1242 1243 PetscFunctionBegin; 1244 if (jac->splitdefined) { 1245 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1246 PetscFunctionReturn(0); 1247 } 1248 ierr = PetscNew(&ilink);CHKERRQ(ierr); 1249 if (splitname) { 1250 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1251 } else { 1252 ierr = PetscMalloc1(8,&ilink->splitname);CHKERRQ(ierr); 1253 ierr = PetscSNPrintf(ilink->splitname,7,"%D",jac->nsplits);CHKERRQ(ierr); 1254 } 1255 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1256 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 1257 ilink->is = is; 1258 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1259 ierr = ISDestroy(&ilink->is_col);CHKERRQ(ierr); 1260 ilink->is_col = is; 1261 ilink->next = NULL; 1262 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&ilink->ksp);CHKERRQ(ierr); 1263 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1264 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1265 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1266 1267 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sfieldsplit_%s_",((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "",ilink->splitname);CHKERRQ(ierr); 1268 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1269 1270 if (!next) { 1271 jac->head = ilink; 1272 ilink->previous = NULL; 1273 } else { 1274 while (next->next) { 1275 next = next->next; 1276 } 1277 next->next = ilink; 1278 ilink->previous = next; 1279 } 1280 jac->nsplits++; 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "PCFieldSplitSetFields" 1286 /*@ 1287 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1288 1289 Logically Collective on PC 1290 1291 Input Parameters: 1292 + pc - the preconditioner context 1293 . splitname - name of this split, if NULL the number of the split is used 1294 . n - the number of fields in this split 1295 - fields - the fields in this split 1296 1297 Level: intermediate 1298 1299 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1300 1301 The PCFieldSplitSetFields() is for defining fields as strided blocks. For example, if the block 1302 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 1303 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1304 where the numbered entries indicate what is in the field. 1305 1306 This function is called once per split (it creates a new split each time). Solve options 1307 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1308 1309 Developer Note: This routine does not actually create the IS representing the split, that is delayed 1310 until PCSetUp_FieldSplit(), because information about the vector/matrix layouts may not be 1311 available when this routine is called. 1312 1313 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1314 1315 @*/ 1316 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields,const PetscInt *fields_col) 1317 { 1318 PetscErrorCode ierr; 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1322 PetscValidCharPointer(splitname,2); 1323 if (n < 1) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1324 PetscValidIntPointer(fields,3); 1325 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt*,const PetscInt*),(pc,splitname,n,fields,fields_col));CHKERRQ(ierr); 1326 PetscFunctionReturn(0); 1327 } 1328 1329 #undef __FUNCT__ 1330 #define __FUNCT__ "PCFieldSplitSetIS" 1331 /*@ 1332 PCFieldSplitSetIS - Sets the exact elements for field 1333 1334 Logically Collective on PC 1335 1336 Input Parameters: 1337 + pc - the preconditioner context 1338 . splitname - name of this split, if NULL the number of the split is used 1339 - is - the index set that defines the vector elements in this field 1340 1341 1342 Notes: 1343 Use PCFieldSplitSetFields(), for fields defined by strided types. 1344 1345 This function is called once per split (it creates a new split each time). Solve options 1346 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1347 1348 Level: intermediate 1349 1350 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1351 1352 @*/ 1353 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1354 { 1355 PetscErrorCode ierr; 1356 1357 PetscFunctionBegin; 1358 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1359 if (splitname) PetscValidCharPointer(splitname,2); 1360 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1361 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1362 PetscFunctionReturn(0); 1363 } 1364 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "PCFieldSplitGetIS" 1367 /*@ 1368 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1369 1370 Logically Collective on PC 1371 1372 Input Parameters: 1373 + pc - the preconditioner context 1374 - splitname - name of this split 1375 1376 Output Parameter: 1377 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1378 1379 Level: intermediate 1380 1381 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1382 1383 @*/ 1384 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1385 { 1386 PetscErrorCode ierr; 1387 1388 PetscFunctionBegin; 1389 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1390 PetscValidCharPointer(splitname,2); 1391 PetscValidPointer(is,3); 1392 { 1393 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1394 PC_FieldSplitLink ilink = jac->head; 1395 PetscBool found; 1396 1397 *is = NULL; 1398 while (ilink) { 1399 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1400 if (found) { 1401 *is = ilink->is; 1402 break; 1403 } 1404 ilink = ilink->next; 1405 } 1406 } 1407 PetscFunctionReturn(0); 1408 } 1409 1410 #undef __FUNCT__ 1411 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1412 /*@ 1413 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1414 fieldsplit preconditioner. If not set the matrix block size is used. 1415 1416 Logically Collective on PC 1417 1418 Input Parameters: 1419 + pc - the preconditioner context 1420 - bs - the block size 1421 1422 Level: intermediate 1423 1424 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1425 1426 @*/ 1427 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1428 { 1429 PetscErrorCode ierr; 1430 1431 PetscFunctionBegin; 1432 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1433 PetscValidLogicalCollectiveInt(pc,bs,2); 1434 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1435 PetscFunctionReturn(0); 1436 } 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1440 /*@C 1441 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1442 1443 Collective on KSP 1444 1445 Input Parameter: 1446 . pc - the preconditioner context 1447 1448 Output Parameters: 1449 + n - the number of splits 1450 - pc - the array of KSP contexts 1451 1452 Note: 1453 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1454 (not the KSP just the array that contains them). 1455 1456 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1457 1458 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1459 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1460 KSP array must be. 1461 1462 1463 Level: advanced 1464 1465 .seealso: PCFIELDSPLIT 1466 @*/ 1467 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1468 { 1469 PetscErrorCode ierr; 1470 1471 PetscFunctionBegin; 1472 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1473 if (n) PetscValidIntPointer(n,2); 1474 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1475 PetscFunctionReturn(0); 1476 } 1477 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "PCFieldSplitSetSchurPre" 1480 /*@ 1481 PCFieldSplitSetSchurPre - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1482 A11 matrix. Otherwise no preconditioner is used. 1483 1484 Collective on PC 1485 1486 Input Parameters: 1487 + pc - the preconditioner context 1488 . 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 1489 - userpre - matrix to use for preconditioning, or NULL 1490 1491 Options Database: 1492 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> default is a11 1493 1494 Notes: 1495 $ If ptype is 1496 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1497 $ to this function). 1498 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1499 $ matrix associated with the Schur complement (i.e. A11) 1500 $ full then the preconditioner uses the exact Schur complement (this is expensive) 1501 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1502 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1503 $ preconditioner 1504 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1505 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. 1506 1507 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1508 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1509 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1510 1511 Level: intermediate 1512 1513 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1514 1515 @*/ 1516 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1517 { 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1522 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1523 PetscFunctionReturn(0); 1524 } 1525 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1526 1527 #undef __FUNCT__ 1528 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1529 /*@ 1530 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1531 preconditioned. See PCFieldSplitSetSchurPre() for details. 1532 1533 Logically Collective on PC 1534 1535 Input Parameters: 1536 . pc - the preconditioner context 1537 1538 Output Parameters: 1539 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1540 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1541 1542 Level: intermediate 1543 1544 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1545 1546 @*/ 1547 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1548 { 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1553 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "PCFieldSplitSchurGetS" 1559 /*@ 1560 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1561 1562 Not collective 1563 1564 Input Parameter: 1565 . pc - the preconditioner context 1566 1567 Output Parameter: 1568 . S - the Schur complement matrix 1569 1570 Notes: 1571 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1572 1573 Level: advanced 1574 1575 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1576 1577 @*/ 1578 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1579 { 1580 PetscErrorCode ierr; 1581 const char* t; 1582 PetscBool isfs; 1583 PC_FieldSplit *jac; 1584 1585 PetscFunctionBegin; 1586 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1587 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1588 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1589 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1590 jac = (PC_FieldSplit*)pc->data; 1591 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1592 if (S) *S = jac->schur; 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1598 /*@ 1599 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1600 1601 Not collective 1602 1603 Input Parameters: 1604 + pc - the preconditioner context 1605 . S - the Schur complement matrix 1606 1607 Level: advanced 1608 1609 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 1610 1611 @*/ 1612 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1613 { 1614 PetscErrorCode ierr; 1615 const char* t; 1616 PetscBool isfs; 1617 PC_FieldSplit *jac; 1618 1619 PetscFunctionBegin; 1620 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1621 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1622 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1623 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1624 jac = (PC_FieldSplit*)pc->data; 1625 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1626 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1627 PetscFunctionReturn(0); 1628 } 1629 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 1633 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1634 { 1635 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1636 PetscErrorCode ierr; 1637 1638 PetscFunctionBegin; 1639 jac->schurpre = ptype; 1640 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1641 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1642 jac->schur_user = pre; 1643 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1644 } 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #undef __FUNCT__ 1649 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 1650 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1651 { 1652 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1653 1654 PetscFunctionBegin; 1655 *ptype = jac->schurpre; 1656 *pre = jac->schur_user; 1657 PetscFunctionReturn(0); 1658 } 1659 1660 #undef __FUNCT__ 1661 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1662 /*@ 1663 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1664 1665 Collective on PC 1666 1667 Input Parameters: 1668 + pc - the preconditioner context 1669 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1670 1671 Options Database: 1672 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1673 1674 1675 Level: intermediate 1676 1677 Notes: 1678 The FULL factorization is 1679 1680 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1681 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1682 1683 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, 1684 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). 1685 1686 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 1687 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 1688 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, 1689 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1690 1691 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 1692 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). 1693 1694 References: 1695 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1696 1697 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1698 1699 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1700 @*/ 1701 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1702 { 1703 PetscErrorCode ierr; 1704 1705 PetscFunctionBegin; 1706 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1707 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1713 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1714 { 1715 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1716 1717 PetscFunctionBegin; 1718 jac->schurfactorization = ftype; 1719 PetscFunctionReturn(0); 1720 } 1721 1722 #undef __FUNCT__ 1723 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1724 /*@C 1725 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1726 1727 Collective on KSP 1728 1729 Input Parameter: 1730 . pc - the preconditioner context 1731 1732 Output Parameters: 1733 + A00 - the (0,0) block 1734 . A01 - the (0,1) block 1735 . A10 - the (1,0) block 1736 - A11 - the (1,1) block 1737 1738 Level: advanced 1739 1740 .seealso: PCFIELDSPLIT 1741 @*/ 1742 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1743 { 1744 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1745 1746 PetscFunctionBegin; 1747 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1748 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1749 if (A00) *A00 = jac->pmat[0]; 1750 if (A01) *A01 = jac->B; 1751 if (A10) *A10 = jac->C; 1752 if (A11) *A11 = jac->pmat[1]; 1753 PetscFunctionReturn(0); 1754 } 1755 1756 #undef __FUNCT__ 1757 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1758 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1759 { 1760 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1761 PetscErrorCode ierr; 1762 1763 PetscFunctionBegin; 1764 jac->type = type; 1765 if (type == PC_COMPOSITE_SCHUR) { 1766 pc->ops->apply = PCApply_FieldSplit_Schur; 1767 pc->ops->view = PCView_FieldSplit_Schur; 1768 1769 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1770 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 1771 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 1772 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1773 1774 } else { 1775 pc->ops->apply = PCApply_FieldSplit; 1776 pc->ops->view = PCView_FieldSplit; 1777 1778 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1779 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 1781 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1782 } 1783 PetscFunctionReturn(0); 1784 } 1785 1786 #undef __FUNCT__ 1787 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1788 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1789 { 1790 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1791 1792 PetscFunctionBegin; 1793 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1794 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); 1795 jac->bs = bs; 1796 PetscFunctionReturn(0); 1797 } 1798 1799 #undef __FUNCT__ 1800 #define __FUNCT__ "PCFieldSplitSetType" 1801 /*@ 1802 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1803 1804 Collective on PC 1805 1806 Input Parameter: 1807 . pc - the preconditioner context 1808 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1809 1810 Options Database Key: 1811 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1812 1813 Level: Intermediate 1814 1815 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1816 1817 .seealso: PCCompositeSetType() 1818 1819 @*/ 1820 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1821 { 1822 PetscErrorCode ierr; 1823 1824 PetscFunctionBegin; 1825 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1826 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1827 PetscFunctionReturn(0); 1828 } 1829 1830 #undef __FUNCT__ 1831 #define __FUNCT__ "PCFieldSplitGetType" 1832 /*@ 1833 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1834 1835 Not collective 1836 1837 Input Parameter: 1838 . pc - the preconditioner context 1839 1840 Output Parameter: 1841 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1842 1843 Level: Intermediate 1844 1845 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1846 .seealso: PCCompositeSetType() 1847 @*/ 1848 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1849 { 1850 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1851 1852 PetscFunctionBegin; 1853 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1854 PetscValidIntPointer(type,2); 1855 *type = jac->type; 1856 PetscFunctionReturn(0); 1857 } 1858 1859 #undef __FUNCT__ 1860 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1861 /*@ 1862 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1863 1864 Logically Collective 1865 1866 Input Parameters: 1867 + pc - the preconditioner context 1868 - flg - boolean indicating whether to use field splits defined by the DM 1869 1870 Options Database Key: 1871 . -pc_fieldsplit_dm_splits 1872 1873 Level: Intermediate 1874 1875 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1876 1877 .seealso: PCFieldSplitGetDMSplits() 1878 1879 @*/ 1880 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1881 { 1882 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1883 PetscBool isfs; 1884 PetscErrorCode ierr; 1885 1886 PetscFunctionBegin; 1887 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1888 PetscValidLogicalCollectiveBool(pc,flg,2); 1889 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1890 if (isfs) { 1891 jac->dm_splits = flg; 1892 } 1893 PetscFunctionReturn(0); 1894 } 1895 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1899 /*@ 1900 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1901 1902 Logically Collective 1903 1904 Input Parameter: 1905 . pc - the preconditioner context 1906 1907 Output Parameter: 1908 . flg - boolean indicating whether to use field splits defined by the DM 1909 1910 Level: Intermediate 1911 1912 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1913 1914 .seealso: PCFieldSplitSetDMSplits() 1915 1916 @*/ 1917 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 1918 { 1919 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1920 PetscBool isfs; 1921 PetscErrorCode ierr; 1922 1923 PetscFunctionBegin; 1924 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1925 PetscValidPointer(flg,2); 1926 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1927 if (isfs) { 1928 if(flg) *flg = jac->dm_splits; 1929 } 1930 PetscFunctionReturn(0); 1931 } 1932 1933 /* -------------------------------------------------------------------------------------*/ 1934 /*MC 1935 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1936 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1937 1938 To set options on the solvers for each block append -fieldsplit_ to all the PC 1939 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1940 1941 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1942 and set the options directly on the resulting KSP object 1943 1944 Level: intermediate 1945 1946 Options Database Keys: 1947 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1948 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1949 been supplied explicitly by -pc_fieldsplit_%d_fields 1950 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1951 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 1952 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11 1953 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1954 1955 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 1956 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 1957 1958 Notes: 1959 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1960 to define a field by an arbitrary collection of entries. 1961 1962 If no fields are set the default is used. The fields are defined by entries strided by bs, 1963 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1964 if this is not called the block size defaults to the blocksize of the second matrix passed 1965 to KSPSetOperators()/PCSetOperators(). 1966 1967 $ For the Schur complement preconditioner if J = ( A00 A01 ) 1968 $ ( A10 A11 ) 1969 $ the preconditioner using full factorization is 1970 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 1971 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 1972 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 1973 $ S = A11 - A10 ksp(A00) A01 1974 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 1975 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 1976 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1977 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() to turn on or off this 1978 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 1979 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 1980 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 1981 (e.g., -fieldsplit_1_pc_type asm). 1982 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 1983 diag gives 1984 $ ( inv(A00) 0 ) 1985 $ ( 0 -ksp(S) ) 1986 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 1987 $ ( A00 0 ) 1988 $ ( A10 S ) 1989 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 1990 $ ( A00 A01 ) 1991 $ ( 0 S ) 1992 where again the inverses of A00 and S are applied using KSPs. 1993 1994 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1995 is used automatically for a second block. 1996 1997 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1998 Generally it should be used with the AIJ format. 1999 2000 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2001 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2002 inside a smoother resulting in "Distributive Smoothers". 2003 2004 Concepts: physics based preconditioners, block preconditioners 2005 2006 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2007 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre() 2008 M*/ 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "PCCreate_FieldSplit" 2012 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2013 { 2014 PetscErrorCode ierr; 2015 PC_FieldSplit *jac; 2016 2017 PetscFunctionBegin; 2018 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2019 2020 jac->bs = -1; 2021 jac->nsplits = 0; 2022 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2023 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2024 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2025 jac->dm_splits = PETSC_TRUE; 2026 2027 pc->data = (void*)jac; 2028 2029 pc->ops->apply = PCApply_FieldSplit; 2030 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2031 pc->ops->setup = PCSetUp_FieldSplit; 2032 pc->ops->reset = PCReset_FieldSplit; 2033 pc->ops->destroy = PCDestroy_FieldSplit; 2034 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2035 pc->ops->view = PCView_FieldSplit; 2036 pc->ops->applyrichardson = 0; 2037 2038 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2039 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2040 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2041 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2042 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 2047 2048