Actual source code: lsc.c

      1: #include <petsc/private/pcimpl.h>

      3: typedef struct {
      4:   PetscBool allocated, commute, scalediag;
      5:   KSP       kspL, kspMass;
      6:   Vec       Avec0, Avec1, Svec0, scale;
      7:   Mat       L;
      8: } PC_LSC;

     10: static PetscErrorCode PCLSCAllocate_Private(PC pc)
     11: {
     12:   PC_LSC *lsc = (PC_LSC *)pc->data;
     13:   Mat     A;

     15:   PetscFunctionBegin;
     16:   if (lsc->allocated) PetscFunctionReturn(PETSC_SUCCESS);
     17:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL));
     18:   PetscCall(KSPSetNestLevel(lsc->kspL, pc->kspnestlevel));
     19:   PetscCall(KSPSetErrorIfNotConverged(lsc->kspL, pc->erroriffailure));
     20:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspL, (PetscObject)pc, 1));
     21:   PetscCall(KSPSetType(lsc->kspL, KSPPREONLY));
     22:   PetscCall(KSPSetOptionsPrefix(lsc->kspL, ((PetscObject)pc)->prefix));
     23:   PetscCall(KSPAppendOptionsPrefix(lsc->kspL, "lsc_"));
     24:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, NULL, NULL, NULL));
     25:   PetscCall(MatCreateVecs(A, &lsc->Avec0, &lsc->Avec1));
     26:   PetscCall(MatCreateVecs(pc->pmat, &lsc->Svec0, NULL));
     27:   if (lsc->scalediag) PetscCall(VecDuplicate(lsc->Avec0, &lsc->scale));

     29:   if (lsc->commute) {
     30:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspMass));
     31:     PetscCall(KSPSetErrorIfNotConverged(lsc->kspMass, pc->erroriffailure));
     32:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspMass, (PetscObject)pc, 1));
     33:     PetscCall(KSPSetType(lsc->kspMass, KSPPREONLY));
     34:     PetscCall(KSPSetOptionsPrefix(lsc->kspMass, ((PetscObject)pc)->prefix));
     35:     PetscCall(KSPAppendOptionsPrefix(lsc->kspMass, "lsc_mass_"));
     36:   } else lsc->kspMass = NULL;

     38:   lsc->allocated = PETSC_TRUE;
     39:   PetscFunctionReturn(PETSC_SUCCESS);
     40: }

     42: static PetscErrorCode PCSetUp_LSC(PC pc)
     43: {
     44:   PC_LSC *lsc = (PC_LSC *)pc->data;
     45:   Mat     L, Lp, Qscale;

     47:   PetscFunctionBegin;
     48:   PetscCall(PCLSCAllocate_Private(pc));

     50:   /* Query for L operator */
     51:   PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_L", (PetscObject *)&L));
     52:   if (!L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_L", (PetscObject *)&L));
     53:   PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Lp", (PetscObject *)&Lp));
     54:   if (!Lp) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Lp", (PetscObject *)&Lp));

     56:   /* Query for mass operator */
     57:   PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Qscale", (PetscObject *)&Qscale));
     58:   if (!Qscale) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Qscale", (PetscObject *)&Qscale));

     60:   if (lsc->commute) {
     61:     PetscCheck(L || Lp, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide an L operator for LSC preconditioning when commuting");
     62:     if (!L && Lp) L = Lp;
     63:     else if (L && !Lp) Lp = L;

     65:     PetscCheck(Qscale, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide a Qscale matrix for LSC preconditioning when commuting");
     66:   } else {
     67:     if (lsc->scale) {
     68:       if (!Qscale) PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, &Qscale, NULL, NULL, NULL));
     69:       PetscCall(MatGetDiagonal(Qscale, lsc->scale));
     70:       PetscCall(VecReciprocal(lsc->scale));
     71:     }
     72:     if (!L) {
     73:       Mat B, C;
     74:       PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, NULL, &B, &C, NULL));
     75:       if (lsc->scale) {
     76:         Mat CAdiaginv;
     77:         PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &CAdiaginv));
     78:         PetscCall(MatDiagonalScale(CAdiaginv, NULL, lsc->scale));
     79:         if (!lsc->L) PetscCall(MatMatMult(CAdiaginv, B, MAT_INITIAL_MATRIX, PETSC_CURRENT, &lsc->L));
     80:         else PetscCall(MatMatMult(CAdiaginv, B, MAT_REUSE_MATRIX, PETSC_CURRENT, &lsc->L));
     81:         PetscCall(MatDestroy(&CAdiaginv));
     82:       } else {
     83:         if (!lsc->L) {
     84:           PetscCall(MatProductCreate(C, B, NULL, &lsc->L));
     85:           PetscCall(MatProductSetType(lsc->L, MATPRODUCT_AB));
     86:           PetscCall(MatProductSetFromOptions(lsc->L));
     87:           PetscCall(MatProductSymbolic(lsc->L));
     88:         }
     89:         PetscCall(MatProductNumeric(lsc->L));
     90:       }
     91:       Lp = L = lsc->L;
     92:     }
     93:   }

     95:   PetscCall(KSPSetOperators(lsc->kspL, L, Lp));
     96:   PetscCall(KSPSetFromOptions(lsc->kspL));
     97:   if (lsc->commute) {
     98:     PetscCall(KSPSetOperators(lsc->kspMass, Qscale, Qscale));
     99:     PetscCall(KSPSetFromOptions(lsc->kspMass));
    100:   }
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode PCApply_LSC(PC pc, Vec x, Vec y)
    105: {
    106:   PC_LSC *lsc = (PC_LSC *)pc->data;
    107:   Mat     A, B, C;

    109:   PetscFunctionBegin;
    110:   PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL));
    111:   if (lsc->commute) {
    112:     PetscCall(KSPSolve(lsc->kspMass, x, lsc->Svec0));
    113:     PetscCall(KSPCheckSolve(lsc->kspMass, pc, lsc->Svec0));
    114:     PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
    115:     PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
    116:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
    117:     PetscCall(MatMult(A, lsc->Avec1, lsc->Avec0));
    118:     PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
    119:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
    120:     PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
    121:     PetscCall(KSPSolve(lsc->kspMass, lsc->Svec0, y));
    122:     PetscCall(KSPCheckSolve(lsc->kspMass, pc, y));
    123:   } else {
    124:     PetscCall(KSPSolve(lsc->kspL, x, lsc->Svec0));
    125:     PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Svec0));
    126:     PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
    127:     if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec0, lsc->Avec0, lsc->scale));
    128:     PetscCall(MatMult(A, lsc->Avec0, lsc->Avec1));
    129:     if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec1, lsc->Avec1, lsc->scale));
    130:     PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
    131:     PetscCall(KSPSolve(lsc->kspL, lsc->Svec0, y));
    132:     PetscCall(KSPCheckSolve(lsc->kspL, pc, y));
    133:   }
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: static PetscErrorCode PCReset_LSC(PC pc)
    138: {
    139:   PC_LSC *lsc = (PC_LSC *)pc->data;

    141:   PetscFunctionBegin;
    142:   PetscCall(VecDestroy(&lsc->Avec0));
    143:   PetscCall(VecDestroy(&lsc->Avec1));
    144:   PetscCall(VecDestroy(&lsc->Svec0));
    145:   PetscCall(KSPDestroy(&lsc->kspL));
    146:   if (lsc->commute) PetscCall(KSPDestroy(&lsc->kspMass));
    147:   if (lsc->L) PetscCall(MatDestroy(&lsc->L));
    148:   if (lsc->scale) PetscCall(VecDestroy(&lsc->scale));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: static PetscErrorCode PCDestroy_LSC(PC pc)
    153: {
    154:   PetscFunctionBegin;
    155:   PetscCall(PCReset_LSC(pc));
    156:   PetscCall(PetscFree(pc->data));
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

    160: static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject)
    161: {
    162:   PC_LSC *lsc = (PC_LSC *)pc->data;

    164:   PetscFunctionBegin;
    165:   PetscOptionsHeadBegin(PetscOptionsObject, "LSC options");
    166:   {
    167:     PetscCall(PetscOptionsBool("-pc_lsc_commute", "Whether to commute the LSC preconditioner in the style of Olshanskii", "None", lsc->commute, &lsc->commute, NULL));
    168:     PetscCall(PetscOptionsBool("-pc_lsc_scale_diag", "Whether to scale BBt products. Will use the inverse of the diagonal of Qscale or A if the former is not provided.", "None", lsc->scalediag, &lsc->scalediag, NULL));
    169:     PetscCheck(!lsc->scalediag || !lsc->commute, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Diagonal-based scaling is not used when doing a commuted LSC. Either do not ask for diagonal-based scaling or use non-commuted LSC in the original style of Elman");
    170:   }
    171:   PetscOptionsHeadEnd();
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }

    175: static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
    176: {
    177:   PC_LSC   *jac = (PC_LSC *)pc->data;
    178:   PetscBool iascii;

    180:   PetscFunctionBegin;
    181:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    182:   if (iascii) {
    183:     PetscCall(PetscViewerASCIIPushTab(viewer));
    184:     if (jac->kspL) {
    185:       PetscCall(KSPView(jac->kspL, viewer));
    186:     } else {
    187:       PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display"));
    188:     }
    189:     if (jac->commute) {
    190:       if (jac->kspMass) {
    191:         PetscCall(KSPView(jac->kspMass, viewer));
    192:       } else {
    193:         PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC Mass KSP object not yet created, hence cannot display"));
    194:       }
    195:     }
    196:     PetscCall(PetscViewerASCIIPopTab(viewer));
    197:   }
    198:   PetscFunctionReturn(PETSC_SUCCESS);
    199: }

    201: /*MC
    202:    PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators {cite}`elmanhowleshadidshuttleworthtuminaro2006` {cite}`silvester2001efficient`

    204:    Options Database Key:
    205: +    -pc_lsc_commute    - Whether to commute the LSC preconditioner in the style of Olshanskii
    206: -    -pc_lsc_scale_diag - Whether to scale $BB^T$ products. Will use the inverse of the diagonal of $Qscale$ or $A$ if the former is not provided

    208:    Level: intermediate

    210:    Notes:
    211:    This preconditioner will normally be used with `PCFIELDSPLIT` to precondition the Schur complement, but
    212:    it can be used for any Schur complement system.  Consider the Schur complement

    214:    $$
    215:    S = A11 - A10 A00^{-1} A01
    216:    $$

    218:    `PCLSC` currently doesn't do anything with $A11$, so let's assume it is 0.  The idea is that a good approximation to
    219:    $S^{-1}$ is given by

    221:    $$
    222:    (A10 A01)^{-1} A10 A00 A01 (A10 A01)^{-1}
    223:    $$

    225:    The product $A10 A01$ can be computed for you, but you can provide it (this is
    226:    usually more efficient anyway).  In the case of incompressible flow, $A10 A01$ is a Laplacian; call it $L$.  The current
    227:    interface is to compose $L$ and a matrix from which to construct its preconditioning  $Lp$ on the preconditioning matrix.

    229:    If you had called `KSPSetOperators`(ksp,S,Sp), $S$ should have type `MATSCHURCOMPLEMENT` and $Sp$ can be any type you
    230:    like (`PCLSC` doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
    231:    For example, you might have setup code like this

    233: .vb
    234:    PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
    235:    PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
    236: .ve

    238:    And then your Jacobian assembly would look like

    240: .vb
    241:    PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
    242:    PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
    243:    if (L) { assembly L }
    244:    if (Lp) { assemble Lp }
    245: .ve

    247:    With this, you should be able to choose LSC preconditioning, using e.g. the `PCML` algebraic multigrid to solve with L
    248: .vb
    249:    -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
    250: .ve

    252:    Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.

    254: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`,
    255:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`,
    256:           `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`,
    257:           `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()`
    258: M*/

    260: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
    261: {
    262:   PC_LSC *lsc;

    264:   PetscFunctionBegin;
    265:   PetscCall(PetscNew(&lsc));
    266:   pc->data = (void *)lsc;

    268:   pc->ops->apply           = PCApply_LSC;
    269:   pc->ops->applytranspose  = NULL;
    270:   pc->ops->setup           = PCSetUp_LSC;
    271:   pc->ops->reset           = PCReset_LSC;
    272:   pc->ops->destroy         = PCDestroy_LSC;
    273:   pc->ops->setfromoptions  = PCSetFromOptions_LSC;
    274:   pc->ops->view            = PCView_LSC;
    275:   pc->ops->applyrichardson = NULL;
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }