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