1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* This is needed to provide the appropriate PETSC_EXTERN for KSP_Solve_FS ....*/ 3 #include <petscdm.h> 4 5 const char *const PCFieldSplitSchurPreTypes[] = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL}; 6 const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL}; 7 8 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; 9 10 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 11 struct _PC_FieldSplitLink { 12 KSP ksp; 13 Vec x, y, z; 14 char *splitname; 15 PetscInt nfields; 16 PetscInt *fields, *fields_col; 17 VecScatter sctx; 18 IS is, is_col; 19 PC_FieldSplitLink next, previous; 20 PetscLogEvent event; 21 22 /* Used only when setting coordinates with PCSetCoordinates */ 23 PetscInt dim; 24 PetscInt ndofs; 25 PetscReal *coords; 26 }; 27 28 typedef struct { 29 PCCompositeType type; 30 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 31 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 32 PetscInt bs; /* Block size for IS and Mat structures */ 33 PetscInt nsplits; /* Number of field divisions defined */ 34 Vec *x, *y, w1, w2; 35 Mat *mat; /* The diagonal block for each split */ 36 Mat *pmat; /* The preconditioning diagonal block for each split */ 37 Mat *Afield; /* The rows of the matrix associated with each split */ 38 PetscBool issetup; 39 40 /* Only used when Schur complement preconditioning is used */ 41 Mat B; /* The (0,1) block */ 42 Mat C; /* The (1,0) block */ 43 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */ 44 Mat schurp; /* Assembled approximation to S built by MatSchurComplement to be used as a preconditioning matrix when solving with S */ 45 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 46 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 47 PCFieldSplitSchurFactType schurfactorization; 48 KSP kspschur; /* The solver for S */ 49 KSP kspupper; /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */ 50 PetscScalar schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */ 51 52 /* Only used when Golub-Kahan bidiagonalization preconditioning is used */ 53 Mat H; /* The modified matrix H = A00 + nu*A01*A01' */ 54 PetscReal gkbtol; /* Stopping tolerance for lower bound estimate */ 55 PetscInt gkbdelay; /* The delay window for the stopping criterion */ 56 PetscReal gkbnu; /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */ 57 PetscInt gkbmaxit; /* Maximum number of iterations for outer loop */ 58 PetscBool gkbmonitor; /* Monitor for gkb iterations and the lower bound error */ 59 PetscViewer gkbviewer; /* Viewer context for gkbmonitor */ 60 Vec u, v, d, Hu; /* Work vectors for the GKB algorithm */ 61 PetscScalar *vecz; /* Contains intermediate values, eg for lower bound */ 62 63 PC_FieldSplitLink head; 64 PetscBool isrestrict; /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */ 65 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 66 PetscBool dm_splits; /* Whether to use DMCreateFieldDecomposition() whenever possible */ 67 PetscBool diag_use_amat; /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 68 PetscBool offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */ 69 PetscBool detect; /* Whether to form 2-way split by finding zero diagonal entries */ 70 PetscBool coordinates_set; /* Whether PCSetCoordinates has been called */ 71 } PC_FieldSplit; 72 73 /* 74 Note: 75 there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 76 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 77 PC you could change this. 78 */ 79 80 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 81 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 82 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 83 { 84 switch (jac->schurpre) { 85 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 86 return jac->schur; 87 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 88 return jac->schurp; 89 case PC_FIELDSPLIT_SCHUR_PRE_A11: 90 return jac->pmat[1]; 91 case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */ 92 case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */ 93 default: 94 return jac->schur_user ? jac->schur_user : jac->pmat[1]; 95 } 96 } 97 98 #include <petscdraw.h> 99 static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer) 100 { 101 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 102 PetscBool iascii, isdraw; 103 PetscInt i, j; 104 PC_FieldSplitLink ilink = jac->head; 105 106 PetscFunctionBegin; 107 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 108 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 109 if (iascii) { 110 if (jac->bs > 0) { 111 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs)); 112 } else { 113 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits)); 114 } 115 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 116 if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n")); 117 if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n")); 118 PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for each split is in the following KSP objects:\n")); 119 for (i = 0; i < jac->nsplits; i++) { 120 if (ilink->fields) { 121 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i)); 122 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 123 for (j = 0; j < ilink->nfields; j++) { 124 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 125 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 126 } 127 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 128 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 129 } else { 130 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i)); 131 } 132 PetscCall(KSPView(ilink->ksp, viewer)); 133 ilink = ilink->next; 134 } 135 } 136 137 if (isdraw) { 138 PetscDraw draw; 139 PetscReal x, y, w, wd; 140 141 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 142 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 143 w = 2 * PetscMin(1.0 - x, x); 144 wd = w / (jac->nsplits + 1); 145 x = x - wd * (jac->nsplits - 1) / 2.0; 146 for (i = 0; i < jac->nsplits; i++) { 147 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 148 PetscCall(KSPView(ilink->ksp, viewer)); 149 PetscCall(PetscDrawPopCurrentPoint(draw)); 150 x += wd; 151 ilink = ilink->next; 152 } 153 } 154 PetscFunctionReturn(PETSC_SUCCESS); 155 } 156 157 static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer) 158 { 159 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 160 PetscBool iascii, isdraw; 161 PetscInt i, j; 162 PC_FieldSplitLink ilink = jac->head; 163 MatSchurComplementAinvType atype; 164 165 PetscFunctionBegin; 166 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 167 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 168 if (iascii) { 169 if (jac->bs > 0) { 170 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization])); 171 } else { 172 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization])); 173 } 174 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 175 switch (jac->schurpre) { 176 case PC_FIELDSPLIT_SCHUR_PRE_SELF: 177 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from S itself\n")); 178 break; 179 case PC_FIELDSPLIT_SCHUR_PRE_SELFP: 180 if (jac->schur) { 181 PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype)); 182 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's ")))); 183 } 184 break; 185 case PC_FIELDSPLIT_SCHUR_PRE_A11: 186 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n")); 187 break; 188 case PC_FIELDSPLIT_SCHUR_PRE_FULL: 189 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from the exact Schur complement\n")); 190 break; 191 case PC_FIELDSPLIT_SCHUR_PRE_USER: 192 if (jac->schur_user) { 193 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from user provided matrix\n")); 194 } else { 195 PetscCall(PetscViewerASCIIPrintf(viewer, " Preconditioner for the Schur complement formed from A11\n")); 196 } 197 break; 198 default: 199 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre); 200 } 201 PetscCall(PetscViewerASCIIPrintf(viewer, " Split info:\n")); 202 PetscCall(PetscViewerASCIIPushTab(viewer)); 203 for (i = 0; i < jac->nsplits; i++) { 204 if (ilink->fields) { 205 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i)); 206 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 207 for (j = 0; j < ilink->nfields; j++) { 208 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 209 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 210 } 211 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 212 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 213 } else { 214 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i)); 215 } 216 ilink = ilink->next; 217 } 218 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n")); 219 PetscCall(PetscViewerASCIIPushTab(viewer)); 220 if (jac->head) { 221 PetscCall(KSPView(jac->head->ksp, viewer)); 222 } else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 223 PetscCall(PetscViewerASCIIPopTab(viewer)); 224 if (jac->head && jac->kspupper != jac->head->ksp) { 225 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n")); 226 PetscCall(PetscViewerASCIIPushTab(viewer)); 227 if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer)); 228 else PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 229 PetscCall(PetscViewerASCIIPopTab(viewer)); 230 } 231 PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n")); 232 PetscCall(PetscViewerASCIIPushTab(viewer)); 233 if (jac->kspschur) { 234 PetscCall(KSPView(jac->kspschur, viewer)); 235 } else { 236 PetscCall(PetscViewerASCIIPrintf(viewer, " not yet available\n")); 237 } 238 PetscCall(PetscViewerASCIIPopTab(viewer)); 239 PetscCall(PetscViewerASCIIPopTab(viewer)); 240 } else if (isdraw && jac->head) { 241 PetscDraw draw; 242 PetscReal x, y, w, wd, h; 243 PetscInt cnt = 2; 244 char str[32]; 245 246 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 247 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 248 if (jac->kspupper != jac->head->ksp) cnt++; 249 w = 2 * PetscMin(1.0 - x, x); 250 wd = w / (cnt + 1); 251 252 PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization])); 253 PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h)); 254 y -= h; 255 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) { 256 PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11])); 257 } else { 258 PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre])); 259 } 260 PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h)); 261 y -= h; 262 x = x - wd * (cnt - 1) / 2.0; 263 264 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 265 PetscCall(KSPView(jac->head->ksp, viewer)); 266 PetscCall(PetscDrawPopCurrentPoint(draw)); 267 if (jac->kspupper != jac->head->ksp) { 268 x += wd; 269 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 270 PetscCall(KSPView(jac->kspupper, viewer)); 271 PetscCall(PetscDrawPopCurrentPoint(draw)); 272 } 273 x += wd; 274 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 275 PetscCall(KSPView(jac->kspschur, viewer)); 276 PetscCall(PetscDrawPopCurrentPoint(draw)); 277 } 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer) 282 { 283 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 284 PetscBool iascii, isdraw; 285 PetscInt i, j; 286 PC_FieldSplitLink ilink = jac->head; 287 288 PetscFunctionBegin; 289 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 290 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 291 if (iascii) { 292 if (jac->bs > 0) { 293 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs)); 294 } else { 295 PetscCall(PetscViewerASCIIPrintf(viewer, " FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits)); 296 } 297 if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for blocks\n")); 298 if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for diagonal blocks\n")); 299 if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat (not Pmat) as operator for off-diagonal blocks\n")); 300 301 PetscCall(PetscViewerASCIIPrintf(viewer, " Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit)); 302 PetscCall(PetscViewerASCIIPrintf(viewer, " Solver info for H = A00 + nu*A01*A01' matrix:\n")); 303 PetscCall(PetscViewerASCIIPushTab(viewer)); 304 305 if (ilink->fields) { 306 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields ")); 307 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 308 for (j = 0; j < ilink->nfields; j++) { 309 if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ",")); 310 PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j])); 311 } 312 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 313 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 314 } else { 315 PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n")); 316 } 317 PetscCall(KSPView(ilink->ksp, viewer)); 318 319 PetscCall(PetscViewerASCIIPopTab(viewer)); 320 } 321 322 if (isdraw) { 323 PetscDraw draw; 324 PetscReal x, y, w, wd; 325 326 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 327 PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y)); 328 w = 2 * PetscMin(1.0 - x, x); 329 wd = w / (jac->nsplits + 1); 330 x = x - wd * (jac->nsplits - 1) / 2.0; 331 for (i = 0; i < jac->nsplits; i++) { 332 PetscCall(PetscDrawPushCurrentPoint(draw, x, y)); 333 PetscCall(KSPView(ilink->ksp, viewer)); 334 PetscCall(PetscDrawPopCurrentPoint(draw)); 335 x += wd; 336 ilink = ilink->next; 337 } 338 } 339 PetscFunctionReturn(PETSC_SUCCESS); 340 } 341 342 /* Precondition: jac->bs is set to a meaningful value */ 343 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 344 { 345 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 346 PetscInt i, nfields, *ifields, nfields_col, *ifields_col; 347 PetscBool flg, flg_col; 348 char optionname[128], splitname[8], optionname_col[128]; 349 350 PetscFunctionBegin; 351 PetscCall(PetscMalloc1(jac->bs, &ifields)); 352 PetscCall(PetscMalloc1(jac->bs, &ifields_col)); 353 for (i = 0, flg = PETSC_TRUE;; i++) { 354 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 355 PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i)); 356 PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i)); 357 nfields = jac->bs; 358 nfields_col = jac->bs; 359 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 360 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col)); 361 if (!flg) break; 362 else if (flg && !flg_col) { 363 PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields"); 364 PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields)); 365 } else { 366 PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields"); 367 PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match"); 368 PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col)); 369 } 370 } 371 if (i > 0) { 372 /* Makes command-line setting of splits take precedence over setting them in code. 373 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 374 create new splits, which would probably not be what the user wanted. */ 375 jac->splitdefined = PETSC_TRUE; 376 } 377 PetscCall(PetscFree(ifields)); 378 PetscCall(PetscFree(ifields_col)); 379 PetscFunctionReturn(PETSC_SUCCESS); 380 } 381 382 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 383 { 384 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 385 PC_FieldSplitLink ilink = jac->head; 386 PetscBool fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE; 387 PetscInt i; 388 389 PetscFunctionBegin; 390 /* 391 Kinda messy, but at least this now uses DMCreateFieldDecomposition(). 392 Should probably be rewritten. 393 */ 394 if (!ilink) { 395 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL)); 396 if (pc->dm && jac->dm_splits && !jac->detect && !coupling) { 397 PetscInt numFields, f, i, j; 398 char **fieldNames; 399 IS *fields; 400 DM *dms; 401 DM subdm[128]; 402 PetscBool flg; 403 404 PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms)); 405 /* Allow the user to prescribe the splits */ 406 for (i = 0, flg = PETSC_TRUE;; i++) { 407 PetscInt ifields[128]; 408 IS compField; 409 char optionname[128], splitname[8]; 410 PetscInt nfields = numFields; 411 412 PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i)); 413 PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg)); 414 if (!flg) break; 415 PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields); 416 PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i])); 417 if (nfields == 1) { 418 PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField)); 419 } else { 420 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 421 PetscCall(PCFieldSplitSetIS(pc, splitname, compField)); 422 } 423 PetscCall(ISDestroy(&compField)); 424 for (j = 0; j < nfields; ++j) { 425 f = ifields[j]; 426 PetscCall(PetscFree(fieldNames[f])); 427 PetscCall(ISDestroy(&fields[f])); 428 } 429 } 430 if (i == 0) { 431 for (f = 0; f < numFields; ++f) { 432 PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f])); 433 PetscCall(PetscFree(fieldNames[f])); 434 PetscCall(ISDestroy(&fields[f])); 435 } 436 } else { 437 for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j)); 438 PetscCall(PetscFree(dms)); 439 PetscCall(PetscMalloc1(i, &dms)); 440 for (j = 0; j < i; ++j) dms[j] = subdm[j]; 441 } 442 PetscCall(PetscFree(fieldNames)); 443 PetscCall(PetscFree(fields)); 444 if (dms) { 445 PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n")); 446 for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) { 447 const char *prefix; 448 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix)); 449 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix)); 450 PetscCall(KSPSetDM(ilink->ksp, dms[i])); 451 PetscCall(KSPSetDMActive(ilink->ksp, PETSC_FALSE)); 452 { 453 PetscErrorCode (*func)(KSP, Mat, Mat, void *); 454 void *ctx; 455 456 PetscCall(DMKSPGetComputeOperators(pc->dm, &func, &ctx)); 457 PetscCall(DMKSPSetComputeOperators(dms[i], func, ctx)); 458 } 459 PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0)); 460 PetscCall(DMDestroy(&dms[i])); 461 } 462 PetscCall(PetscFree(dms)); 463 } 464 } else { 465 if (jac->bs <= 0) { 466 if (pc->pmat) { 467 PetscCall(MatGetBlockSize(pc->pmat, &jac->bs)); 468 } else jac->bs = 1; 469 } 470 471 if (jac->detect) { 472 IS zerodiags, rest; 473 PetscInt nmin, nmax; 474 475 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 476 if (jac->diag_use_amat) { 477 PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags)); 478 } else { 479 PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags)); 480 } 481 PetscCall(ISComplement(zerodiags, nmin, nmax, &rest)); 482 PetscCall(PCFieldSplitSetIS(pc, "0", rest)); 483 PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags)); 484 PetscCall(ISDestroy(&zerodiags)); 485 PetscCall(ISDestroy(&rest)); 486 } else if (coupling) { 487 IS coupling, rest; 488 PetscInt nmin, nmax; 489 490 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 491 if (jac->offdiag_use_amat) { 492 PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling)); 493 } else { 494 PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling)); 495 } 496 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest)); 497 PetscCall(ISSetIdentity(rest)); 498 PetscCall(PCFieldSplitSetIS(pc, "0", rest)); 499 PetscCall(PCFieldSplitSetIS(pc, "1", coupling)); 500 PetscCall(ISDestroy(&coupling)); 501 PetscCall(ISDestroy(&rest)); 502 } else { 503 PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL)); 504 if (!fieldsplit_default) { 505 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 506 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 507 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 508 if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n")); 509 } 510 if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) { 511 Mat M = pc->pmat; 512 PetscBool isnest; 513 514 PetscCall(PetscInfo(pc, "Using default splitting of fields\n")); 515 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest)); 516 if (!isnest) { 517 M = pc->mat; 518 PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest)); 519 } 520 if (isnest) { 521 IS *fields; 522 PetscInt nf; 523 524 PetscCall(MatNestGetSize(M, &nf, NULL)); 525 PetscCall(PetscMalloc1(nf, &fields)); 526 PetscCall(MatNestGetISs(M, fields, NULL)); 527 for (i = 0; i < nf; i++) PetscCall(PCFieldSplitSetIS(pc, NULL, fields[i])); 528 PetscCall(PetscFree(fields)); 529 } else { 530 for (i = 0; i < jac->bs; i++) { 531 char splitname[8]; 532 PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i)); 533 PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i)); 534 } 535 jac->defaultsplit = PETSC_TRUE; 536 } 537 } 538 } 539 } 540 } else if (jac->nsplits == 1) { 541 IS is2; 542 PetscInt nmin, nmax; 543 544 PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()"); 545 PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax)); 546 PetscCall(ISComplement(ilink->is, nmin, nmax, &is2)); 547 PetscCall(PCFieldSplitSetIS(pc, "1", is2)); 548 PetscCall(ISDestroy(&is2)); 549 } 550 551 PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits); 552 PetscFunctionReturn(PETSC_SUCCESS); 553 } 554 555 static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu) 556 { 557 Mat BT, T; 558 PetscReal nrmT, nrmB; 559 560 PetscFunctionBegin; 561 PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */ 562 PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN)); 563 PetscCall(MatNorm(T, NORM_1, &nrmT)); 564 PetscCall(MatNorm(B, NORM_1, &nrmB)); 565 PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/hermitian, GKB is not applicable."); 566 567 /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */ 568 /* setting N := 1/nu*I in [Ar13]. */ 569 PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT)); 570 PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, H)); /* H = A01*A01' */ 571 PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN)); /* H = A00 + nu*A01*A01' */ 572 573 PetscCall(MatDestroy(&BT)); 574 PetscCall(MatDestroy(&T)); 575 PetscFunctionReturn(PETSC_SUCCESS); 576 } 577 578 PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg); 579 580 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 581 { 582 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 583 PC_FieldSplitLink ilink; 584 PetscInt i, nsplit; 585 PetscBool sorted, sorted_col; 586 587 PetscFunctionBegin; 588 pc->failedreason = PC_NOERROR; 589 PetscCall(PCFieldSplitSetDefaults(pc)); 590 nsplit = jac->nsplits; 591 ilink = jac->head; 592 593 /* get the matrices for each split */ 594 if (!jac->issetup) { 595 PetscInt rstart, rend, nslots, bs; 596 597 jac->issetup = PETSC_TRUE; 598 599 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 600 if (jac->defaultsplit || !ilink->is) { 601 if (jac->bs <= 0) jac->bs = nsplit; 602 } 603 604 /* MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */ 605 PetscCall(MatGetBlockSize(pc->pmat, &bs)); 606 if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) { 607 PetscBool blk; 608 609 PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL)); 610 PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes"); 611 } 612 613 bs = jac->bs; 614 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 615 nslots = (rend - rstart) / bs; 616 for (i = 0; i < nsplit; i++) { 617 if (jac->defaultsplit) { 618 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is)); 619 PetscCall(ISDuplicate(ilink->is, &ilink->is_col)); 620 } else if (!ilink->is) { 621 if (ilink->nfields > 1) { 622 PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col; 623 PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii)); 624 PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj)); 625 for (j = 0; j < nslots; j++) { 626 for (k = 0; k < nfields; k++) { 627 ii[nfields * j + k] = rstart + bs * j + fields[k]; 628 jj[nfields * j + k] = rstart + bs * j + fields_col[k]; 629 } 630 } 631 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is)); 632 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col)); 633 PetscCall(ISSetBlockSize(ilink->is, nfields)); 634 PetscCall(ISSetBlockSize(ilink->is_col, nfields)); 635 } else { 636 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is)); 637 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col)); 638 } 639 } 640 PetscCall(ISSorted(ilink->is, &sorted)); 641 if (ilink->is_col) PetscCall(ISSorted(ilink->is_col, &sorted_col)); 642 PetscCheck(sorted && sorted_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split"); 643 ilink = ilink->next; 644 } 645 } 646 647 ilink = jac->head; 648 if (!jac->pmat) { 649 Vec xtmp; 650 651 PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL)); 652 PetscCall(PetscMalloc1(nsplit, &jac->pmat)); 653 PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y)); 654 for (i = 0; i < nsplit; i++) { 655 MatNullSpace sp; 656 657 /* Check for preconditioning matrix attached to IS */ 658 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i])); 659 if (jac->pmat[i]) { 660 PetscCall(PetscObjectReference((PetscObject)jac->pmat[i])); 661 if (jac->type == PC_COMPOSITE_SCHUR) { 662 jac->schur_user = jac->pmat[i]; 663 664 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 665 } 666 } else { 667 const char *prefix; 668 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i])); 669 PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix)); 670 if (!prefix) { 671 PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix)); 672 PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix)); 673 } 674 PetscCall(MatSetFromOptions(jac->pmat[i])); 675 PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view")); 676 } 677 /* create work vectors for each split */ 678 PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i])); 679 ilink->x = jac->x[i]; 680 ilink->y = jac->y[i]; 681 ilink->z = NULL; 682 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 683 PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx)); 684 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp)); 685 if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp)); 686 ilink = ilink->next; 687 } 688 PetscCall(VecDestroy(&xtmp)); 689 } else { 690 MatReuse scall; 691 MatNullSpace *nullsp = NULL; 692 693 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 694 PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp)); 695 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i])); 696 scall = MAT_INITIAL_MATRIX; 697 } else scall = MAT_REUSE_MATRIX; 698 699 for (i = 0; i < nsplit; i++) { 700 Mat pmat; 701 702 /* Check for preconditioning matrix attached to IS */ 703 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat)); 704 if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i])); 705 ilink = ilink->next; 706 } 707 if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp)); 708 } 709 if (jac->diag_use_amat) { 710 ilink = jac->head; 711 if (!jac->mat) { 712 PetscCall(PetscMalloc1(nsplit, &jac->mat)); 713 for (i = 0; i < nsplit; i++) { 714 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i])); 715 ilink = ilink->next; 716 } 717 } else { 718 MatReuse scall; 719 MatNullSpace *nullsp = NULL; 720 721 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 722 PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp)); 723 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i])); 724 scall = MAT_INITIAL_MATRIX; 725 } else scall = MAT_REUSE_MATRIX; 726 727 for (i = 0; i < nsplit; i++) { 728 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i])); 729 ilink = ilink->next; 730 } 731 if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp)); 732 } 733 } else { 734 jac->mat = jac->pmat; 735 } 736 737 /* Check for null space attached to IS */ 738 ilink = jac->head; 739 for (i = 0; i < nsplit; i++) { 740 MatNullSpace sp; 741 742 PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp)); 743 if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp)); 744 ilink = ilink->next; 745 } 746 747 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) { 748 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 749 /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */ 750 ilink = jac->head; 751 if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) { 752 /* 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 */ 753 if (!jac->Afield) { 754 PetscCall(PetscCalloc1(nsplit, &jac->Afield)); 755 if (jac->offdiag_use_amat) { 756 PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1])); 757 } else { 758 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1])); 759 } 760 } else { 761 MatReuse scall; 762 763 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 764 PetscCall(MatDestroy(&jac->Afield[1])); 765 scall = MAT_INITIAL_MATRIX; 766 } else scall = MAT_REUSE_MATRIX; 767 768 if (jac->offdiag_use_amat) { 769 PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1])); 770 } else { 771 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1])); 772 } 773 } 774 } else { 775 if (!jac->Afield) { 776 PetscCall(PetscMalloc1(nsplit, &jac->Afield)); 777 for (i = 0; i < nsplit; i++) { 778 if (jac->offdiag_use_amat) { 779 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i])); 780 } else { 781 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i])); 782 } 783 ilink = ilink->next; 784 } 785 } else { 786 MatReuse scall; 787 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 788 for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i])); 789 scall = MAT_INITIAL_MATRIX; 790 } else scall = MAT_REUSE_MATRIX; 791 792 for (i = 0; i < nsplit; i++) { 793 if (jac->offdiag_use_amat) { 794 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i])); 795 } else { 796 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i])); 797 } 798 ilink = ilink->next; 799 } 800 } 801 } 802 } 803 804 if (jac->type == PC_COMPOSITE_SCHUR) { 805 IS ccis; 806 PetscBool isset, isspd; 807 PetscInt rstart, rend; 808 char lscname[256]; 809 PetscObject LSC_L; 810 PetscBool set, flg; 811 812 PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields"); 813 814 /* If pc->mat is SPD, don't scale by -1 the Schur complement */ 815 if (jac->schurscale == (PetscScalar)-1.0) { 816 PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd)); 817 jac->schurscale = (isset && isspd) ? 1.0 : -1.0; 818 } 819 820 /* When extracting off-diagonal submatrices, we take complements from this range */ 821 PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend)); 822 PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 823 824 if (jac->schur) { 825 KSP kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper; 826 MatReuse scall; 827 828 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 829 scall = MAT_INITIAL_MATRIX; 830 PetscCall(MatDestroy(&jac->B)); 831 PetscCall(MatDestroy(&jac->C)); 832 } else scall = MAT_REUSE_MATRIX; 833 834 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 835 ilink = jac->head; 836 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 837 if (jac->offdiag_use_amat) { 838 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->B)); 839 } else { 840 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->B)); 841 } 842 PetscCall(ISDestroy(&ccis)); 843 if (!flg) { 844 ilink = ilink->next; 845 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 846 if (jac->offdiag_use_amat) { 847 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->C)); 848 } else { 849 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->C)); 850 } 851 PetscCall(ISDestroy(&ccis)); 852 } else { 853 PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg)); 854 if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C)); 855 else PetscCall(MatCreateTranspose(jac->B, &jac->C)); 856 } 857 PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1])); 858 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) { 859 PetscCall(MatDestroy(&jac->schurp)); 860 PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp)); 861 } 862 if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0])); 863 if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0])); 864 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac))); 865 } else { 866 const char *Dprefix; 867 char schurprefix[256], schurmatprefix[256]; 868 char schurtestoption[256]; 869 MatNullSpace sp; 870 KSP kspt; 871 872 /* extract the A01 and A10 matrices */ 873 ilink = jac->head; 874 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 875 if (jac->offdiag_use_amat) { 876 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); 877 } else { 878 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); 879 } 880 PetscCall(ISDestroy(&ccis)); 881 ilink = ilink->next; 882 if (!flg) { 883 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 884 if (jac->offdiag_use_amat) { 885 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C)); 886 } else { 887 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C)); 888 } 889 PetscCall(ISDestroy(&ccis)); 890 } else { 891 PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg)); 892 if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C)); 893 else PetscCall(MatCreateTranspose(jac->B, &jac->C)); 894 } 895 /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */ 896 PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur)); 897 PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT)); 898 PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1])); 899 PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 900 PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix)); 901 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt)); 902 PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix)); 903 904 /* Note: this is not true in general */ 905 PetscCall(MatGetNullSpace(jac->mat[1], &sp)); 906 if (sp) PetscCall(MatSetNullSpace(jac->schur, sp)); 907 908 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname)); 909 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg)); 910 if (flg) { 911 DM dmInner; 912 KSP kspInner; 913 PC pcInner; 914 915 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 916 PetscCall(KSPReset(kspInner)); 917 PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0])); 918 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 919 /* Indent this deeper to emphasize the "inner" nature of this solver. */ 920 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2)); 921 PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2)); 922 PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix)); 923 924 /* Set DM for new solver */ 925 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 926 PetscCall(KSPSetDM(kspInner, dmInner)); 927 PetscCall(KSPSetDMActive(kspInner, PETSC_FALSE)); 928 929 /* Defaults to PCKSP as preconditioner */ 930 PetscCall(KSPGetPC(kspInner, &pcInner)); 931 PetscCall(PCSetType(pcInner, PCKSP)); 932 PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp)); 933 } else { 934 /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or 935 * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact, 936 * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for 937 * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make 938 * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used 939 * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */ 940 PetscCall(KSPSetType(jac->head->ksp, KSPGMRES)); 941 PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp)); 942 } 943 PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0])); 944 PetscCall(KSPSetFromOptions(jac->head->ksp)); 945 PetscCall(MatSetFromOptions(jac->schur)); 946 947 PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg)); 948 if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */ 949 KSP kspInner; 950 PC pcInner; 951 952 PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner)); 953 PetscCall(KSPGetPC(kspInner, &pcInner)); 954 PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg)); 955 if (flg) { 956 KSP ksp; 957 958 PetscCall(PCKSPGetKSP(pcInner, &ksp)); 959 if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE)); 960 } 961 } 962 PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname)); 963 PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg)); 964 if (flg) { 965 DM dmInner; 966 967 PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 968 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper)); 969 PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel)); 970 PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure)); 971 PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix)); 972 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1)); 973 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1)); 974 PetscCall(KSPGetDM(jac->head->ksp, &dmInner)); 975 PetscCall(KSPSetDM(jac->kspupper, dmInner)); 976 PetscCall(KSPSetDMActive(jac->kspupper, PETSC_FALSE)); 977 PetscCall(KSPSetFromOptions(jac->kspupper)); 978 PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0])); 979 PetscCall(VecDuplicate(jac->head->x, &jac->head->z)); 980 } else { 981 jac->kspupper = jac->head->ksp; 982 PetscCall(PetscObjectReference((PetscObject)jac->head->ksp)); 983 } 984 985 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp)); 986 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur)); 987 PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel)); 988 PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure)); 989 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1)); 990 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 991 PC pcschur; 992 PetscCall(KSPGetPC(jac->kspschur, &pcschur)); 993 PetscCall(PCSetType(pcschur, PCNONE)); 994 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 995 } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) { 996 PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user)); 997 } 998 PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac))); 999 PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix)); 1000 PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix)); 1001 /* propagate DM */ 1002 { 1003 DM sdm; 1004 PetscCall(KSPGetDM(jac->head->next->ksp, &sdm)); 1005 if (sdm) { 1006 PetscCall(KSPSetDM(jac->kspschur, sdm)); 1007 PetscCall(KSPSetDMActive(jac->kspschur, PETSC_FALSE)); 1008 } 1009 } 1010 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1011 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 1012 PetscCall(KSPSetFromOptions(jac->kspschur)); 1013 } 1014 PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY)); 1015 PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY)); 1016 1017 /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */ 1018 PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname)); 1019 PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L)); 1020 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L)); 1021 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L)); 1022 PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname)); 1023 PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L)); 1024 if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L)); 1025 if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L)); 1026 } else if (jac->type == PC_COMPOSITE_GKB) { 1027 IS ccis; 1028 PetscInt rstart, rend; 1029 1030 PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields"); 1031 1032 ilink = jac->head; 1033 1034 /* When extracting off-diagonal submatrices, we take complements from this range */ 1035 PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend)); 1036 1037 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 1038 if (jac->offdiag_use_amat) { 1039 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); 1040 } else { 1041 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); 1042 } 1043 PetscCall(ISDestroy(&ccis)); 1044 /* Create work vectors for GKB algorithm */ 1045 PetscCall(VecDuplicate(ilink->x, &jac->u)); 1046 PetscCall(VecDuplicate(ilink->x, &jac->Hu)); 1047 PetscCall(VecDuplicate(ilink->x, &jac->w2)); 1048 ilink = ilink->next; 1049 PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); 1050 if (jac->offdiag_use_amat) { 1051 PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C)); 1052 } else { 1053 PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C)); 1054 } 1055 PetscCall(ISDestroy(&ccis)); 1056 /* Create work vectors for GKB algorithm */ 1057 PetscCall(VecDuplicate(ilink->x, &jac->v)); 1058 PetscCall(VecDuplicate(ilink->x, &jac->d)); 1059 PetscCall(VecDuplicate(ilink->x, &jac->w1)); 1060 PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu)); 1061 PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz)); 1062 1063 ilink = jac->head; 1064 PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H)); 1065 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1066 /* Create gkb_monitor context */ 1067 if (jac->gkbmonitor) { 1068 PetscInt tablevel; 1069 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer)); 1070 PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII)); 1071 PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel)); 1072 PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel)); 1073 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1)); 1074 } 1075 } else { 1076 /* set up the individual splits' PCs */ 1077 i = 0; 1078 ilink = jac->head; 1079 while (ilink) { 1080 PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i])); 1081 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 1082 if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp)); 1083 i++; 1084 ilink = ilink->next; 1085 } 1086 } 1087 1088 /* Set coordinates to the sub PC objects whenever these are set */ 1089 if (jac->coordinates_set) { 1090 PC pc_coords; 1091 if (jac->type == PC_COMPOSITE_SCHUR) { 1092 // Head is first block. 1093 PetscCall(KSPGetPC(jac->head->ksp, &pc_coords)); 1094 PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords)); 1095 // Second one is Schur block, but its KSP object is in kspschur. 1096 PetscCall(KSPGetPC(jac->kspschur, &pc_coords)); 1097 PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords)); 1098 } else if (jac->type == PC_COMPOSITE_GKB) { 1099 PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n")); 1100 } else { 1101 ilink = jac->head; 1102 while (ilink) { 1103 PetscCall(KSPGetPC(ilink->ksp, &pc_coords)); 1104 PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords)); 1105 ilink = ilink->next; 1106 } 1107 } 1108 } 1109 1110 jac->suboptionsset = PETSC_TRUE; 1111 PetscFunctionReturn(PETSC_SUCCESS); 1112 } 1113 1114 #define FieldSplitSplitSolveAdd(ilink, xx, yy) \ 1115 ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \ 1116 KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \ 1117 VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE))) 1118 1119 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y) 1120 { 1121 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1122 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1123 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1124 1125 PetscFunctionBegin; 1126 switch (jac->schurfactorization) { 1127 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1128 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1129 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1130 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1131 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1132 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1133 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1134 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1135 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1136 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1137 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1138 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1139 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1140 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1141 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1142 PetscCall(VecScale(ilinkD->y, jac->schurscale)); 1143 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1144 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1145 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1146 break; 1147 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1148 /* [A00 0; A10 S], suitable for left preconditioning */ 1149 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1150 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1151 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1152 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1153 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1154 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1155 PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x)); 1156 PetscCall(VecScale(ilinkD->x, -1.)); 1157 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1158 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1159 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1160 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1161 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1162 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1163 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1164 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1165 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1166 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1167 break; 1168 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1169 /* [A00 A01; 0 S], suitable for right preconditioning */ 1170 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1171 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1172 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1173 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1174 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1175 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1176 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x)); 1177 PetscCall(VecScale(ilinkA->x, -1.)); 1178 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1179 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1180 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1181 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1182 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1183 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1184 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1185 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1186 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1187 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1188 break; 1189 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1190 /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */ 1191 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1192 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1193 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL)); 1194 PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y)); 1195 PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y)); 1196 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL)); 1197 PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x)); 1198 PetscCall(VecScale(ilinkD->x, -1.0)); 1199 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1200 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1201 1202 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1203 PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y)); 1204 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1205 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1206 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1207 1208 if (kspUpper == kspA) { 1209 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y)); 1210 PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y)); 1211 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1212 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1213 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1214 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1215 } else { 1216 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1217 PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y)); 1218 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1219 PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x)); 1220 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL)); 1221 PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z)); 1222 PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z)); 1223 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL)); 1224 PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z)); 1225 } 1226 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1227 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1228 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1229 } 1230 PetscFunctionReturn(PETSC_SUCCESS); 1231 } 1232 1233 static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y) 1234 { 1235 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1236 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1237 KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper; 1238 1239 PetscFunctionBegin; 1240 switch (jac->schurfactorization) { 1241 case PC_FIELDSPLIT_SCHUR_FACT_DIAG: 1242 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 1243 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1244 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1245 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1246 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1247 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1248 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1249 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1250 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1251 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1252 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1253 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1254 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1255 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1256 PetscCall(VecScale(ilinkD->y, jac->schurscale)); 1257 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1258 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1259 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1260 break; 1261 case PC_FIELDSPLIT_SCHUR_FACT_UPPER: 1262 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1263 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1264 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1265 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1266 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1267 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1268 PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x)); 1269 PetscCall(VecScale(ilinkD->x, -1.)); 1270 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1271 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1272 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1273 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1274 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1275 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1276 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1277 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1278 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1279 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1280 break; 1281 case PC_FIELDSPLIT_SCHUR_FACT_LOWER: 1282 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1283 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1284 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1285 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1286 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1287 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1288 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x)); 1289 PetscCall(VecScale(ilinkA->x, -1.)); 1290 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1291 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1292 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD)); 1293 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1294 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1295 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1296 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1297 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1298 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1299 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1300 break; 1301 case PC_FIELDSPLIT_SCHUR_FACT_FULL: 1302 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1303 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1304 PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL)); 1305 PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y)); 1306 PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y)); 1307 PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL)); 1308 PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x)); 1309 PetscCall(VecScale(ilinkD->x, -1.0)); 1310 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1311 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD)); 1312 1313 PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1314 PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y)); 1315 PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y)); 1316 PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL)); 1317 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1318 1319 if (kspLower == kspA) { 1320 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y)); 1321 PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y)); 1322 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1323 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1324 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1325 PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1326 } else { 1327 PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL)); 1328 PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y)); 1329 PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y)); 1330 PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x)); 1331 PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL)); 1332 PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z)); 1333 PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z)); 1334 PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL)); 1335 PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z)); 1336 } 1337 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1338 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1339 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1340 } 1341 PetscFunctionReturn(PETSC_SUCCESS); 1342 } 1343 1344 static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y) 1345 { 1346 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1347 PC_FieldSplitLink ilink = jac->head; 1348 PetscInt cnt, bs; 1349 1350 PetscFunctionBegin; 1351 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1352 if (jac->defaultsplit) { 1353 PetscCall(VecGetBlockSize(x, &bs)); 1354 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1355 PetscCall(VecGetBlockSize(y, &bs)); 1356 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1357 PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES)); 1358 while (ilink) { 1359 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1360 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1361 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1362 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1363 ilink = ilink->next; 1364 } 1365 PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES)); 1366 } else { 1367 PetscCall(VecSet(y, 0.0)); 1368 while (ilink) { 1369 PetscCall(FieldSplitSplitSolveAdd(ilink, x, y)); 1370 ilink = ilink->next; 1371 } 1372 } 1373 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) { 1374 PetscCall(VecSet(y, 0.0)); 1375 /* solve on first block for first block variables */ 1376 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD)); 1377 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD)); 1378 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1379 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1380 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1381 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1382 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1383 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1384 1385 /* compute the residual only onto second block variables using first block variables */ 1386 PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x)); 1387 ilink = ilink->next; 1388 PetscCall(VecScale(ilink->x, -1.0)); 1389 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1390 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1391 1392 /* solve on second block variables */ 1393 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1394 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1395 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1396 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1397 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1398 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1399 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1400 if (!jac->w1) { 1401 PetscCall(VecDuplicate(x, &jac->w1)); 1402 PetscCall(VecDuplicate(x, &jac->w2)); 1403 } 1404 PetscCall(VecSet(y, 0.0)); 1405 PetscCall(FieldSplitSplitSolveAdd(ilink, x, y)); 1406 cnt = 1; 1407 while (ilink->next) { 1408 ilink = ilink->next; 1409 /* compute the residual only over the part of the vector needed */ 1410 PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x)); 1411 PetscCall(VecScale(ilink->x, -1.0)); 1412 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1413 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1414 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1415 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1416 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1417 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1418 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1419 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1420 } 1421 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1422 cnt -= 2; 1423 while (ilink->previous) { 1424 ilink = ilink->previous; 1425 /* compute the residual only over the part of the vector needed */ 1426 PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x)); 1427 PetscCall(VecScale(ilink->x, -1.0)); 1428 PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1429 PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD)); 1430 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1431 PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y)); 1432 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1433 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1434 PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1435 PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE)); 1436 } 1437 } 1438 } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type); 1439 PetscFunctionReturn(PETSC_SUCCESS); 1440 } 1441 1442 static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y) 1443 { 1444 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1445 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 1446 KSP ksp = ilinkA->ksp; 1447 Vec u, v, Hu, d, work1, work2; 1448 PetscScalar alpha, z, nrmz2, *vecz; 1449 PetscReal lowbnd, nu, beta; 1450 PetscInt j, iterGKB; 1451 1452 PetscFunctionBegin; 1453 PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1454 PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1455 PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD)); 1456 PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD)); 1457 1458 u = jac->u; 1459 v = jac->v; 1460 Hu = jac->Hu; 1461 d = jac->d; 1462 work1 = jac->w1; 1463 work2 = jac->w2; 1464 vecz = jac->vecz; 1465 1466 /* Change RHS to comply with matrix regularization H = A + nu*B*B' */ 1467 /* Add q = q + nu*B*b */ 1468 if (jac->gkbnu) { 1469 nu = jac->gkbnu; 1470 PetscCall(VecScale(ilinkD->x, jac->gkbnu)); 1471 PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */ 1472 } else { 1473 /* Situation when no augmented Lagrangian is used. Then we set inner */ 1474 /* matrix N = I in [Ar13], and thus nu = 1. */ 1475 nu = 1; 1476 } 1477 1478 /* Transform rhs from [q,tilde{b}] to [0,b] */ 1479 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL)); 1480 PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y)); 1481 PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y)); 1482 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL)); 1483 PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1)); 1484 PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x */ 1485 1486 /* First step of algorithm */ 1487 PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/ 1488 KSPCheckDot(ksp, beta); 1489 beta = PetscSqrtReal(nu) * beta; 1490 PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c */ 1491 PetscCall(MatMult(jac->B, v, work2)); /* u = H^{-1}*B*v */ 1492 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL)); 1493 PetscCall(KSPSolve(ksp, work2, u)); 1494 PetscCall(KSPCheckSolve(ksp, pc, u)); 1495 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL)); 1496 PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */ 1497 PetscCall(VecDot(Hu, u, &alpha)); 1498 KSPCheckDot(ksp, alpha); 1499 PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite"); 1500 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1501 PetscCall(VecScale(u, 1.0 / alpha)); 1502 PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c */ 1503 1504 z = beta / alpha; 1505 vecz[1] = z; 1506 1507 /* Computation of first iterate x(1) and p(1) */ 1508 PetscCall(VecAXPY(ilinkA->y, z, u)); 1509 PetscCall(VecCopy(d, ilinkD->y)); 1510 PetscCall(VecScale(ilinkD->y, -z)); 1511 1512 iterGKB = 1; 1513 lowbnd = 2 * jac->gkbtol; 1514 if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd)); 1515 1516 while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) { 1517 iterGKB += 1; 1518 PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */ 1519 PetscCall(VecAXPBY(v, nu, -alpha, work1)); 1520 PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v */ 1521 beta = beta / PetscSqrtReal(nu); 1522 PetscCall(VecScale(v, 1.0 / beta)); 1523 PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */ 1524 PetscCall(MatMult(jac->H, u, Hu)); 1525 PetscCall(VecAXPY(work2, -beta, Hu)); 1526 PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL)); 1527 PetscCall(KSPSolve(ksp, work2, u)); 1528 PetscCall(KSPCheckSolve(ksp, pc, u)); 1529 PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL)); 1530 PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u */ 1531 PetscCall(VecDot(Hu, u, &alpha)); 1532 KSPCheckDot(ksp, alpha); 1533 PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite"); 1534 alpha = PetscSqrtReal(PetscAbsScalar(alpha)); 1535 PetscCall(VecScale(u, 1.0 / alpha)); 1536 1537 z = -beta / alpha * z; /* z <- beta/alpha*z */ 1538 vecz[0] = z; 1539 1540 /* Computation of new iterate x(i+1) and p(i+1) */ 1541 PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */ 1542 PetscCall(VecAXPY(ilinkA->y, z, u)); /* r = r + z*u */ 1543 PetscCall(VecAXPY(ilinkD->y, -z, d)); /* p = p - z*d */ 1544 PetscCall(MatMult(jac->H, ilinkA->y, Hu)); /* ||u||_H = u'*H*u */ 1545 PetscCall(VecDot(Hu, ilinkA->y, &nrmz2)); 1546 1547 /* Compute Lower Bound estimate */ 1548 if (iterGKB > jac->gkbdelay) { 1549 lowbnd = 0.0; 1550 for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]); 1551 lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2)); 1552 } 1553 1554 for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2]; 1555 if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd)); 1556 } 1557 1558 PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1559 PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1560 PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1561 PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE)); 1562 PetscFunctionReturn(PETSC_SUCCESS); 1563 } 1564 1565 #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \ 1566 ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \ 1567 KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \ 1568 VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE))) 1569 1570 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y) 1571 { 1572 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1573 PC_FieldSplitLink ilink = jac->head; 1574 PetscInt bs; 1575 1576 PetscFunctionBegin; 1577 if (jac->type == PC_COMPOSITE_ADDITIVE) { 1578 if (jac->defaultsplit) { 1579 PetscCall(VecGetBlockSize(x, &bs)); 1580 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1581 PetscCall(VecGetBlockSize(y, &bs)); 1582 PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs); 1583 PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES)); 1584 while (ilink) { 1585 PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1586 PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y)); 1587 PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y)); 1588 PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL)); 1589 ilink = ilink->next; 1590 } 1591 PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES)); 1592 } else { 1593 PetscCall(VecSet(y, 0.0)); 1594 while (ilink) { 1595 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 1596 ilink = ilink->next; 1597 } 1598 } 1599 } else { 1600 if (!jac->w1) { 1601 PetscCall(VecDuplicate(x, &jac->w1)); 1602 PetscCall(VecDuplicate(x, &jac->w2)); 1603 } 1604 PetscCall(VecSet(y, 0.0)); 1605 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 1606 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 1607 while (ilink->next) { 1608 ilink = ilink->next; 1609 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 1610 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 1611 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 1612 } 1613 while (ilink->previous) { 1614 ilink = ilink->previous; 1615 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 1616 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 1617 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 1618 } 1619 } else { 1620 while (ilink->next) { /* get to last entry in linked list */ 1621 ilink = ilink->next; 1622 } 1623 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y)); 1624 while (ilink->previous) { 1625 ilink = ilink->previous; 1626 PetscCall(MatMultTranspose(pc->mat, y, jac->w1)); 1627 PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x)); 1628 PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y)); 1629 } 1630 } 1631 } 1632 PetscFunctionReturn(PETSC_SUCCESS); 1633 } 1634 1635 static PetscErrorCode PCReset_FieldSplit(PC pc) 1636 { 1637 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1638 PC_FieldSplitLink ilink = jac->head, next; 1639 1640 PetscFunctionBegin; 1641 while (ilink) { 1642 PetscCall(KSPDestroy(&ilink->ksp)); 1643 PetscCall(VecDestroy(&ilink->x)); 1644 PetscCall(VecDestroy(&ilink->y)); 1645 PetscCall(VecDestroy(&ilink->z)); 1646 PetscCall(VecScatterDestroy(&ilink->sctx)); 1647 PetscCall(ISDestroy(&ilink->is)); 1648 PetscCall(ISDestroy(&ilink->is_col)); 1649 PetscCall(PetscFree(ilink->splitname)); 1650 PetscCall(PetscFree(ilink->fields)); 1651 PetscCall(PetscFree(ilink->fields_col)); 1652 next = ilink->next; 1653 PetscCall(PetscFree(ilink)); 1654 ilink = next; 1655 } 1656 jac->head = NULL; 1657 PetscCall(PetscFree2(jac->x, jac->y)); 1658 if (jac->mat && jac->mat != jac->pmat) { 1659 PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat)); 1660 } else if (jac->mat) { 1661 jac->mat = NULL; 1662 } 1663 if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat)); 1664 if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield)); 1665 jac->nsplits = 0; 1666 PetscCall(VecDestroy(&jac->w1)); 1667 PetscCall(VecDestroy(&jac->w2)); 1668 PetscCall(MatDestroy(&jac->schur)); 1669 PetscCall(MatDestroy(&jac->schurp)); 1670 PetscCall(MatDestroy(&jac->schur_user)); 1671 PetscCall(KSPDestroy(&jac->kspschur)); 1672 PetscCall(KSPDestroy(&jac->kspupper)); 1673 PetscCall(MatDestroy(&jac->B)); 1674 PetscCall(MatDestroy(&jac->C)); 1675 PetscCall(MatDestroy(&jac->H)); 1676 PetscCall(VecDestroy(&jac->u)); 1677 PetscCall(VecDestroy(&jac->v)); 1678 PetscCall(VecDestroy(&jac->Hu)); 1679 PetscCall(VecDestroy(&jac->d)); 1680 PetscCall(PetscFree(jac->vecz)); 1681 PetscCall(PetscViewerDestroy(&jac->gkbviewer)); 1682 jac->isrestrict = PETSC_FALSE; 1683 PetscFunctionReturn(PETSC_SUCCESS); 1684 } 1685 1686 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 1687 { 1688 PetscFunctionBegin; 1689 PetscCall(PCReset_FieldSplit(pc)); 1690 PetscCall(PetscFree(pc->data)); 1691 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL)); 1692 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL)); 1693 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL)); 1694 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL)); 1695 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL)); 1696 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL)); 1697 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL)); 1698 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL)); 1699 1700 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL)); 1701 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL)); 1702 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL)); 1703 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL)); 1704 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL)); 1705 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL)); 1706 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL)); 1707 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL)); 1708 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL)); 1709 PetscFunctionReturn(PETSC_SUCCESS); 1710 } 1711 1712 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject) 1713 { 1714 PetscInt bs; 1715 PetscBool flg; 1716 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1717 PCCompositeType ctype; 1718 1719 PetscFunctionBegin; 1720 PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options"); 1721 PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL)); 1722 PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg)); 1723 if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs)); 1724 jac->diag_use_amat = pc->useAmat; 1725 PetscCall(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)); 1726 jac->offdiag_use_amat = pc->useAmat; 1727 PetscCall(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)); 1728 PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL)); 1729 PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */ 1730 PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg)); 1731 if (flg) PetscCall(PCFieldSplitSetType(pc, ctype)); 1732 /* Only setup fields once */ 1733 if ((jac->bs > 0) && (jac->nsplits == 0)) { 1734 /* only allow user to set fields from command line if bs is already known. 1735 otherwise user can set them in PCFieldSplitSetDefaults() */ 1736 PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc)); 1737 if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n")); 1738 } 1739 if (jac->type == PC_COMPOSITE_SCHUR) { 1740 PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg)); 1741 if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n")); 1742 PetscCall(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)); 1743 PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL)); 1744 PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL)); 1745 } else if (jac->type == PC_COMPOSITE_GKB) { 1746 PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL)); 1747 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL)); 1748 PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0)); 1749 PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL)); 1750 PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL)); 1751 } 1752 /* 1753 In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet. 1754 But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it 1755 is called on the outer solver in case changes were made in the options database 1756 1757 But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions() 1758 if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete. 1759 Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types. 1760 1761 There could be a negative side effect of calling the KSPSetFromOptions() below. 1762 1763 If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call 1764 */ 1765 if (jac->issetup) { 1766 PC_FieldSplitLink ilink = jac->head; 1767 if (jac->type == PC_COMPOSITE_SCHUR) { 1768 if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper)); 1769 if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur)); 1770 } 1771 while (ilink) { 1772 if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp)); 1773 ilink = ilink->next; 1774 } 1775 } 1776 PetscOptionsHeadEnd(); 1777 PetscFunctionReturn(PETSC_SUCCESS); 1778 } 1779 1780 static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col) 1781 { 1782 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1783 PC_FieldSplitLink ilink, next = jac->head; 1784 char prefix[128]; 1785 PetscInt i; 1786 1787 PetscFunctionBegin; 1788 if (jac->splitdefined) { 1789 PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname)); 1790 PetscFunctionReturn(PETSC_SUCCESS); 1791 } 1792 for (i = 0; i < n; i++) { 1793 PetscCheck(fields[i] < jac->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", fields[i], jac->bs); 1794 PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]); 1795 } 1796 PetscCall(PetscNew(&ilink)); 1797 if (splitname) { 1798 PetscCall(PetscStrallocpy(splitname, &ilink->splitname)); 1799 } else { 1800 PetscCall(PetscMalloc1(3, &ilink->splitname)); 1801 PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits)); 1802 } 1803 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 */ 1804 PetscCall(PetscMalloc1(n, &ilink->fields)); 1805 PetscCall(PetscArraycpy(ilink->fields, fields, n)); 1806 PetscCall(PetscMalloc1(n, &ilink->fields_col)); 1807 PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n)); 1808 1809 ilink->nfields = n; 1810 ilink->next = NULL; 1811 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp)); 1812 PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel)); 1813 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure)); 1814 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1)); 1815 PetscCall(KSPSetType(ilink->ksp, KSPPREONLY)); 1816 1817 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 1818 PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix)); 1819 1820 if (!next) { 1821 jac->head = ilink; 1822 ilink->previous = NULL; 1823 } else { 1824 while (next->next) next = next->next; 1825 next->next = ilink; 1826 ilink->previous = next; 1827 } 1828 jac->nsplits++; 1829 PetscFunctionReturn(PETSC_SUCCESS); 1830 } 1831 1832 static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp) 1833 { 1834 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1835 1836 PetscFunctionBegin; 1837 *subksp = NULL; 1838 if (n) *n = 0; 1839 if (jac->type == PC_COMPOSITE_SCHUR) { 1840 PetscInt nn; 1841 1842 PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()"); 1843 PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits); 1844 nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0); 1845 PetscCall(PetscMalloc1(nn, subksp)); 1846 (*subksp)[0] = jac->head->ksp; 1847 (*subksp)[1] = jac->kspschur; 1848 if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper; 1849 if (n) *n = nn; 1850 } 1851 PetscFunctionReturn(PETSC_SUCCESS); 1852 } 1853 1854 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp) 1855 { 1856 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1857 1858 PetscFunctionBegin; 1859 PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()"); 1860 PetscCall(PetscMalloc1(jac->nsplits, subksp)); 1861 PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp)); 1862 1863 (*subksp)[1] = jac->kspschur; 1864 if (n) *n = jac->nsplits; 1865 PetscFunctionReturn(PETSC_SUCCESS); 1866 } 1867 1868 static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp) 1869 { 1870 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1871 PetscInt cnt = 0; 1872 PC_FieldSplitLink ilink = jac->head; 1873 1874 PetscFunctionBegin; 1875 PetscCall(PetscMalloc1(jac->nsplits, subksp)); 1876 while (ilink) { 1877 (*subksp)[cnt++] = ilink->ksp; 1878 ilink = ilink->next; 1879 } 1880 PetscCheck(cnt == jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt PCFIELDSPLIT object: number of splits in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, jac->nsplits); 1881 if (n) *n = jac->nsplits; 1882 PetscFunctionReturn(PETSC_SUCCESS); 1883 } 1884 1885 /*@C 1886 PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`. 1887 1888 Input Parameters: 1889 + pc - the preconditioner context 1890 - isy - the index set that defines the indices to which the fieldsplit is to be restricted 1891 1892 Level: advanced 1893 1894 Developer Notes: 1895 It seems the resulting `IS`s will not cover the entire space, so 1896 how can they define a convergent preconditioner? Needs explaining. 1897 1898 .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 1899 @*/ 1900 PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy) 1901 { 1902 PetscFunctionBegin; 1903 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1904 PetscValidHeaderSpecific(isy, IS_CLASSID, 2); 1905 PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy)); 1906 PetscFunctionReturn(PETSC_SUCCESS); 1907 } 1908 1909 static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy) 1910 { 1911 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1912 PC_FieldSplitLink ilink = jac->head, next; 1913 PetscInt localsize, size, sizez, i; 1914 const PetscInt *ind, *indz; 1915 PetscInt *indc, *indcz; 1916 PetscBool flg; 1917 1918 PetscFunctionBegin; 1919 PetscCall(ISGetLocalSize(isy, &localsize)); 1920 PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy))); 1921 size -= localsize; 1922 while (ilink) { 1923 IS isrl, isr; 1924 PC subpc; 1925 PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl)); 1926 PetscCall(ISGetLocalSize(isrl, &localsize)); 1927 PetscCall(PetscMalloc1(localsize, &indc)); 1928 PetscCall(ISGetIndices(isrl, &ind)); 1929 PetscCall(PetscArraycpy(indc, ind, localsize)); 1930 PetscCall(ISRestoreIndices(isrl, &ind)); 1931 PetscCall(ISDestroy(&isrl)); 1932 for (i = 0; i < localsize; i++) *(indc + i) += size; 1933 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr)); 1934 PetscCall(PetscObjectReference((PetscObject)isr)); 1935 PetscCall(ISDestroy(&ilink->is)); 1936 ilink->is = isr; 1937 PetscCall(PetscObjectReference((PetscObject)isr)); 1938 PetscCall(ISDestroy(&ilink->is_col)); 1939 ilink->is_col = isr; 1940 PetscCall(ISDestroy(&isr)); 1941 PetscCall(KSPGetPC(ilink->ksp, &subpc)); 1942 PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg)); 1943 if (flg) { 1944 IS iszl, isz; 1945 MPI_Comm comm; 1946 PetscCall(ISGetLocalSize(ilink->is, &localsize)); 1947 comm = PetscObjectComm((PetscObject)ilink->is); 1948 PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl)); 1949 PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm)); 1950 sizez -= localsize; 1951 PetscCall(ISGetLocalSize(iszl, &localsize)); 1952 PetscCall(PetscMalloc1(localsize, &indcz)); 1953 PetscCall(ISGetIndices(iszl, &indz)); 1954 PetscCall(PetscArraycpy(indcz, indz, localsize)); 1955 PetscCall(ISRestoreIndices(iszl, &indz)); 1956 PetscCall(ISDestroy(&iszl)); 1957 for (i = 0; i < localsize; i++) *(indcz + i) += sizez; 1958 PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz)); 1959 PetscCall(PCFieldSplitRestrictIS(subpc, isz)); 1960 PetscCall(ISDestroy(&isz)); 1961 } 1962 next = ilink->next; 1963 ilink = next; 1964 } 1965 jac->isrestrict = PETSC_TRUE; 1966 PetscFunctionReturn(PETSC_SUCCESS); 1967 } 1968 1969 static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is) 1970 { 1971 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 1972 PC_FieldSplitLink ilink, next = jac->head; 1973 char prefix[128]; 1974 1975 PetscFunctionBegin; 1976 if (jac->splitdefined) { 1977 PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname)); 1978 PetscFunctionReturn(PETSC_SUCCESS); 1979 } 1980 PetscCall(PetscNew(&ilink)); 1981 if (splitname) { 1982 PetscCall(PetscStrallocpy(splitname, &ilink->splitname)); 1983 } else { 1984 PetscCall(PetscMalloc1(8, &ilink->splitname)); 1985 PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits)); 1986 } 1987 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 */ 1988 PetscCall(PetscObjectReference((PetscObject)is)); 1989 PetscCall(ISDestroy(&ilink->is)); 1990 ilink->is = is; 1991 PetscCall(PetscObjectReference((PetscObject)is)); 1992 PetscCall(ISDestroy(&ilink->is_col)); 1993 ilink->is_col = is; 1994 ilink->next = NULL; 1995 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp)); 1996 PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel)); 1997 PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure)); 1998 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1)); 1999 PetscCall(KSPSetType(ilink->ksp, KSPPREONLY)); 2000 2001 PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname)); 2002 PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix)); 2003 2004 if (!next) { 2005 jac->head = ilink; 2006 ilink->previous = NULL; 2007 } else { 2008 while (next->next) next = next->next; 2009 next->next = ilink; 2010 ilink->previous = next; 2011 } 2012 jac->nsplits++; 2013 PetscFunctionReturn(PETSC_SUCCESS); 2014 } 2015 2016 /*@C 2017 PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT` 2018 2019 Logically Collective 2020 2021 Input Parameters: 2022 + pc - the preconditioner context 2023 . splitname - name of this split, if `NULL` the number of the split is used 2024 . n - the number of fields in this split 2025 . fields - the fields in this split 2026 - fields_col - generally the same as fields, if it does not match fields then the matrix block that is solved for this set of fields comes from an off-diagonal block 2027 of the matrix and fields_col provides the column indices for that block 2028 2029 Level: intermediate 2030 2031 Notes: 2032 Use `PCFieldSplitSetIS()` to set a general set of indices as a split. 2033 2034 `PCFieldSplitSetFields()` is for defining fields as strided blocks. For example, if the block 2035 size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean 2036 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 2037 where the numbered entries indicate what is in the split. 2038 2039 This function is called once per split (it creates a new split each time). Solve options 2040 for this split will be available under the prefix `-fieldsplit_SPLITNAME_`. 2041 2042 `PCFieldSplitSetIS()` does not support having a fields_col different from fields 2043 2044 Developer Notes: 2045 This routine does not actually create the `IS` representing the split, that is delayed 2046 until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be 2047 available when this routine is called. 2048 2049 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()` 2050 @*/ 2051 PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col) 2052 { 2053 PetscFunctionBegin; 2054 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2055 PetscAssertPointer(splitname, 2); 2056 PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname); 2057 PetscAssertPointer(fields, 4); 2058 PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col)); 2059 PetscFunctionReturn(PETSC_SUCCESS); 2060 } 2061 2062 /*@ 2063 PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build 2064 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2065 2066 Logically Collective 2067 2068 Input Parameters: 2069 + pc - the preconditioner object 2070 - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 2071 2072 Options Database Key: 2073 . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks 2074 2075 Level: intermediate 2076 2077 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT` 2078 @*/ 2079 PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg) 2080 { 2081 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2082 PetscBool isfs; 2083 2084 PetscFunctionBegin; 2085 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2086 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2087 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2088 jac->diag_use_amat = flg; 2089 PetscFunctionReturn(PETSC_SUCCESS); 2090 } 2091 2092 /*@ 2093 PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build 2094 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2095 2096 Logically Collective 2097 2098 Input Parameter: 2099 . pc - the preconditioner object 2100 2101 Output Parameter: 2102 . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from 2103 2104 Level: intermediate 2105 2106 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT` 2107 @*/ 2108 PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg) 2109 { 2110 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2111 PetscBool isfs; 2112 2113 PetscFunctionBegin; 2114 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2115 PetscAssertPointer(flg, 2); 2116 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2117 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2118 *flg = jac->diag_use_amat; 2119 PetscFunctionReturn(PETSC_SUCCESS); 2120 } 2121 2122 /*@ 2123 PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build 2124 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2125 2126 Logically Collective 2127 2128 Input Parameters: 2129 + pc - the preconditioner object 2130 - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2131 2132 Options Database Key: 2133 . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks 2134 2135 Level: intermediate 2136 2137 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT` 2138 @*/ 2139 PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg) 2140 { 2141 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2142 PetscBool isfs; 2143 2144 PetscFunctionBegin; 2145 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2146 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2147 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2148 jac->offdiag_use_amat = flg; 2149 PetscFunctionReturn(PETSC_SUCCESS); 2150 } 2151 2152 /*@ 2153 PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build 2154 the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators. 2155 2156 Logically Collective 2157 2158 Input Parameter: 2159 . pc - the preconditioner object 2160 2161 Output Parameter: 2162 . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from 2163 2164 Level: intermediate 2165 2166 .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT` 2167 @*/ 2168 PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg) 2169 { 2170 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2171 PetscBool isfs; 2172 2173 PetscFunctionBegin; 2174 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2175 PetscAssertPointer(flg, 2); 2176 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 2177 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT); 2178 *flg = jac->offdiag_use_amat; 2179 PetscFunctionReturn(PETSC_SUCCESS); 2180 } 2181 2182 /*@C 2183 PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT` 2184 2185 Logically Collective 2186 2187 Input Parameters: 2188 + pc - the preconditioner context 2189 . splitname - name of this split, if `NULL` the number of the split is used 2190 - is - the index set that defines the elements in this split 2191 2192 Level: intermediate 2193 2194 Notes: 2195 Use `PCFieldSplitSetFields()`, for splits defined by strided types. 2196 2197 This function is called once per split (it creates a new split each time). Solve options 2198 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 2199 2200 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()` 2201 @*/ 2202 PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is) 2203 { 2204 PetscFunctionBegin; 2205 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2206 if (splitname) PetscAssertPointer(splitname, 2); 2207 PetscValidHeaderSpecific(is, IS_CLASSID, 3); 2208 PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is)); 2209 PetscFunctionReturn(PETSC_SUCCESS); 2210 } 2211 2212 /*@C 2213 PCFieldSplitGetIS - Retrieves the elements for a split as an `IS` 2214 2215 Logically Collective 2216 2217 Input Parameters: 2218 + pc - the preconditioner context 2219 - splitname - name of this split 2220 2221 Output Parameter: 2222 . is - the index set that defines the elements in this split, or `NULL` if the split is not found 2223 2224 Level: intermediate 2225 2226 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()` 2227 @*/ 2228 PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is) 2229 { 2230 PetscFunctionBegin; 2231 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2232 PetscAssertPointer(splitname, 2); 2233 PetscAssertPointer(is, 3); 2234 { 2235 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2236 PC_FieldSplitLink ilink = jac->head; 2237 PetscBool found; 2238 2239 *is = NULL; 2240 while (ilink) { 2241 PetscCall(PetscStrcmp(ilink->splitname, splitname, &found)); 2242 if (found) { 2243 *is = ilink->is; 2244 break; 2245 } 2246 ilink = ilink->next; 2247 } 2248 } 2249 PetscFunctionReturn(PETSC_SUCCESS); 2250 } 2251 2252 /*@C 2253 PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS` 2254 2255 Logically Collective 2256 2257 Input Parameters: 2258 + pc - the preconditioner context 2259 - index - index of this split 2260 2261 Output Parameter: 2262 . is - the index set that defines the elements in this split 2263 2264 Level: intermediate 2265 2266 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()` 2267 @*/ 2268 PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is) 2269 { 2270 PetscFunctionBegin; 2271 PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index); 2272 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2273 PetscAssertPointer(is, 3); 2274 { 2275 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2276 PC_FieldSplitLink ilink = jac->head; 2277 PetscInt i = 0; 2278 PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits); 2279 2280 while (i < index) { 2281 ilink = ilink->next; 2282 ++i; 2283 } 2284 PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is)); 2285 } 2286 PetscFunctionReturn(PETSC_SUCCESS); 2287 } 2288 2289 /*@ 2290 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 2291 fieldsplit preconditioner when calling `PCFieldSplitSetIS()`. If not set the matrix block size is used. 2292 2293 Logically Collective 2294 2295 Input Parameters: 2296 + pc - the preconditioner context 2297 - bs - the block size 2298 2299 Level: intermediate 2300 2301 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()` 2302 @*/ 2303 PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs) 2304 { 2305 PetscFunctionBegin; 2306 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2307 PetscValidLogicalCollectiveInt(pc, bs, 2); 2308 PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs)); 2309 PetscFunctionReturn(PETSC_SUCCESS); 2310 } 2311 2312 /*@C 2313 PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits 2314 2315 Collective 2316 2317 Input Parameter: 2318 . pc - the preconditioner context 2319 2320 Output Parameters: 2321 + n - the number of splits 2322 - subksp - the array of `KSP` contexts 2323 2324 Level: advanced 2325 2326 Notes: 2327 After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()` 2328 (not the `KSP`, just the array that contains them). 2329 2330 You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`. 2331 2332 If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the 2333 Schur complement and the `KSP` object used to iterate over the Schur complement. 2334 To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`. 2335 2336 If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the 2337 inner linear system defined by the matrix H in each loop. 2338 2339 Fortran Notes: 2340 You must pass in a `KSP` array that is large enough to contain all the `KSP`s. 2341 You can call `PCFieldSplitGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the 2342 `KSP` array must be. 2343 2344 Developer Notes: 2345 There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()` 2346 2347 The Fortran interface should be modernized to return directly the array of values. 2348 2349 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()` 2350 @*/ 2351 PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[]) 2352 { 2353 PetscFunctionBegin; 2354 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2355 if (n) PetscAssertPointer(n, 2); 2356 PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp)); 2357 PetscFunctionReturn(PETSC_SUCCESS); 2358 } 2359 2360 /*@C 2361 PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT` 2362 2363 Collective 2364 2365 Input Parameter: 2366 . pc - the preconditioner context 2367 2368 Output Parameters: 2369 + n - the number of splits 2370 - subksp - the array of `KSP` contexts 2371 2372 Level: advanced 2373 2374 Notes: 2375 After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()` 2376 (not the `KSP` just the array that contains them). 2377 2378 You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`. 2379 2380 If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order) 2381 + 1 - the `KSP` used for the (1,1) block 2382 . 2 - the `KSP` used for the Schur complement (not the one used for the interior Schur solver) 2383 - 3 - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block). 2384 2385 It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`. 2386 2387 Fortran Notes: 2388 You must pass in a `KSP` array that is large enough to contain all the local `KSP`s. 2389 You can call `PCFieldSplitSchurGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the 2390 `KSP` array must be. 2391 2392 Developer Notes: 2393 There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()` 2394 2395 Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged? 2396 2397 The Fortran interface should be modernized to return directly the array of values. 2398 2399 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()` 2400 @*/ 2401 PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[]) 2402 { 2403 PetscFunctionBegin; 2404 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2405 if (n) PetscAssertPointer(n, 2); 2406 PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp)); 2407 PetscFunctionReturn(PETSC_SUCCESS); 2408 } 2409 2410 /*@ 2411 PCFieldSplitSetSchurPre - Indicates from what operator the preconditioner is constructed for the Schur complement. 2412 The default is the A11 matrix. 2413 2414 Collective 2415 2416 Input Parameters: 2417 + pc - the preconditioner context 2418 . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default), 2419 `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`, 2420 `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL` 2421 - pre - matrix to use for preconditioning, or `NULL` 2422 2423 Options Database Keys: 2424 + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments 2425 - -fieldsplit_1_pc_type <pctype> - the preconditioner algorithm that is used to construct the preconditioner from the operator 2426 2427 Level: intermediate 2428 2429 Notes: 2430 If ptype is 2431 + a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner 2432 matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix 2433 . self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix: 2434 The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM` 2435 . user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument 2436 to this function). 2437 . selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $ 2438 This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be 2439 lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump` 2440 - full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation 2441 computed internally by `PCFIELDSPLIT` (this is expensive) 2442 useful mostly as a test that the Schur complement approach can work for your problem 2443 2444 When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense 2445 with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and 2446 `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement. 2447 2448 Developer Note: 2449 The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere. 2450 2451 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, 2452 `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()` 2453 @*/ 2454 PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2455 { 2456 PetscFunctionBegin; 2457 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2458 PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre)); 2459 PetscFunctionReturn(PETSC_SUCCESS); 2460 } 2461 2462 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2463 { 2464 return PCFieldSplitSetSchurPre(pc, ptype, pre); 2465 } /* Deprecated name */ 2466 2467 /*@ 2468 PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be 2469 preconditioned. See `PCFieldSplitSetSchurPre()` for details. 2470 2471 Logically Collective 2472 2473 Input Parameter: 2474 . pc - the preconditioner context 2475 2476 Output Parameters: 2477 + ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER` 2478 - pre - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL` 2479 2480 Level: intermediate 2481 2482 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC` 2483 @*/ 2484 PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre) 2485 { 2486 PetscFunctionBegin; 2487 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2488 PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre)); 2489 PetscFunctionReturn(PETSC_SUCCESS); 2490 } 2491 2492 /*@ 2493 PCFieldSplitSchurGetS - extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately 2494 2495 Not Collective 2496 2497 Input Parameter: 2498 . pc - the preconditioner context 2499 2500 Output Parameter: 2501 . S - the Schur complement matrix 2502 2503 Level: advanced 2504 2505 Note: 2506 This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`. 2507 2508 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`, 2509 `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()` 2510 @*/ 2511 PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S) 2512 { 2513 const char *t; 2514 PetscBool isfs; 2515 PC_FieldSplit *jac; 2516 2517 PetscFunctionBegin; 2518 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2519 PetscCall(PetscObjectGetType((PetscObject)pc, &t)); 2520 PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs)); 2521 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t); 2522 jac = (PC_FieldSplit *)pc->data; 2523 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type); 2524 if (S) *S = jac->schur; 2525 PetscFunctionReturn(PETSC_SUCCESS); 2526 } 2527 2528 /*@ 2529 PCFieldSplitSchurRestoreS - returns the `MATSCHURCOMPLEMENT` matrix used by this `PC` 2530 2531 Not Collective 2532 2533 Input Parameters: 2534 + pc - the preconditioner context 2535 - S - the Schur complement matrix 2536 2537 Level: advanced 2538 2539 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()` 2540 @*/ 2541 PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S) 2542 { 2543 const char *t; 2544 PetscBool isfs; 2545 PC_FieldSplit *jac; 2546 2547 PetscFunctionBegin; 2548 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2549 PetscCall(PetscObjectGetType((PetscObject)pc, &t)); 2550 PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs)); 2551 PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t); 2552 jac = (PC_FieldSplit *)pc->data; 2553 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type); 2554 PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten"); 2555 PetscFunctionReturn(PETSC_SUCCESS); 2556 } 2557 2558 static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre) 2559 { 2560 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2561 2562 PetscFunctionBegin; 2563 jac->schurpre = ptype; 2564 if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) { 2565 PetscCall(MatDestroy(&jac->schur_user)); 2566 jac->schur_user = pre; 2567 PetscCall(PetscObjectReference((PetscObject)jac->schur_user)); 2568 } 2569 PetscFunctionReturn(PETSC_SUCCESS); 2570 } 2571 2572 static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre) 2573 { 2574 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2575 2576 PetscFunctionBegin; 2577 if (ptype) *ptype = jac->schurpre; 2578 if (pre) *pre = jac->schur_user; 2579 PetscFunctionReturn(PETSC_SUCCESS); 2580 } 2581 2582 /*@ 2583 PCFieldSplitSetSchurFactType - sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note` 2584 2585 Collective 2586 2587 Input Parameters: 2588 + pc - the preconditioner context 2589 - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default 2590 2591 Options Database Key: 2592 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full` 2593 2594 Level: intermediate 2595 2596 Notes: 2597 The FULL factorization is 2598 2599 ```{math} 2600 \left(\begin{array}{cc} A & B \\ 2601 C & E \\ 2602 \end{array}\right) = 2603 \left(\begin{array}{cc} 1 & 0 \\ 2604 C*A^{-1} & I \\ 2605 \end{array}\right) 2606 \left(\begin{array}{cc} A & 0 \\ 2607 0 & S \\ 2608 \end{array}\right) 2609 \left(\begin{array}{cc} I & A^{-1}B \\ 2610 0 & I \\ 2611 \end{array}\right) = L D U. 2612 ``` 2613 2614 where $ S = E - C*A^{-1}*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$, 2615 and DIAG is the diagonal part with the sign of $ S $ flipped (because this makes the preconditioner positive definite for many formulations, 2616 thus allowing the use of `KSPMINRES)`. Sign flipping of $ S $ can be turned off with `PCFieldSplitSetSchurScale()`. 2617 2618 If $A$ and $S$ are solved exactly 2619 + 1 - FULL factorization is a direct solver. 2620 . 2 - 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. 2621 - 3 - 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. 2622 2623 If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner 2624 application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice. 2625 2626 For symmetric problems in which $A$ is positive definite and $S$ is negative definite, DIAG can be used with `KSPMINRES`. 2627 2628 A flexible method like `KSPFGMRES` or `KSPGCR`, [](sec_flexibleksp), 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). 2629 2630 .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`, 2631 [](sec_flexibleksp), `PCFieldSplitSetSchurPre()` 2632 @*/ 2633 PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype) 2634 { 2635 PetscFunctionBegin; 2636 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2637 PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype)); 2638 PetscFunctionReturn(PETSC_SUCCESS); 2639 } 2640 2641 static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype) 2642 { 2643 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2644 2645 PetscFunctionBegin; 2646 jac->schurfactorization = ftype; 2647 PetscFunctionReturn(PETSC_SUCCESS); 2648 } 2649 2650 /*@ 2651 PCFieldSplitSetSchurScale - Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`. 2652 2653 Collective 2654 2655 Input Parameters: 2656 + pc - the preconditioner context 2657 - scale - scaling factor for the Schur complement 2658 2659 Options Database Key: 2660 . -pc_fieldsplit_schur_scale <scale> - default is -1.0 2661 2662 Level: intermediate 2663 2664 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()` 2665 @*/ 2666 PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale) 2667 { 2668 PetscFunctionBegin; 2669 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2670 PetscValidLogicalCollectiveScalar(pc, scale, 2); 2671 PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale)); 2672 PetscFunctionReturn(PETSC_SUCCESS); 2673 } 2674 2675 static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale) 2676 { 2677 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2678 2679 PetscFunctionBegin; 2680 jac->schurscale = scale; 2681 PetscFunctionReturn(PETSC_SUCCESS); 2682 } 2683 2684 /*@C 2685 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 2686 2687 Collective 2688 2689 Input Parameter: 2690 . pc - the preconditioner context 2691 2692 Output Parameters: 2693 + A00 - the (0,0) block 2694 . A01 - the (0,1) block 2695 . A10 - the (1,0) block 2696 - A11 - the (1,1) block 2697 2698 Level: advanced 2699 2700 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()` 2701 @*/ 2702 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11) 2703 { 2704 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2705 2706 PetscFunctionBegin; 2707 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2708 PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 2709 if (A00) *A00 = jac->pmat[0]; 2710 if (A01) *A01 = jac->B; 2711 if (A10) *A10 = jac->C; 2712 if (A11) *A11 = jac->pmat[1]; 2713 PetscFunctionReturn(PETSC_SUCCESS); 2714 } 2715 2716 /*@ 2717 PCFieldSplitSetGKBTol - Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT` 2718 2719 Collective 2720 2721 Input Parameters: 2722 + pc - the preconditioner context 2723 - tolerance - the solver tolerance 2724 2725 Options Database Key: 2726 . -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5 2727 2728 Level: intermediate 2729 2730 Note: 2731 The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion. 2732 It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than 2733 this estimate, the stopping criterion is satisfactory in practical cases. 2734 2735 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()` 2736 @*/ 2737 PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance) 2738 { 2739 PetscFunctionBegin; 2740 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2741 PetscValidLogicalCollectiveReal(pc, tolerance, 2); 2742 PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance)); 2743 PetscFunctionReturn(PETSC_SUCCESS); 2744 } 2745 2746 static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance) 2747 { 2748 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2749 2750 PetscFunctionBegin; 2751 jac->gkbtol = tolerance; 2752 PetscFunctionReturn(PETSC_SUCCESS); 2753 } 2754 2755 /*@ 2756 PCFieldSplitSetGKBMaxit - Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner in `PCFIELDSPLIT` 2757 2758 Collective 2759 2760 Input Parameters: 2761 + pc - the preconditioner context 2762 - maxit - the maximum number of iterations 2763 2764 Options Database Key: 2765 . -pc_fieldsplit_gkb_maxit <maxit> - default is 100 2766 2767 Level: intermediate 2768 2769 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()` 2770 @*/ 2771 PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit) 2772 { 2773 PetscFunctionBegin; 2774 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2775 PetscValidLogicalCollectiveInt(pc, maxit, 2); 2776 PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit)); 2777 PetscFunctionReturn(PETSC_SUCCESS); 2778 } 2779 2780 static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit) 2781 { 2782 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2783 2784 PetscFunctionBegin; 2785 jac->gkbmaxit = maxit; 2786 PetscFunctionReturn(PETSC_SUCCESS); 2787 } 2788 2789 /*@ 2790 PCFieldSplitSetGKBDelay - Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT` 2791 preconditioner. 2792 2793 Collective 2794 2795 Input Parameters: 2796 + pc - the preconditioner context 2797 - delay - the delay window in the lower bound estimate 2798 2799 Options Database Key: 2800 . -pc_fieldsplit_gkb_delay <delay> - default is 5 2801 2802 Level: intermediate 2803 2804 Notes: 2805 The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error $ ||u-u^k||_H $ 2806 is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs 2807 at least (`delay` + 1) iterations to stop. 2808 2809 For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013` 2810 2811 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2812 @*/ 2813 PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay) 2814 { 2815 PetscFunctionBegin; 2816 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2817 PetscValidLogicalCollectiveInt(pc, delay, 2); 2818 PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay)); 2819 PetscFunctionReturn(PETSC_SUCCESS); 2820 } 2821 2822 static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay) 2823 { 2824 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2825 2826 PetscFunctionBegin; 2827 jac->gkbdelay = delay; 2828 PetscFunctionReturn(PETSC_SUCCESS); 2829 } 2830 2831 /*@ 2832 PCFieldSplitSetGKBNu - Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the 2833 Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT` 2834 2835 Collective 2836 2837 Input Parameters: 2838 + pc - the preconditioner context 2839 - nu - the shift parameter 2840 2841 Options Database Key: 2842 . -pc_fieldsplit_gkb_nu <nu> - default is 1 2843 2844 Level: intermediate 2845 2846 Notes: 2847 This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing `nu` sufficiently large. However, 2848 if `nu` is chosen too large, the matrix H might be badly conditioned and the solution of the linear system $Hx = b$ in the inner loop becomes difficult. It is therefore 2849 necessary to find a good balance in between the convergence of the inner and outer loop. 2850 2851 For `nu` = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in {cite}`arioli2013` is then chosen as identity. 2852 2853 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()` 2854 @*/ 2855 PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu) 2856 { 2857 PetscFunctionBegin; 2858 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2859 PetscValidLogicalCollectiveReal(pc, nu, 2); 2860 PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu)); 2861 PetscFunctionReturn(PETSC_SUCCESS); 2862 } 2863 2864 static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu) 2865 { 2866 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2867 2868 PetscFunctionBegin; 2869 jac->gkbnu = nu; 2870 PetscFunctionReturn(PETSC_SUCCESS); 2871 } 2872 2873 static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type) 2874 { 2875 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2876 2877 PetscFunctionBegin; 2878 jac->type = type; 2879 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL)); 2880 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL)); 2881 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL)); 2882 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL)); 2883 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL)); 2884 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL)); 2885 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL)); 2886 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL)); 2887 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL)); 2888 2889 if (type == PC_COMPOSITE_SCHUR) { 2890 pc->ops->apply = PCApply_FieldSplit_Schur; 2891 pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur; 2892 pc->ops->view = PCView_FieldSplit_Schur; 2893 2894 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur)); 2895 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit)); 2896 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit)); 2897 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit)); 2898 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit)); 2899 } else if (type == PC_COMPOSITE_GKB) { 2900 pc->ops->apply = PCApply_FieldSplit_GKB; 2901 pc->ops->view = PCView_FieldSplit_GKB; 2902 2903 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit)); 2904 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit)); 2905 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit)); 2906 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit)); 2907 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit)); 2908 } else { 2909 pc->ops->apply = PCApply_FieldSplit; 2910 pc->ops->view = PCView_FieldSplit; 2911 2912 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit)); 2913 } 2914 PetscFunctionReturn(PETSC_SUCCESS); 2915 } 2916 2917 static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs) 2918 { 2919 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2920 2921 PetscFunctionBegin; 2922 PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs); 2923 PetscCheck(jac->bs <= 0 || jac->bs == bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change fieldsplit blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", jac->bs, bs); 2924 jac->bs = bs; 2925 PetscFunctionReturn(PETSC_SUCCESS); 2926 } 2927 2928 static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[]) 2929 { 2930 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 2931 PC_FieldSplitLink ilink_current = jac->head; 2932 IS is_owned; 2933 2934 PetscFunctionBegin; 2935 jac->coordinates_set = PETSC_TRUE; // Internal flag 2936 PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL)); 2937 2938 while (ilink_current) { 2939 // For each IS, embed it to get local coords indces 2940 IS is_coords; 2941 PetscInt ndofs_block; 2942 const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block 2943 2944 // Setting drop to true for safety. It should make no difference. 2945 PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords)); 2946 PetscCall(ISGetLocalSize(is_coords, &ndofs_block)); 2947 PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration)); 2948 2949 // Allocate coordinates vector and set it directly 2950 PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords)); 2951 for (PetscInt dof = 0; dof < ndofs_block; ++dof) { 2952 for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d]; 2953 } 2954 ilink_current->dim = dim; 2955 ilink_current->ndofs = ndofs_block; 2956 PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration)); 2957 PetscCall(ISDestroy(&is_coords)); 2958 ilink_current = ilink_current->next; 2959 } 2960 PetscCall(ISDestroy(&is_owned)); 2961 PetscFunctionReturn(PETSC_SUCCESS); 2962 } 2963 2964 /*@ 2965 PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT` 2966 2967 Collective 2968 2969 Input Parameters: 2970 + pc - the preconditioner context 2971 - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, 2972 `PC_COMPOSITE_GKB` 2973 2974 Options Database Key: 2975 . -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 2976 2977 Level: intermediate 2978 2979 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, 2980 `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()` 2981 @*/ 2982 PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type) 2983 { 2984 PetscFunctionBegin; 2985 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2986 PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type)); 2987 PetscFunctionReturn(PETSC_SUCCESS); 2988 } 2989 2990 /*@ 2991 PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT` 2992 2993 Not collective 2994 2995 Input Parameter: 2996 . pc - the preconditioner context 2997 2998 Output Parameter: 2999 . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR` 3000 3001 Level: intermediate 3002 3003 .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`, 3004 `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR` 3005 @*/ 3006 PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type) 3007 { 3008 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3009 3010 PetscFunctionBegin; 3011 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3012 PetscAssertPointer(type, 2); 3013 *type = jac->type; 3014 PetscFunctionReturn(PETSC_SUCCESS); 3015 } 3016 3017 /*@ 3018 PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible. 3019 3020 Logically Collective 3021 3022 Input Parameters: 3023 + pc - the preconditioner context 3024 - flg - boolean indicating whether to use field splits defined by the `DM` 3025 3026 Options Database Key: 3027 . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM` 3028 3029 Level: intermediate 3030 3031 Developer Note: 3032 The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database 3033 3034 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()` 3035 @*/ 3036 PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg) 3037 { 3038 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3039 PetscBool isfs; 3040 3041 PetscFunctionBegin; 3042 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3043 PetscValidLogicalCollectiveBool(pc, flg, 2); 3044 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 3045 if (isfs) jac->dm_splits = flg; 3046 PetscFunctionReturn(PETSC_SUCCESS); 3047 } 3048 3049 /*@ 3050 PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible. 3051 3052 Logically Collective 3053 3054 Input Parameter: 3055 . pc - the preconditioner context 3056 3057 Output Parameter: 3058 . flg - boolean indicating whether to use field splits defined by the `DM` 3059 3060 Level: intermediate 3061 3062 Developer Note: 3063 The name should be `PCFieldSplitGetUseDMSplits()` 3064 3065 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldsplitSetIS()` 3066 @*/ 3067 PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg) 3068 { 3069 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3070 PetscBool isfs; 3071 3072 PetscFunctionBegin; 3073 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 3074 PetscAssertPointer(flg, 2); 3075 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs)); 3076 if (isfs) { 3077 if (flg) *flg = jac->dm_splits; 3078 } 3079 PetscFunctionReturn(PETSC_SUCCESS); 3080 } 3081 3082 /*@ 3083 PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries. 3084 3085 Logically Collective 3086 3087 Input Parameter: 3088 . pc - the preconditioner context 3089 3090 Output Parameter: 3091 . flg - boolean indicating whether to detect fields or not 3092 3093 Level: intermediate 3094 3095 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()` 3096 @*/ 3097 PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg) 3098 { 3099 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3100 3101 PetscFunctionBegin; 3102 *flg = jac->detect; 3103 PetscFunctionReturn(PETSC_SUCCESS); 3104 } 3105 3106 /*@ 3107 PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries. 3108 3109 Logically Collective 3110 3111 Input Parameter: 3112 . pc - the preconditioner context 3113 3114 Output Parameter: 3115 . flg - boolean indicating whether to detect fields or not 3116 3117 Options Database Key: 3118 . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point 3119 3120 Level: intermediate 3121 3122 Note: 3123 Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`). 3124 3125 .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF` 3126 @*/ 3127 PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg) 3128 { 3129 PC_FieldSplit *jac = (PC_FieldSplit *)pc->data; 3130 3131 PetscFunctionBegin; 3132 jac->detect = flg; 3133 if (jac->detect) { 3134 PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR)); 3135 PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL)); 3136 } 3137 PetscFunctionReturn(PETSC_SUCCESS); 3138 } 3139 3140 /*MC 3141 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 3142 collections of variables (that may overlap) called splits. See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details. 3143 3144 Options Database Keys: 3145 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the `%d`'th split 3146 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 3147 been supplied explicitly by `-pc_fieldsplit_%d_fields` 3148 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 3149 . -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting 3150 . -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`; see `PCFieldSplitSetSchurPre()` 3151 . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - set factorization type when using `-pc_fieldsplit_type schur`; 3152 see `PCFieldSplitSetSchurFactType()` 3153 . -pc_fieldsplit_dm_splits <true,false> (default is true) - Whether to use `DMCreateFieldDecomposition()` for splits 3154 - -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver 3155 3156 Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` . 3157 The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_` 3158 For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields. 3159 3160 To set options on the solvers for each block append `-fieldsplit_` to all the `PC` 3161 options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1` 3162 3163 To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()` 3164 and set the options directly on the resulting `KSP` object 3165 3166 Level: intermediate 3167 3168 Notes: 3169 Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries and `PCFieldSplitSetIS()` 3170 to define a split by an arbitrary collection of entries. 3171 3172 If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports 3173 `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise, the splits are defined by entries strided by bs, 3174 beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`, 3175 if this is not called the block size defaults to the blocksize of the second matrix passed 3176 to `KSPSetOperators()`/`PCSetOperators()`. 3177 3178 For the Schur complement preconditioner if 3179 3180 ```{math} 3181 J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right] 3182 ``` 3183 3184 the preconditioner using `full` factorization is logically 3185 ```{math} 3186 \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right] 3187 ``` 3188 where the action of $\text{inv}(A_{00})$ is applied using the KSP solver with prefix `-fieldsplit_0_`. $S$ is the Schur complement 3189 ```{math} 3190 S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01} 3191 ``` 3192 which is usually dense and not stored explicitly. The action of $\text{ksp}(S)$ is computed using the KSP solver with prefix `-fieldsplit_splitname_` (where `splitname` was given 3193 in providing the SECOND split or 1 if not given). For `PCFieldSplitGetSubKSP()` when field number is 0, 3194 it returns the KSP associated with `-fieldsplit_0_` while field number 1 gives `-fieldsplit_1_` KSP. By default 3195 $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$. 3196 3197 The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above, 3198 `diag` gives 3199 ```{math} 3200 \left[\begin{array}{cc} \text{inv}(A_{00}) & 0 \\ 0 & -\text{ksp}(S) \end{array}\right] 3201 ``` 3202 Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$ so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip 3203 can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of 3204 ```{math} 3205 \left[\begin{array}{cc} A_{00} & 0 \\ A_{10} & S \end{array}\right] 3206 ``` 3207 where the inverses of A_{00} and S are applied using KSPs. The upper factorization is the inverse of 3208 ```{math} 3209 \left[\begin{array}{cc} A_{00} & A_{01} \\ 0 & S \end{array}\right] 3210 ``` 3211 where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s. 3212 3213 If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS` 3214 is used automatically for a second block. 3215 3216 The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1. 3217 Generally it should be used with the `MATAIJ` format. 3218 3219 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 3220 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`. 3221 One can also use `PCFIELDSPLIT` 3222 inside a smoother resulting in "Distributive Smoothers". 3223 3224 See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`. 3225 3226 The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the 3227 residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables. 3228 3229 The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape 3230 ```{math} 3231 \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right] 3232 ``` 3233 with $A_{00}$ positive semi-definite. The implementation follows {cite}`arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$. 3234 A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`. 3235 3236 Developer Note: 3237 The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their 3238 user API. 3239 3240 .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`, 3241 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`, 3242 `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`, 3243 `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()` 3244 M*/ 3245 3246 PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc) 3247 { 3248 PC_FieldSplit *jac; 3249 3250 PetscFunctionBegin; 3251 PetscCall(PetscNew(&jac)); 3252 3253 jac->bs = -1; 3254 jac->nsplits = 0; 3255 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 3256 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 3257 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL; 3258 jac->schurscale = -1.0; 3259 jac->dm_splits = PETSC_TRUE; 3260 jac->detect = PETSC_FALSE; 3261 jac->gkbtol = 1e-5; 3262 jac->gkbdelay = 5; 3263 jac->gkbnu = 1; 3264 jac->gkbmaxit = 100; 3265 jac->gkbmonitor = PETSC_FALSE; 3266 jac->coordinates_set = PETSC_FALSE; 3267 3268 pc->data = (void *)jac; 3269 3270 pc->ops->apply = PCApply_FieldSplit; 3271 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 3272 pc->ops->setup = PCSetUp_FieldSplit; 3273 pc->ops->reset = PCReset_FieldSplit; 3274 pc->ops->destroy = PCDestroy_FieldSplit; 3275 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 3276 pc->ops->view = PCView_FieldSplit; 3277 pc->ops->applyrichardson = NULL; 3278 3279 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit)); 3280 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit)); 3281 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit)); 3282 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit)); 3283 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit)); 3284 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit)); 3285 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit)); 3286 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit)); 3287 PetscFunctionReturn(PETSC_SUCCESS); 3288 } 3289