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