Actual source code: galerkin.c

      1: /*
      2:       Defines a preconditioner defined by R^T S R
      3: */
      4: #include <petsc/private/pcimpl.h>
      5: #include <petscksp.h>

      7: typedef struct {
      8:   KSP ksp;
      9:   Mat R, P;
     10:   Vec b, x;
     11:   PetscErrorCode (*computeasub)(PC, Mat, Mat, Mat *, void *);
     12:   void *computeasub_ctx;
     13: } PC_Galerkin;

     15: static PetscErrorCode PCApply_Galerkin(PC pc, Vec x, Vec y)
     16: {
     17:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

     19:   PetscFunctionBegin;
     20:   if (jac->R) {
     21:     PetscCall(MatRestrict(jac->R, x, jac->b));
     22:   } else {
     23:     PetscCall(MatRestrict(jac->P, x, jac->b));
     24:   }
     25:   PetscCall(KSPSolve(jac->ksp, jac->b, jac->x));
     26:   PetscCall(KSPCheckSolve(jac->ksp, pc, jac->x));
     27:   if (jac->P) {
     28:     PetscCall(MatInterpolate(jac->P, jac->x, y));
     29:   } else {
     30:     PetscCall(MatInterpolate(jac->R, jac->x, y));
     31:   }
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: static PetscErrorCode PCSetUp_Galerkin(PC pc)
     36: {
     37:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;
     38:   PetscBool    a;
     39:   Vec         *xx, *yy;

     41:   PetscFunctionBegin;
     42:   if (jac->computeasub) {
     43:     Mat Ap;
     44:     if (!pc->setupcalled) {
     45:       PetscCall((*jac->computeasub)(pc, pc->pmat, NULL, &Ap, jac->computeasub_ctx));
     46:       PetscCall(KSPSetOperators(jac->ksp, Ap, Ap));
     47:       PetscCall(MatDestroy(&Ap));
     48:     } else {
     49:       PetscCall(KSPGetOperators(jac->ksp, NULL, &Ap));
     50:       PetscCall((*jac->computeasub)(pc, pc->pmat, Ap, NULL, jac->computeasub_ctx));
     51:     }
     52:   }

     54:   if (!jac->x) {
     55:     PetscCall(KSPGetOperatorsSet(jac->ksp, &a, NULL));
     56:     PetscCheck(a, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
     57:     PetscCall(KSPCreateVecs(jac->ksp, 1, &xx, 1, &yy));
     58:     jac->x = *xx;
     59:     jac->b = *yy;
     60:     PetscCall(PetscFree(xx));
     61:     PetscCall(PetscFree(yy));
     62:   }
     63:   PetscCheck(jac->R || jac->P, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()");
     64:   /* should check here that sizes of R/P match size of a */
     65:   PetscFunctionReturn(PETSC_SUCCESS);
     66: }

     68: static PetscErrorCode PCReset_Galerkin(PC pc)
     69: {
     70:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

     72:   PetscFunctionBegin;
     73:   PetscCall(MatDestroy(&jac->R));
     74:   PetscCall(MatDestroy(&jac->P));
     75:   PetscCall(VecDestroy(&jac->x));
     76:   PetscCall(VecDestroy(&jac->b));
     77:   PetscCall(KSPReset(jac->ksp));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: static PetscErrorCode PCDestroy_Galerkin(PC pc)
     82: {
     83:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

     85:   PetscFunctionBegin;
     86:   PetscCall(PCReset_Galerkin(pc));
     87:   PetscCall(KSPDestroy(&jac->ksp));
     88:   PetscCall(PetscFree(pc->data));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
     93: {
     94:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;
     95:   PetscBool    iascii;

     97:   PetscFunctionBegin;
     98:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     99:   if (iascii) {
    100:     PetscCall(PetscViewerASCIIPrintf(viewer, "  KSP on Galerkin follow\n"));
    101:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
    102:   }
    103:   PetscCall(KSPView(jac->ksp, viewer));
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
    108: {
    109:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

    111:   PetscFunctionBegin;
    112:   *ksp = jac->ksp;
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
    117: {
    118:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

    120:   PetscFunctionBegin;
    121:   PetscCall(PetscObjectReference((PetscObject)R));
    122:   PetscCall(MatDestroy(&jac->R));
    123:   jac->R = R;
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
    128: {
    129:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

    131:   PetscFunctionBegin;
    132:   PetscCall(PetscObjectReference((PetscObject)P));
    133:   PetscCall(MatDestroy(&jac->P));
    134:   jac->P = P;
    135:   PetscFunctionReturn(PETSC_SUCCESS);
    136: }

    138: static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
    139: {
    140:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

    142:   PetscFunctionBegin;
    143:   jac->computeasub     = computeAsub;
    144:   jac->computeasub_ctx = ctx;
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: /*@
    149:   PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner

    151:   Logically Collective

    153:   Input Parameters:
    154: + pc - the preconditioner context
    155: - R  - the restriction operator

    157:   Level: intermediate

    159:   Note:
    160:   Either this or `PCGalerkinSetInterpolation()` or both must be called

    162: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
    163:           `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
    164: @*/
    165: PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
    166: {
    167:   PetscFunctionBegin;
    169:   PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*@
    174:   PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner

    176:   Logically Collective

    178:   Input Parameters:
    179: + pc - the preconditioner context
    180: - P  - the interpolation operator

    182:   Level: intermediate

    184:   Note:
    185:   Either this or `PCGalerkinSetRestriction()` or both must be called

    187: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
    188:           `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
    189: @*/
    190: PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
    191: {
    192:   PetscFunctionBegin;
    194:   PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
    195:   PetscFunctionReturn(PETSC_SUCCESS);
    196: }

    198: /*@C
    199:   PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix

    201:   Logically Collective

    203:   Input Parameters:
    204: + pc          - the preconditioner context
    205: . computeAsub - routine that computes the submatrix from the global matrix
    206: - ctx         - context used by the routine, or `NULL`

    208:   Calling sequence of `computeAsub`:
    209: + pc  - the `PCGALERKIN` preconditioner
    210: . A   - the matrix in the `PCGALERKIN`
    211: . Ap  - the computed submatrix from any previous computation, if `NULL` it has not previously been computed
    212: . cAp - the submatrix computed by this routine
    213: - ctx - optional user-defined function context

    215:   Level: intermediate

    217:   Notes:
    218:   Instead of providing this routine you can call `PCGalerkinGetKSP()` and then `KSPSetOperators()` to provide the submatrix,
    219:   but that will not work for multiple `KSPSolve()`s with different matrices unless you call it for each solve.

    221:   This routine is called each time the outer matrix is changed. In the first call the Ap argument is `NULL` and the routine should create the
    222:   matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.

    224:   Developer Notes:
    225:   If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN`
    226:   could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()`

    228: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
    229:           `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
    230: @*/
    231: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC pc, Mat A, Mat Ap, Mat *cAp, void *ctx), void *ctx)
    232: {
    233:   PetscFunctionBegin;
    235:   PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode (*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /*@
    240:   PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`

    242:   Not Collective

    244:   Input Parameter:
    245: . pc - the preconditioner context

    247:   Output Parameter:
    248: . ksp - the `KSP` object

    250:   Level: intermediate

    252:   Note:
    253:   Once you have called this routine you can call `KSPSetOperators()` on the resulting `KSP` to provide the operator for the Galerkin problem,
    254:   an alternative is to use `PCGalerkinSetComputeSubmatrix()` to provide a routine that computes the submatrix as needed.

    256: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
    257:           `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
    258: @*/
    259: PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
    260: {
    261:   PetscFunctionBegin;
    263:   PetscAssertPointer(ksp, 2);
    264:   PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems *PetscOptionsObject)
    269: {
    270:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;
    271:   const char  *prefix;
    272:   PetscBool    flg;

    274:   PetscFunctionBegin;
    275:   PetscCall(KSPGetOptionsPrefix(jac->ksp, &prefix));
    276:   PetscCall(PetscStrendswith(prefix, "galerkin_", &flg));
    277:   if (!flg) {
    278:     PetscCall(PCGetOptionsPrefix(pc, &prefix));
    279:     PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
    280:     PetscCall(KSPAppendOptionsPrefix(jac->ksp, "galerkin_"));
    281:   }

    283:   PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
    284:   if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
    285:   PetscOptionsHeadEnd();
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: /*MC
    290:      PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)

    292:     Level: intermediate

    294:     Note:
    295:     Use
    296: .vb
    297:      `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P)
    298:      `PCGalerkinGetKSP`(pc,&ksp);
    299:      `KSPSetOperators`(ksp,A,....)
    300:      ...
    301: .ve

    303:     Developer Notes:
    304:     If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute
    305:     the operators automatically.

    307:     Should there be a prefix for the inner `KSP`?

    309:     There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`

    311: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
    312:           `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
    313: M*/

    315: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
    316: {
    317:   PC_Galerkin *jac;

    319:   PetscFunctionBegin;
    320:   PetscCall(PetscNew(&jac));

    322:   pc->ops->apply           = PCApply_Galerkin;
    323:   pc->ops->setup           = PCSetUp_Galerkin;
    324:   pc->ops->reset           = PCReset_Galerkin;
    325:   pc->ops->destroy         = PCDestroy_Galerkin;
    326:   pc->ops->view            = PCView_Galerkin;
    327:   pc->ops->setfromoptions  = PCSetFromOptions_Galerkin;
    328:   pc->ops->applyrichardson = NULL;

    330:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
    331:   PetscCall(KSPSetNestLevel(jac->ksp, pc->kspnestlevel));
    332:   PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
    333:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));

    335:   pc->data = (void *)jac;

    337:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin));
    338:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin));
    339:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin));
    340:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin));
    341:   PetscFunctionReturn(PETSC_SUCCESS);
    342: }