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