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