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