1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 3 typedef struct { 4 PetscBool allocated; 5 PetscBool scalediag; 6 KSP kspL; 7 Vec scale; 8 Vec x0, y0, x1; 9 Mat L; /* keep a copy to reuse when obtained with L = A10*A01 */ 10 } PC_LSC; 11 12 static PetscErrorCode PCLSCAllocate_Private(PC pc) { 13 PC_LSC *lsc = (PC_LSC *)pc->data; 14 Mat A; 15 16 PetscFunctionBegin; 17 if (lsc->allocated) PetscFunctionReturn(0); 18 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL)); 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->x0, &lsc->y0)); 26 PetscCall(MatCreateVecs(pc->pmat, &lsc->x1, NULL)); 27 if (lsc->scalediag) { PetscCall(VecDuplicate(lsc->x0, &lsc->scale)); } 28 lsc->allocated = PETSC_TRUE; 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode PCSetUp_LSC(PC pc) { 33 PC_LSC *lsc = (PC_LSC *)pc->data; 34 Mat L, Lp, B, C; 35 36 PetscFunctionBegin; 37 PetscCall(PCLSCAllocate_Private(pc)); 38 PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_L", (PetscObject *)&L)); 39 if (!L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_L", (PetscObject *)&L)); 40 PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Lp", (PetscObject *)&Lp)); 41 if (!Lp) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Lp", (PetscObject *)&Lp)); 42 if (!L) { 43 PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, NULL, &B, &C, NULL)); 44 if (!lsc->L) { 45 PetscCall(MatMatMult(C, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &lsc->L)); 46 } else { 47 PetscCall(MatMatMult(C, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &lsc->L)); 48 } 49 Lp = L = lsc->L; 50 } 51 if (lsc->scale) { 52 Mat Ap; 53 PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, &Ap, NULL, NULL, NULL)); 54 PetscCall(MatGetDiagonal(Ap, lsc->scale)); /* Should be the mass matrix, but we don't have plumbing for that yet */ 55 PetscCall(VecReciprocal(lsc->scale)); 56 } 57 PetscCall(KSPSetOperators(lsc->kspL, L, Lp)); 58 PetscCall(KSPSetFromOptions(lsc->kspL)); 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode PCApply_LSC(PC pc, Vec x, Vec y) { 63 PC_LSC *lsc = (PC_LSC *)pc->data; 64 Mat A, B, C; 65 66 PetscFunctionBegin; 67 PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL)); 68 PetscCall(KSPSolve(lsc->kspL, x, lsc->x1)); 69 PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->x1)); 70 PetscCall(MatMult(B, lsc->x1, lsc->x0)); 71 if (lsc->scale) PetscCall(VecPointwiseMult(lsc->x0, lsc->x0, lsc->scale)); 72 PetscCall(MatMult(A, lsc->x0, lsc->y0)); 73 if (lsc->scale) PetscCall(VecPointwiseMult(lsc->y0, lsc->y0, lsc->scale)); 74 PetscCall(MatMult(C, lsc->y0, lsc->x1)); 75 PetscCall(KSPSolve(lsc->kspL, lsc->x1, y)); 76 PetscCall(KSPCheckSolve(lsc->kspL, pc, y)); 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode PCReset_LSC(PC pc) { 81 PC_LSC *lsc = (PC_LSC *)pc->data; 82 83 PetscFunctionBegin; 84 PetscCall(VecDestroy(&lsc->x0)); 85 PetscCall(VecDestroy(&lsc->y0)); 86 PetscCall(VecDestroy(&lsc->x1)); 87 PetscCall(VecDestroy(&lsc->scale)); 88 PetscCall(KSPDestroy(&lsc->kspL)); 89 PetscCall(MatDestroy(&lsc->L)); 90 PetscFunctionReturn(0); 91 } 92 93 static PetscErrorCode PCDestroy_LSC(PC pc) { 94 PetscFunctionBegin; 95 PetscCall(PCReset_LSC(pc)); 96 PetscCall(PetscFree(pc->data)); 97 PetscFunctionReturn(0); 98 } 99 100 static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject) { 101 PC_LSC *lsc = (PC_LSC *)pc->data; 102 103 PetscFunctionBegin; 104 PetscOptionsHeadBegin(PetscOptionsObject, "LSC options"); 105 { PetscCall(PetscOptionsBool("-pc_lsc_scale_diag", "Use diagonal of velocity block (A) for scaling", "None", lsc->scalediag, &lsc->scalediag, NULL)); } 106 PetscOptionsHeadEnd(); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer) { 111 PC_LSC *jac = (PC_LSC *)pc->data; 112 PetscBool iascii; 113 114 PetscFunctionBegin; 115 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 116 if (iascii) { 117 PetscCall(PetscViewerASCIIPushTab(viewer)); 118 if (jac->kspL) { 119 PetscCall(KSPView(jac->kspL, viewer)); 120 } else { 121 PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display")); 122 } 123 PetscCall(PetscViewerASCIIPopTab(viewer)); 124 } 125 PetscFunctionReturn(0); 126 } 127 128 /*MC 129 PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators 130 131 Options Database Key: 132 . -pc_lsc_scale_diag - Use the diagonal of A for scaling 133 134 Level: intermediate 135 136 Notes: 137 This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but 138 it can be used for any Schur complement system. Consider the Schur complement 139 140 .vb 141 S = A11 - A10 inv(A00) A01 142 .ve 143 144 PCLSC currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to 145 inv(S) is given by 146 147 .vb 148 inv(A10 A01) A10 A00 A01 inv(A10 A01) 149 .ve 150 151 The product A10 A01 can be computed for you, but you can provide it (this is 152 usually more efficient anyway). In the case of incompressible flow, A10 A01 is a Laplacian; call it L. The current 153 interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix. 154 155 If you had called KSPSetOperators(ksp,S,Sp), S should have type MATSCHURCOMPLEMENT and Sp can be any type you 156 like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp". 157 For example, you might have setup code like this 158 159 .vb 160 PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L); 161 PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp); 162 .ve 163 164 And then your Jacobian assembly would look like 165 166 .vb 167 PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L); 168 PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp); 169 if (L) { assembly L } 170 if (Lp) { assemble Lp } 171 .ve 172 173 With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L 174 175 .vb 176 -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml 177 .ve 178 179 Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners. 180 181 References: 182 + * - Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006. 183 - * - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier Stokes equations for incompressible flow, 2001. 184 185 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`, 186 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, 187 `MatCreateSchurComplement()` 188 M*/ 189 190 PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc) { 191 PC_LSC *lsc; 192 193 PetscFunctionBegin; 194 PetscCall(PetscNewLog(pc, &lsc)); 195 pc->data = (void *)lsc; 196 197 pc->ops->apply = PCApply_LSC; 198 pc->ops->applytranspose = NULL; 199 pc->ops->setup = PCSetUp_LSC; 200 pc->ops->reset = PCReset_LSC; 201 pc->ops->destroy = PCDestroy_LSC; 202 pc->ops->setfromoptions = PCSetFromOptions_LSC; 203 pc->ops->view = PCView_LSC; 204 pc->ops->applyrichardson = NULL; 205 PetscFunctionReturn(0); 206 } 207