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