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