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 (the lumped) 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 Options Database: 1339 . -pc_fieldsplit_diag_use_amat 1340 1341 Level: intermediate 1342 1343 .seealso: PCFieldSplitGetDiagUseAmat(), PCFieldSplitSetOffDiagUseAmat(), PCFIELDSPLIT 1344 1345 @*/ 1346 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc,PetscBool flg) 1347 { 1348 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1349 PetscBool isfs; 1350 PetscErrorCode ierr; 1351 1352 PetscFunctionBegin; 1353 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1354 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1355 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1356 jac->diag_use_amat = flg; 1357 PetscFunctionReturn(0); 1358 } 1359 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "PCFieldSplitGetDiagUseAmat" 1362 /*@ 1363 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) 1364 1365 Logically Collective on PC 1366 1367 Input Parameters: 1368 . pc - the preconditioner object 1369 1370 Output Parameters: 1371 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 1372 1373 1374 Level: intermediate 1375 1376 .seealso: PCFieldSplitSetDiagUseAmat(), PCFieldSplitGetOffDiagUseAmat(), PCFIELDSPLIT 1377 1378 @*/ 1379 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc,PetscBool *flg) 1380 { 1381 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1382 PetscBool isfs; 1383 PetscErrorCode ierr; 1384 1385 PetscFunctionBegin; 1386 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1387 PetscValidPointer(flg,2); 1388 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1389 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1390 *flg = jac->diag_use_amat; 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "PCFieldSplitSetOffDiagUseAmat" 1396 /*@ 1397 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1398 1399 Logically Collective on PC 1400 1401 Input Parameters: 1402 + pc - the preconditioner object 1403 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1404 1405 Options Database: 1406 . -pc_fieldsplit_off_diag_use_amat 1407 1408 Level: intermediate 1409 1410 .seealso: PCFieldSplitGetOffDiagUseAmat(), PCFieldSplitSetDiagUseAmat(), PCFIELDSPLIT 1411 1412 @*/ 1413 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc,PetscBool flg) 1414 { 1415 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1416 PetscBool isfs; 1417 PetscErrorCode ierr; 1418 1419 PetscFunctionBegin; 1420 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1421 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1422 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1423 jac->offdiag_use_amat = flg; 1424 PetscFunctionReturn(0); 1425 } 1426 1427 #undef __FUNCT__ 1428 #define __FUNCT__ "PCFieldSplitGetOffDiagUseAmat" 1429 /*@ 1430 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) 1431 1432 Logically Collective on PC 1433 1434 Input Parameters: 1435 . pc - the preconditioner object 1436 1437 Output Parameters: 1438 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 1439 1440 1441 Level: intermediate 1442 1443 .seealso: PCFieldSplitSetOffDiagUseAmat(), PCFieldSplitGetDiagUseAmat(), PCFIELDSPLIT 1444 1445 @*/ 1446 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc,PetscBool *flg) 1447 { 1448 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1449 PetscBool isfs; 1450 PetscErrorCode ierr; 1451 1452 PetscFunctionBegin; 1453 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1454 PetscValidPointer(flg,2); 1455 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1456 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"PC not of type %s",PCFIELDSPLIT); 1457 *flg = jac->offdiag_use_amat; 1458 PetscFunctionReturn(0); 1459 } 1460 1461 1462 1463 #undef __FUNCT__ 1464 #define __FUNCT__ "PCFieldSplitSetIS" 1465 /*@ 1466 PCFieldSplitSetIS - Sets the exact elements for field 1467 1468 Logically Collective on PC 1469 1470 Input Parameters: 1471 + pc - the preconditioner context 1472 . splitname - name of this split, if NULL the number of the split is used 1473 - is - the index set that defines the vector elements in this field 1474 1475 1476 Notes: 1477 Use PCFieldSplitSetFields(), for fields defined by strided types. 1478 1479 This function is called once per split (it creates a new split each time). Solve options 1480 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1481 1482 Level: intermediate 1483 1484 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1485 1486 @*/ 1487 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1488 { 1489 PetscErrorCode ierr; 1490 1491 PetscFunctionBegin; 1492 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1493 if (splitname) PetscValidCharPointer(splitname,2); 1494 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1495 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1496 PetscFunctionReturn(0); 1497 } 1498 1499 #undef __FUNCT__ 1500 #define __FUNCT__ "PCFieldSplitGetIS" 1501 /*@ 1502 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1503 1504 Logically Collective on PC 1505 1506 Input Parameters: 1507 + pc - the preconditioner context 1508 - splitname - name of this split 1509 1510 Output Parameter: 1511 - is - the index set that defines the vector elements in this field, or NULL if the field is not found 1512 1513 Level: intermediate 1514 1515 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1516 1517 @*/ 1518 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1519 { 1520 PetscErrorCode ierr; 1521 1522 PetscFunctionBegin; 1523 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1524 PetscValidCharPointer(splitname,2); 1525 PetscValidPointer(is,3); 1526 { 1527 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1528 PC_FieldSplitLink ilink = jac->head; 1529 PetscBool found; 1530 1531 *is = NULL; 1532 while (ilink) { 1533 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1534 if (found) { 1535 *is = ilink->is; 1536 break; 1537 } 1538 ilink = ilink->next; 1539 } 1540 } 1541 PetscFunctionReturn(0); 1542 } 1543 1544 #undef __FUNCT__ 1545 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1546 /*@ 1547 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1548 fieldsplit preconditioner. If not set the matrix block size is used. 1549 1550 Logically Collective on PC 1551 1552 Input Parameters: 1553 + pc - the preconditioner context 1554 - bs - the block size 1555 1556 Level: intermediate 1557 1558 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1559 1560 @*/ 1561 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1562 { 1563 PetscErrorCode ierr; 1564 1565 PetscFunctionBegin; 1566 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1567 PetscValidLogicalCollectiveInt(pc,bs,2); 1568 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 #undef __FUNCT__ 1573 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1574 /*@C 1575 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1576 1577 Collective on KSP 1578 1579 Input Parameter: 1580 . pc - the preconditioner context 1581 1582 Output Parameters: 1583 + n - the number of splits 1584 - pc - the array of KSP contexts 1585 1586 Note: 1587 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1588 (not the KSP just the array that contains them). 1589 1590 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1591 1592 Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs. 1593 You can call PCFieldSplitGetSubKSP(pc,n,NULL_OBJECT,ierr) to determine how large the 1594 KSP array must be. 1595 1596 1597 Level: advanced 1598 1599 .seealso: PCFIELDSPLIT 1600 @*/ 1601 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1602 { 1603 PetscErrorCode ierr; 1604 1605 PetscFunctionBegin; 1606 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1607 if (n) PetscValidIntPointer(n,2); 1608 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1609 PetscFunctionReturn(0); 1610 } 1611 1612 #undef __FUNCT__ 1613 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1614 /*@ 1615 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1616 A11 matrix. Otherwise no preconditioner is used. 1617 1618 Collective on PC 1619 1620 Input Parameters: 1621 + pc - the preconditioner context 1622 . 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 1623 - userpre - matrix to use for preconditioning, or NULL 1624 1625 Options Database: 1626 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11> default is a11 1627 1628 Notes: 1629 $ If ptype is 1630 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1631 $ to this function). 1632 $ a11 then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1633 $ matrix associated with the Schur complement (i.e. A11) 1634 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1635 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1636 $ preconditioner 1637 $ selfp then the preconditioning matrix is an explicitly-assembled approximation Sp = A11 - A10 inv(diag(A00)) A01 1638 $ This is only a good preconditioner when diag(A00) is a good preconditioner for A00. 1639 1640 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1641 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1642 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1643 1644 Developer Notes: This is a terrible name, which gives no good indication of what the function does. 1645 Among other things, it should have 'Set' in the name, since it sets the type of matrix to use. 1646 1647 Level: intermediate 1648 1649 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1650 1651 @*/ 1652 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1653 { 1654 PetscErrorCode ierr; 1655 1656 PetscFunctionBegin; 1657 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1658 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1659 PetscFunctionReturn(0); 1660 } 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "PCFieldSplitSchurGetS" 1664 /*@ 1665 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1666 1667 Not collective 1668 1669 Input Parameter: 1670 . pc - the preconditioner context 1671 1672 Output Parameter: 1673 . S - the Schur complement matrix 1674 1675 Notes: 1676 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1677 1678 Level: advanced 1679 1680 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1681 1682 @*/ 1683 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1684 { 1685 PetscErrorCode ierr; 1686 const char* t; 1687 PetscBool isfs; 1688 PC_FieldSplit *jac; 1689 1690 PetscFunctionBegin; 1691 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1692 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1693 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1694 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1695 jac = (PC_FieldSplit*)pc->data; 1696 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1697 if (S) *S = jac->schur; 1698 PetscFunctionReturn(0); 1699 } 1700 1701 #undef __FUNCT__ 1702 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1703 /*@ 1704 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1705 1706 Not collective 1707 1708 Input Parameters: 1709 + pc - the preconditioner context 1710 . S - the Schur complement matrix 1711 1712 Level: advanced 1713 1714 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSchurPrecondition(), MatSchurComplement, PCFieldSplitSchurGetS() 1715 1716 @*/ 1717 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 1718 { 1719 PetscErrorCode ierr; 1720 const char* t; 1721 PetscBool isfs; 1722 PC_FieldSplit *jac; 1723 1724 PetscFunctionBegin; 1725 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1726 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1727 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1728 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1729 jac = (PC_FieldSplit*)pc->data; 1730 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1731 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 1732 PetscFunctionReturn(0); 1733 } 1734 1735 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1738 static PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1739 { 1740 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1741 PetscErrorCode ierr; 1742 1743 PetscFunctionBegin; 1744 jac->schurpre = ptype; 1745 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 1746 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1747 jac->schur_user = pre; 1748 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1749 } 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 1755 /*@ 1756 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 1757 1758 Collective on PC 1759 1760 Input Parameters: 1761 + pc - the preconditioner context 1762 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 1763 1764 Options Database: 1765 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 1766 1767 1768 Level: intermediate 1769 1770 Notes: 1771 The FULL factorization is 1772 1773 $ (A B) = (1 0) (A 0) (1 Ainv*B) 1774 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 1775 1776 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, 1777 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). 1778 1779 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 1780 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 1781 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, 1782 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 1783 1784 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 1785 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). 1786 1787 References: 1788 Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000) pp. 1969-1972. 1789 1790 Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001), pp. 1050-1051. 1791 1792 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 1793 @*/ 1794 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 1795 { 1796 PetscErrorCode ierr; 1797 1798 PetscFunctionBegin; 1799 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1800 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 1801 PetscFunctionReturn(0); 1802 } 1803 1804 #undef __FUNCT__ 1805 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 1806 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 1807 { 1808 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1809 1810 PetscFunctionBegin; 1811 jac->schurfactorization = ftype; 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1817 /*@C 1818 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1819 1820 Collective on KSP 1821 1822 Input Parameter: 1823 . pc - the preconditioner context 1824 1825 Output Parameters: 1826 + A00 - the (0,0) block 1827 . A01 - the (0,1) block 1828 . A10 - the (1,0) block 1829 - A11 - the (1,1) block 1830 1831 Level: advanced 1832 1833 .seealso: PCFIELDSPLIT 1834 @*/ 1835 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1836 { 1837 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1838 1839 PetscFunctionBegin; 1840 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1841 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1842 if (A00) *A00 = jac->pmat[0]; 1843 if (A01) *A01 = jac->B; 1844 if (A10) *A10 = jac->C; 1845 if (A11) *A11 = jac->pmat[1]; 1846 PetscFunctionReturn(0); 1847 } 1848 1849 #undef __FUNCT__ 1850 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1851 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1852 { 1853 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1854 PetscErrorCode ierr; 1855 1856 PetscFunctionBegin; 1857 jac->type = type; 1858 if (type == PC_COMPOSITE_SCHUR) { 1859 pc->ops->apply = PCApply_FieldSplit_Schur; 1860 pc->ops->view = PCView_FieldSplit_Schur; 1861 1862 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1863 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1864 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 1865 1866 } else { 1867 pc->ops->apply = PCApply_FieldSplit; 1868 pc->ops->view = PCView_FieldSplit; 1869 1870 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1871 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSchurPrecondition_C",0);CHKERRQ(ierr); 1872 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 1873 } 1874 PetscFunctionReturn(0); 1875 } 1876 1877 #undef __FUNCT__ 1878 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1879 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1880 { 1881 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1882 1883 PetscFunctionBegin; 1884 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1885 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); 1886 jac->bs = bs; 1887 PetscFunctionReturn(0); 1888 } 1889 1890 #undef __FUNCT__ 1891 #define __FUNCT__ "PCFieldSplitSetType" 1892 /*@ 1893 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1894 1895 Collective on PC 1896 1897 Input Parameter: 1898 . pc - the preconditioner context 1899 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1900 1901 Options Database Key: 1902 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1903 1904 Level: Intermediate 1905 1906 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1907 1908 .seealso: PCCompositeSetType() 1909 1910 @*/ 1911 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1912 { 1913 PetscErrorCode ierr; 1914 1915 PetscFunctionBegin; 1916 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1917 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 1921 #undef __FUNCT__ 1922 #define __FUNCT__ "PCFieldSplitGetType" 1923 /*@ 1924 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 1925 1926 Not collective 1927 1928 Input Parameter: 1929 . pc - the preconditioner context 1930 1931 Output Parameter: 1932 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1933 1934 Level: Intermediate 1935 1936 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1937 .seealso: PCCompositeSetType() 1938 @*/ 1939 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 1940 { 1941 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 1942 1943 PetscFunctionBegin; 1944 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1945 PetscValidIntPointer(type,2); 1946 *type = jac->type; 1947 PetscFunctionReturn(0); 1948 } 1949 1950 #undef __FUNCT__ 1951 #define __FUNCT__ "PCFieldSplitSetDMSplits" 1952 /*@ 1953 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1954 1955 Logically Collective 1956 1957 Input Parameters: 1958 + pc - the preconditioner context 1959 - flg - boolean indicating whether to use field splits defined by the DM 1960 1961 Options Database Key: 1962 . -pc_fieldsplit_dm_splits 1963 1964 Level: Intermediate 1965 1966 .keywords: PC, DM, composite preconditioner, additive, multiplicative 1967 1968 .seealso: PCFieldSplitGetDMSplits() 1969 1970 @*/ 1971 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 1972 { 1973 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1974 PetscBool isfs; 1975 PetscErrorCode ierr; 1976 1977 PetscFunctionBegin; 1978 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1979 PetscValidLogicalCollectiveBool(pc,flg,2); 1980 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1981 if (isfs) { 1982 jac->dm_splits = flg; 1983 } 1984 PetscFunctionReturn(0); 1985 } 1986 1987 1988 #undef __FUNCT__ 1989 #define __FUNCT__ "PCFieldSplitGetDMSplits" 1990 /*@ 1991 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 1992 1993 Logically Collective 1994 1995 Input Parameter: 1996 . pc - the preconditioner context 1997 1998 Output Parameter: 1999 . flg - boolean indicating whether to use field splits defined by the DM 2000 2001 Level: Intermediate 2002 2003 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2004 2005 .seealso: PCFieldSplitSetDMSplits() 2006 2007 @*/ 2008 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2009 { 2010 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2011 PetscBool isfs; 2012 PetscErrorCode ierr; 2013 2014 PetscFunctionBegin; 2015 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2016 PetscValidPointer(flg,2); 2017 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2018 if (isfs) { 2019 if(flg) *flg = jac->dm_splits; 2020 } 2021 PetscFunctionReturn(0); 2022 } 2023 2024 /* -------------------------------------------------------------------------------------*/ 2025 /*MC 2026 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2027 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2028 2029 To set options on the solvers for each block append -fieldsplit_ to all the PC 2030 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2031 2032 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2033 and set the options directly on the resulting KSP object 2034 2035 Level: intermediate 2036 2037 Options Database Keys: 2038 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2039 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2040 been supplied explicitly by -pc_fieldsplit_%d_fields 2041 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2042 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2043 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11> - default is a11 2044 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2045 2046 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2047 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2048 2049 Notes: 2050 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2051 to define a field by an arbitrary collection of entries. 2052 2053 If no fields are set the default is used. The fields are defined by entries strided by bs, 2054 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2055 if this is not called the block size defaults to the blocksize of the second matrix passed 2056 to KSPSetOperators()/PCSetOperators(). 2057 2058 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2059 $ ( A10 A11 ) 2060 $ the preconditioner using full factorization is 2061 $ ( I -A10 ksp(A00) ) ( inv(A00) 0 ) ( I 0 ) 2062 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2063 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 2064 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given 2065 in providing the SECOND split or 1 if not given). For PCFieldSplitGetKSP() when field number is 0, 2066 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2067 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 2068 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc. 2069 When option -fieldsplit_schur_precondition selfp is given, an approximation to S is assembled -- 2070 Sp = A11 - A10 inv(diag(A00)) A01, which has type AIJ and can be used with a variety of preconditioners 2071 (e.g., -fieldsplit_1_pc_type asm). 2072 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2073 diag gives 2074 $ ( inv(A00) 0 ) 2075 $ ( 0 -ksp(S) ) 2076 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 2077 $ ( A00 0 ) 2078 $ ( A10 S ) 2079 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2080 $ ( A00 A01 ) 2081 $ ( 0 S ) 2082 where again the inverses of A00 and S are applied using KSPs. 2083 2084 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2085 is used automatically for a second block. 2086 2087 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2088 Generally it should be used with the AIJ format. 2089 2090 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2091 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2092 inside a smoother resulting in "Distributive Smoothers". 2093 2094 Concepts: physics based preconditioners, block preconditioners 2095 2096 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2097 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 2098 M*/ 2099 2100 #undef __FUNCT__ 2101 #define __FUNCT__ "PCCreate_FieldSplit" 2102 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2103 { 2104 PetscErrorCode ierr; 2105 PC_FieldSplit *jac; 2106 2107 PetscFunctionBegin; 2108 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2109 2110 jac->bs = -1; 2111 jac->nsplits = 0; 2112 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2113 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2114 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2115 jac->dm_splits = PETSC_TRUE; 2116 2117 pc->data = (void*)jac; 2118 2119 pc->ops->apply = PCApply_FieldSplit; 2120 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2121 pc->ops->setup = PCSetUp_FieldSplit; 2122 pc->ops->reset = PCReset_FieldSplit; 2123 pc->ops->destroy = PCDestroy_FieldSplit; 2124 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2125 pc->ops->view = PCView_FieldSplit; 2126 pc->ops->applyrichardson = 0; 2127 2128 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2129 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2130 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2131 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2132 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2133 PetscFunctionReturn(0); 2134 } 2135 2136 2137 2138