1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 3 typedef struct { 4 PetscBool allocated, commute, scalediag; 5 KSP kspL, kspMass; 6 Vec Avec0, Avec1, Svec0, scale; 7 Mat L; 8 } PC_LSC; 9 10 static PetscErrorCode PCLSCAllocate_Private(PC pc) 11 { 12 PC_LSC *lsc = (PC_LSC *)pc->data; 13 Mat A; 14 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)); 28 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; 37 38 lsc->allocated = PETSC_TRUE; 39 PetscFunctionReturn(PETSC_SUCCESS); 40 } 41 42 static PetscErrorCode PCSetUp_LSC(PC pc) 43 { 44 PC_LSC *lsc = (PC_LSC *)pc->data; 45 Mat L, Lp, Qscale; 46 47 PetscFunctionBegin; 48 PetscCall(PCLSCAllocate_Private(pc)); 49 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)); 55 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)); 59 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; 64 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_DEFAULT, &lsc->L)); 80 else PetscCall(MatMatMult(CAdiaginv, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &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 } 94 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 } 103 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; 108 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 } 136 137 static PetscErrorCode PCReset_LSC(PC pc) 138 { 139 PC_LSC *lsc = (PC_LSC *)pc->data; 140 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 } 151 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 } 159 160 static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject) 161 { 162 PC_LSC *lsc = (PC_LSC *)pc->data; 163 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 } 174 175 static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer) 176 { 177 PC_LSC *jac = (PC_LSC *)pc->data; 178 PetscBool iascii; 179 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 } 200 201 /*MC 202 PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators 203 204 Options Database Key: 205 . -pc_lsc_scale_diag - Use the diagonal of A for scaling 206 207 Level: intermediate 208 209 Notes: 210 This preconditioner will normally be used with `PCFIELDSPLIT` to precondition the Schur complement, but 211 it can be used for any Schur complement system. Consider the Schur complement 212 213 .vb 214 S = A11 - A10 inv(A00) A01 215 .ve 216 217 `PCLSC` currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to 218 inv(S) is given by 219 220 .vb 221 inv(A10 A01) A10 A00 A01 inv(A10 A01) 222 .ve 223 224 The product A10 A01 can be computed for you, but you can provide it (this is 225 usually more efficient anyway). In the case of incompressible flow, A10 A01 is a Laplacian; call it L. The current 226 interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix. 227 228 If you had called `KSPSetOperators`(ksp,S,Sp), S should have type `MATSCHURCOMPLEMENT` and Sp can be any type you 229 like (`PCLSC` doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp". 230 For example, you might have setup code like this 231 232 .vb 233 PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L); 234 PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp); 235 .ve 236 237 And then your Jacobian assembly would look like 238 239 .vb 240 PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L); 241 PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp); 242 if (L) { assembly L } 243 if (Lp) { assemble Lp } 244 .ve 245 246 With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L 247 248 .vb 249 -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml 250 .ve 251 252 Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners. 253 254 References: 255 + * - Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006. 256 - * - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier Stokes equations for incompressible flow, 2001. 257 258 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`, 259 `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, 260 `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`, 261 `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()` 262 M*/ 263 264 PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc) 265 { 266 PC_LSC *lsc; 267 268 PetscFunctionBegin; 269 PetscCall(PetscNew(&lsc)); 270 pc->data = (void *)lsc; 271 272 pc->ops->apply = PCApply_LSC; 273 pc->ops->applytranspose = NULL; 274 pc->ops->setup = PCSetUp_LSC; 275 pc->ops->reset = PCReset_LSC; 276 pc->ops->destroy = PCDestroy_LSC; 277 pc->ops->setfromoptions = PCSetFromOptions_LSC; 278 pc->ops->view = PCView_LSC; 279 pc->ops->applyrichardson = NULL; 280 PetscFunctionReturn(PETSC_SUCCESS); 281 } 282