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_SCHUR_PRE_USER 1881 PC_FIELDSPLIT_SCHUR_PRE_SELFP, and PC_FIELDSPLIT_SCHUR_PRE_FULL 1882 - userpre - matrix to use for preconditioning, or NULL 1883 1884 Options Database: 1885 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11. See notes for meaning of various arguments 1886 1887 Notes: 1888 $ If ptype is 1889 $ a11 then the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 1890 $ matrix associated with the Schur complement (i.e. A11), not he Schur complement matrix 1891 $ self the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 1892 $ The only preconditioner that currently works with this symbolic respresentation matrix object is the PCLSC 1893 $ preconditioner 1894 $ user then the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 1895 $ to this function). 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 using the additional option -fieldsplit_1_mat_schur_complement_ainv_type lump 1899 $ full then the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation computed internally by PFIELDSPLIT (this is expensive) 1900 $ useful mostly as a test that the Schur complement approach can work for your problem 1901 1902 When solving a saddle point problem, where the A11 block is identically zero, using a11 as the ptype only makes sense 1903 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1904 -fieldsplit_1_pc_type lsc which uses the least squares commutator to compute a preconditioner for the Schur complement. 1905 1906 Level: intermediate 1907 1908 .seealso: PCFieldSplitGetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, 1909 MatSchurComplementSetAinvType(), PCLSC 1910 1911 @*/ 1912 PetscErrorCode PCFieldSplitSetSchurPre(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1913 { 1914 PetscErrorCode ierr; 1915 1916 PetscFunctionBegin; 1917 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1918 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurPre_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1919 PetscFunctionReturn(0); 1920 } 1921 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) {return PCFieldSplitSetSchurPre(pc,ptype,pre);} /* Deprecated name */ 1922 1923 #undef __FUNCT__ 1924 #define __FUNCT__ "PCFieldSplitGetSchurPre" 1925 /*@ 1926 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 1927 preconditioned. See PCFieldSplitSetSchurPre() for details. 1928 1929 Logically Collective on PC 1930 1931 Input Parameters: 1932 . pc - the preconditioner context 1933 1934 Output Parameters: 1935 + ptype - which matrix to use for preconditioning the Schur complement: PC_FIELDSPLIT_SCHUR_PRE_A11, PC_FIELDSPLIT_SCHUR_PRE_SELF, PC_FIELDSPLIT_PRE_USER 1936 - userpre - matrix to use for preconditioning (with PC_FIELDSPLIT_PRE_USER), or NULL 1937 1938 Level: intermediate 1939 1940 .seealso: PCFieldSplitSetSchurPre(), PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1941 1942 @*/ 1943 PetscErrorCode PCFieldSplitGetSchurPre(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 1944 { 1945 PetscErrorCode ierr; 1946 1947 PetscFunctionBegin; 1948 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1949 ierr = PetscUseMethod(pc,"PCFieldSplitGetSchurPre_C",(PC,PCFieldSplitSchurPreType*,Mat*),(pc,ptype,pre));CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 #undef __FUNCT__ 1954 #define __FUNCT__ "PCFieldSplitSchurGetS" 1955 /*@ 1956 PCFieldSplitSchurGetS - extract the MatSchurComplement object used by this PC in case it needs to be configured separately 1957 1958 Not collective 1959 1960 Input Parameter: 1961 . pc - the preconditioner context 1962 1963 Output Parameter: 1964 . S - the Schur complement matrix 1965 1966 Notes: 1967 This matrix should not be destroyed using MatDestroy(); rather, use PCFieldSplitSchurRestoreS(). 1968 1969 Level: advanced 1970 1971 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurRestoreS() 1972 1973 @*/ 1974 PetscErrorCode PCFieldSplitSchurGetS(PC pc,Mat *S) 1975 { 1976 PetscErrorCode ierr; 1977 const char* t; 1978 PetscBool isfs; 1979 PC_FieldSplit *jac; 1980 1981 PetscFunctionBegin; 1982 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1983 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 1984 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 1985 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 1986 jac = (PC_FieldSplit*)pc->data; 1987 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 1988 if (S) *S = jac->schur; 1989 PetscFunctionReturn(0); 1990 } 1991 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "PCFieldSplitSchurRestoreS" 1994 /*@ 1995 PCFieldSplitSchurRestoreS - restores the MatSchurComplement object used by this PC 1996 1997 Not collective 1998 1999 Input Parameters: 2000 + pc - the preconditioner context 2001 . S - the Schur complement matrix 2002 2003 Level: advanced 2004 2005 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSchurPreType, PCFieldSplitSetSchurPre(), MatSchurComplement, PCFieldSplitSchurGetS() 2006 2007 @*/ 2008 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc,Mat *S) 2009 { 2010 PetscErrorCode ierr; 2011 const char* t; 2012 PetscBool isfs; 2013 PC_FieldSplit *jac; 2014 2015 PetscFunctionBegin; 2016 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2017 ierr = PetscObjectGetType((PetscObject)pc,&t);CHKERRQ(ierr); 2018 ierr = PetscStrcmp(t,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2019 if (!isfs) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PC of type PCFIELDSPLIT, got %s instead",t); 2020 jac = (PC_FieldSplit*)pc->data; 2021 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected PCFIELDSPLIT of type SCHUR, got %D instead",jac->type); 2022 if (!S || *S != jac->schur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatSchurComplement restored is not the same as gotten"); 2023 PetscFunctionReturn(0); 2024 } 2025 2026 2027 #undef __FUNCT__ 2028 #define __FUNCT__ "PCFieldSplitSetSchurPre_FieldSplit" 2029 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 2030 { 2031 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2032 PetscErrorCode ierr; 2033 2034 PetscFunctionBegin; 2035 jac->schurpre = ptype; 2036 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2037 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 2038 jac->schur_user = pre; 2039 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 2040 } 2041 PetscFunctionReturn(0); 2042 } 2043 2044 #undef __FUNCT__ 2045 #define __FUNCT__ "PCFieldSplitGetSchurPre_FieldSplit" 2046 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc,PCFieldSplitSchurPreType *ptype,Mat *pre) 2047 { 2048 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2049 2050 PetscFunctionBegin; 2051 *ptype = jac->schurpre; 2052 *pre = jac->schur_user; 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "PCFieldSplitSetSchurFactType" 2058 /*@ 2059 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain 2060 2061 Collective on PC 2062 2063 Input Parameters: 2064 + pc - the preconditioner context 2065 - ftype - which blocks of factorization to retain, PC_FIELDSPLIT_SCHUR_FACT_FULL is default 2066 2067 Options Database: 2068 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> default is full 2069 2070 2071 Level: intermediate 2072 2073 Notes: 2074 The FULL factorization is 2075 2076 $ (A B) = (1 0) (A 0) (1 Ainv*B) 2077 $ (C D) (C*Ainv 1) (0 S) (0 1 ) 2078 2079 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, 2080 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). 2081 2082 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 2083 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 2084 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, 2085 the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so KSPGMRES converges in at most 4 iterations. 2086 2087 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 2088 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). 2089 2090 References: 2091 + 1. - Murphy, Golub, and Wathen, A note on preconditioning indefinite linear systems, SIAM J. Sci. Comput., 21 (2000). 2092 - 2. - Ipsen, A note on preconditioning nonsymmetric matrices, SIAM J. Sci. Comput., 23 (2001). 2093 2094 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType 2095 @*/ 2096 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc,PCFieldSplitSchurFactType ftype) 2097 { 2098 PetscErrorCode ierr; 2099 2100 PetscFunctionBegin; 2101 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2102 ierr = PetscTryMethod(pc,"PCFieldSplitSetSchurFactType_C",(PC,PCFieldSplitSchurFactType),(pc,ftype));CHKERRQ(ierr); 2103 PetscFunctionReturn(0); 2104 } 2105 2106 #undef __FUNCT__ 2107 #define __FUNCT__ "PCFieldSplitSetSchurFactType_FieldSplit" 2108 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc,PCFieldSplitSchurFactType ftype) 2109 { 2110 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2111 2112 PetscFunctionBegin; 2113 jac->schurfactorization = ftype; 2114 PetscFunctionReturn(0); 2115 } 2116 2117 #undef __FUNCT__ 2118 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 2119 /*@C 2120 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2121 2122 Collective on KSP 2123 2124 Input Parameter: 2125 . pc - the preconditioner context 2126 2127 Output Parameters: 2128 + A00 - the (0,0) block 2129 . A01 - the (0,1) block 2130 . A10 - the (1,0) block 2131 - A11 - the (1,1) block 2132 2133 Level: advanced 2134 2135 .seealso: PCFIELDSPLIT 2136 @*/ 2137 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 2138 { 2139 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2140 2141 PetscFunctionBegin; 2142 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2143 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2144 if (A00) *A00 = jac->pmat[0]; 2145 if (A01) *A01 = jac->B; 2146 if (A10) *A10 = jac->C; 2147 if (A11) *A11 = jac->pmat[1]; 2148 PetscFunctionReturn(0); 2149 } 2150 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 2153 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 2154 { 2155 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2156 PetscErrorCode ierr; 2157 2158 PetscFunctionBegin; 2159 jac->type = type; 2160 if (type == PC_COMPOSITE_SCHUR) { 2161 pc->ops->apply = PCApply_FieldSplit_Schur; 2162 pc->ops->view = PCView_FieldSplit_Schur; 2163 2164 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 2165 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",PCFieldSplitSetSchurPre_FieldSplit);CHKERRQ(ierr); 2166 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",PCFieldSplitGetSchurPre_FieldSplit);CHKERRQ(ierr); 2167 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",PCFieldSplitSetSchurFactType_FieldSplit);CHKERRQ(ierr); 2168 2169 } else { 2170 pc->ops->apply = PCApply_FieldSplit; 2171 pc->ops->view = PCView_FieldSplit; 2172 2173 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2174 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurPre_C",0);CHKERRQ(ierr); 2175 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSchurPre_C",0);CHKERRQ(ierr); 2176 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetSchurFactType_C",0);CHKERRQ(ierr); 2177 } 2178 PetscFunctionReturn(0); 2179 } 2180 2181 #undef __FUNCT__ 2182 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 2183 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 2184 { 2185 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2186 2187 PetscFunctionBegin; 2188 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 2189 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); 2190 jac->bs = bs; 2191 PetscFunctionReturn(0); 2192 } 2193 2194 #undef __FUNCT__ 2195 #define __FUNCT__ "PCFieldSplitSetType" 2196 /*@ 2197 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 2198 2199 Collective on PC 2200 2201 Input Parameter: 2202 . pc - the preconditioner context 2203 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2204 2205 Options Database Key: 2206 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2207 2208 Level: Intermediate 2209 2210 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2211 2212 .seealso: PCCompositeSetType() 2213 2214 @*/ 2215 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 2216 { 2217 PetscErrorCode ierr; 2218 2219 PetscFunctionBegin; 2220 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2221 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 2222 PetscFunctionReturn(0); 2223 } 2224 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "PCFieldSplitGetType" 2227 /*@ 2228 PCFieldSplitGetType - Gets the type of fieldsplit preconditioner. 2229 2230 Not collective 2231 2232 Input Parameter: 2233 . pc - the preconditioner context 2234 2235 Output Parameter: 2236 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 2237 2238 Level: Intermediate 2239 2240 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 2241 .seealso: PCCompositeSetType() 2242 @*/ 2243 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 2244 { 2245 PC_FieldSplit *jac = (PC_FieldSplit*) pc->data; 2246 2247 PetscFunctionBegin; 2248 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2249 PetscValidIntPointer(type,2); 2250 *type = jac->type; 2251 PetscFunctionReturn(0); 2252 } 2253 2254 #undef __FUNCT__ 2255 #define __FUNCT__ "PCFieldSplitSetDMSplits" 2256 /*@ 2257 PCFieldSplitSetDMSplits - Flags whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2258 2259 Logically Collective 2260 2261 Input Parameters: 2262 + pc - the preconditioner context 2263 - flg - boolean indicating whether to use field splits defined by the DM 2264 2265 Options Database Key: 2266 . -pc_fieldsplit_dm_splits 2267 2268 Level: Intermediate 2269 2270 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2271 2272 .seealso: PCFieldSplitGetDMSplits() 2273 2274 @*/ 2275 PetscErrorCode PCFieldSplitSetDMSplits(PC pc,PetscBool flg) 2276 { 2277 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2278 PetscBool isfs; 2279 PetscErrorCode ierr; 2280 2281 PetscFunctionBegin; 2282 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2283 PetscValidLogicalCollectiveBool(pc,flg,2); 2284 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2285 if (isfs) { 2286 jac->dm_splits = flg; 2287 } 2288 PetscFunctionReturn(0); 2289 } 2290 2291 2292 #undef __FUNCT__ 2293 #define __FUNCT__ "PCFieldSplitGetDMSplits" 2294 /*@ 2295 PCFieldSplitGetDMSplits - Returns flag indicating whether DMCreateFieldDecomposition() should be used to define the splits, whenever possible. 2296 2297 Logically Collective 2298 2299 Input Parameter: 2300 . pc - the preconditioner context 2301 2302 Output Parameter: 2303 . flg - boolean indicating whether to use field splits defined by the DM 2304 2305 Level: Intermediate 2306 2307 .keywords: PC, DM, composite preconditioner, additive, multiplicative 2308 2309 .seealso: PCFieldSplitSetDMSplits() 2310 2311 @*/ 2312 PetscErrorCode PCFieldSplitGetDMSplits(PC pc,PetscBool* flg) 2313 { 2314 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 2315 PetscBool isfs; 2316 PetscErrorCode ierr; 2317 2318 PetscFunctionBegin; 2319 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2320 PetscValidPointer(flg,2); 2321 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&isfs);CHKERRQ(ierr); 2322 if (isfs) { 2323 if(flg) *flg = jac->dm_splits; 2324 } 2325 PetscFunctionReturn(0); 2326 } 2327 2328 /* -------------------------------------------------------------------------------------*/ 2329 /*MC 2330 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 2331 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 2332 2333 To set options on the solvers for each block append -fieldsplit_ to all the PC 2334 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 2335 2336 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 2337 and set the options directly on the resulting KSP object 2338 2339 Level: intermediate 2340 2341 Options Database Keys: 2342 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 2343 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 2344 been supplied explicitly by -pc_fieldsplit_%d_fields 2345 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 2346 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur> - type of relaxation or factorization splitting 2347 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is a11; see PCFieldSplitSetSchurPre() 2348 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 2349 2350 - Options prefix for inner solvers when using Schur complement preconditioner are -fieldsplit_0_ and -fieldsplit_1_ 2351 for all other solvers they are -fieldsplit_%d_ for the dth field, use -fieldsplit_ for all fields 2352 2353 Notes: 2354 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 2355 to define a field by an arbitrary collection of entries. 2356 2357 If no fields are set the default is used. The fields are defined by entries strided by bs, 2358 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 2359 if this is not called the block size defaults to the blocksize of the second matrix passed 2360 to KSPSetOperators()/PCSetOperators(). 2361 2362 $ For the Schur complement preconditioner if J = ( A00 A01 ) 2363 $ ( A10 A11 ) 2364 $ the preconditioner using full factorization is 2365 $ ( I -ksp(A00) A01 ) ( inv(A00) 0 ) ( I 0 ) 2366 $ ( 0 I ) ( 0 ksp(S) ) ( -A10 ksp(A00) I ) 2367 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. S is the Schur complement 2368 $ S = A11 - A10 ksp(A00) A01 2369 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 2370 in providing the SECOND split or 1 if not give). For PCFieldSplitGetKSP() when field number is 0, 2371 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 2372 A11 is used to construct a preconditioner for S, use PCFieldSplitSetSchurPre() for all the possible ways to construct the preconditioner for S. 2373 2374 The factorization type is set using -pc_fieldsplit_schur_fact_type <diag, lower, upper, full>. The full is shown above, 2375 diag gives 2376 $ ( inv(A00) 0 ) 2377 $ ( 0 -ksp(S) ) 2378 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 2379 $ ( A00 0 ) 2380 $ ( A10 S ) 2381 where the inverses of A00 and S are applied using KSPs. The upper factorization is the inverse of 2382 $ ( A00 A01 ) 2383 $ ( 0 S ) 2384 where again the inverses of A00 and S are applied using KSPs. 2385 2386 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 2387 is used automatically for a second block. 2388 2389 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 2390 Generally it should be used with the AIJ format. 2391 2392 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 2393 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 2394 inside a smoother resulting in "Distributive Smoothers". 2395 2396 Concepts: physics based preconditioners, block preconditioners 2397 2398 There is a nice discussion of block preconditioners in 2399 2400 [El08] A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations 2401 Howard Elman, V.E. Howle, John Shadid, Robert Shuttleworth, Ray Tuminaro, Journal of Computational Physics 227 (2008) 1790--1808 2402 http://chess.cs.umd.edu/~elman/papers/tax.pdf 2403 2404 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 2405 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 2406 of backsolving for the saturations in the last step it solves a full coupled (ILU) system for updates to all the variables. 2407 2408 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 2409 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 2410 MatSchurComplementSetAinvType() 2411 M*/ 2412 2413 #undef __FUNCT__ 2414 #define __FUNCT__ "PCCreate_FieldSplit" 2415 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 2416 { 2417 PetscErrorCode ierr; 2418 PC_FieldSplit *jac; 2419 2420 PetscFunctionBegin; 2421 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2422 2423 jac->bs = -1; 2424 jac->nsplits = 0; 2425 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 2426 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 2427 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 2428 jac->dm_splits = PETSC_TRUE; 2429 2430 pc->data = (void*)jac; 2431 2432 pc->ops->apply = PCApply_FieldSplit; 2433 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 2434 pc->ops->setup = PCSetUp_FieldSplit; 2435 pc->ops->reset = PCReset_FieldSplit; 2436 pc->ops->destroy = PCDestroy_FieldSplit; 2437 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 2438 pc->ops->view = PCView_FieldSplit; 2439 pc->ops->applyrichardson = 0; 2440 2441 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 2442 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetFields_C",PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 2443 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetIS_C",PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 2444 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetType_C",PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 2445 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitSetBlockSize_C",PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 2446 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCFieldSplitRestrictIS_C",PCFieldSplitRestrictIS_FieldSplit);CHKERRQ(ierr); 2447 PetscFunctionReturn(0); 2448 } 2449 2450 2451 2452