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