Actual source code: lsc.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
4: typedef struct {
5: PetscBool allocated;
6: PetscBool scalediag;
7: KSP kspL;
8: Vec scale;
9: Vec x0,y0,x1;
10: Mat L; /* keep a copy to reuse when obtained with L = A10*A01 */
11: } PC_LSC;
15: static PetscErrorCode PCLSCAllocate_Private(PC pc)
16: {
17: PC_LSC *lsc = (PC_LSC*)pc->data;
18: Mat A;
22: if (lsc->allocated) return(0);
23: KSPCreate(PetscObjectComm((PetscObject)pc),&lsc->kspL);
24: PetscObjectIncrementTabLevel((PetscObject)lsc->kspL,(PetscObject)pc,1);
25: KSPSetType(lsc->kspL,KSPPREONLY);
26: KSPSetOptionsPrefix(lsc->kspL,((PetscObject)pc)->prefix);
27: KSPAppendOptionsPrefix(lsc->kspL,"lsc_");
28: KSPSetFromOptions(lsc->kspL);
29: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,NULL,NULL,NULL);
30: MatGetVecs(A,&lsc->x0,&lsc->y0);
31: MatGetVecs(pc->pmat,&lsc->x1,NULL);
32: if (lsc->scalediag) {
33: VecDuplicate(lsc->x0,&lsc->scale);
34: }
35: lsc->allocated = PETSC_TRUE;
36: return(0);
37: }
41: static PetscErrorCode PCSetUp_LSC(PC pc)
42: {
43: PC_LSC *lsc = (PC_LSC*)pc->data;
44: Mat L,Lp,B,C;
48: PCLSCAllocate_Private(pc);
49: PetscObjectQuery((PetscObject)pc->mat,"LSC_L",(PetscObject*)&L);
50: if (!L) {PetscObjectQuery((PetscObject)pc->pmat,"LSC_L",(PetscObject*)&L);}
51: PetscObjectQuery((PetscObject)pc->pmat,"LSC_Lp",(PetscObject*)&Lp);
52: if (!Lp) {PetscObjectQuery((PetscObject)pc->mat,"LSC_Lp",(PetscObject*)&Lp);}
53: if (!L) {
54: MatSchurComplementGetSubMatrices(pc->mat,NULL,NULL,&B,&C,NULL);
55: if (!lsc->L) {
56: MatMatMult(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lsc->L);
57: } else {
58: MatMatMult(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&lsc->L);
59: }
60: Lp = L = lsc->L;
61: }
62: if (lsc->scale) {
63: Mat Ap;
64: MatSchurComplementGetSubMatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);
65: MatGetDiagonal(Ap,lsc->scale); /* Should be the mass matrix, but we don't have plumbing for that yet */
66: VecReciprocal(lsc->scale);
67: }
68: KSPSetOperators(lsc->kspL,L,Lp);
69: return(0);
70: }
74: static PetscErrorCode PCApply_LSC(PC pc,Vec x,Vec y)
75: {
76: PC_LSC *lsc = (PC_LSC*)pc->data;
77: Mat A,B,C;
81: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,&B,&C,NULL);
82: KSPSolve(lsc->kspL,x,lsc->x1);
83: MatMult(B,lsc->x1,lsc->x0);
84: if (lsc->scale) {
85: VecPointwiseMult(lsc->x0,lsc->x0,lsc->scale);
86: }
87: MatMult(A,lsc->x0,lsc->y0);
88: if (lsc->scale) {
89: VecPointwiseMult(lsc->y0,lsc->y0,lsc->scale);
90: }
91: MatMult(C,lsc->y0,lsc->x1);
92: KSPSolve(lsc->kspL,lsc->x1,y);
93: return(0);
94: }
98: static PetscErrorCode PCReset_LSC(PC pc)
99: {
100: PC_LSC *lsc = (PC_LSC*)pc->data;
104: VecDestroy(&lsc->x0);
105: VecDestroy(&lsc->y0);
106: VecDestroy(&lsc->x1);
107: VecDestroy(&lsc->scale);
108: KSPDestroy(&lsc->kspL);
109: MatDestroy(&lsc->L);
110: return(0);
111: }
115: static PetscErrorCode PCDestroy_LSC(PC pc)
116: {
120: PCReset_LSC(pc);
121: PetscFree(pc->data);
122: return(0);
123: }
127: static PetscErrorCode PCSetFromOptions_LSC(PC pc)
128: {
129: PC_LSC *lsc = (PC_LSC*)pc->data;
133: PetscOptionsHead("LSC options");
134: {
135: PetscOptionsBool("-pc_lsc_scale_diag","Use diagonal of velocity block (A) for scaling","None",lsc->scalediag,&lsc->scalediag,NULL);
136: }
137: PetscOptionsTail();
138: return(0);
139: }
143: static PetscErrorCode PCView_LSC(PC pc,PetscViewer viewer)
144: {
145: PC_LSC *jac = (PC_LSC*)pc->data;
147: PetscBool iascii;
150: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
151: if (iascii) {
152: PetscViewerASCIIPushTab(viewer);
153: KSPView(jac->kspL,viewer);
154: PetscViewerASCIIPopTab(viewer);
155: }
156: return(0);
157: }
159: /*MC
160: PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators
162: Options Database Key:
163: . -pc_lsc_scale_diag - Use the diagonal of A for scaling
165: Level: intermediate
167: Notes:
168: This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but
169: it can be used for any Schur complement system. Consider the Schur complement
171: .vb
172: S = A11 - A10 inv(A00) A01
173: .ve
175: PCLSC currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to
176: inv(S) is given by
178: .vb
179: inv(A10 A01) A10 A00 A01 inv(A10 A01)
180: .ve
182: The product A10 A01 can be computed for you, but you can provide it (this is
183: usually more efficient anyway). In the case of incompressible flow, A10 A10 is a Laplacian, call it L. The current
184: interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix.
186: If you had called KSPSetOperators(ksp,S,Sp), S should have type MATSCHURCOMPLEMENT and Sp can be any type you
187: like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
188: For example, you might have setup code like this
190: .vb
191: PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
192: PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
193: .ve
195: And then your Jacobian assembly would look like
197: .vb
198: PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
199: PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
200: if (L) { assembly L }
201: if (Lp) { assemble Lp }
202: .ve
204: With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L
206: .vb
207: -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
208: .ve
210: Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.
212: References:
213: + Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006.
214: - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow, 2001.
216: Concepts: physics based preconditioners, block preconditioners
218: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCFIELDSPLIT,
219: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
220: MatCreateSchurComplement()
221: M*/
225: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
226: {
227: PC_LSC *lsc;
231: PetscNewLog(pc,&lsc);
232: pc->data = (void*)lsc;
234: pc->ops->apply = PCApply_LSC;
235: pc->ops->applytranspose = 0;
236: pc->ops->setup = PCSetUp_LSC;
237: pc->ops->reset = PCReset_LSC;
238: pc->ops->destroy = PCDestroy_LSC;
239: pc->ops->setfromoptions = PCSetFromOptions_LSC;
240: pc->ops->view = PCView_LSC;
241: pc->ops->applyrichardson = 0;
242: return(0);
243: }