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