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