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