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