Actual source code: lsc.c
petsc-3.6.4 2016-04-12
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: KSPSetErrorIfNotConverged(lsc->kspL,pc->erroriffailure);
25: PetscObjectIncrementTabLevel((PetscObject)lsc->kspL,(PetscObject)pc,1);
26: KSPSetType(lsc->kspL,KSPPREONLY);
27: KSPSetOptionsPrefix(lsc->kspL,((PetscObject)pc)->prefix);
28: KSPAppendOptionsPrefix(lsc->kspL,"lsc_");
29: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,NULL,NULL,NULL);
30: MatCreateVecs(A,&lsc->x0,&lsc->y0);
31: MatCreateVecs(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: KSPSetFromOptions(lsc->kspL);
70: return(0);
71: }
75: static PetscErrorCode PCApply_LSC(PC pc,Vec x,Vec y)
76: {
77: PC_LSC *lsc = (PC_LSC*)pc->data;
78: Mat A,B,C;
82: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,&B,&C,NULL);
83: KSPSolve(lsc->kspL,x,lsc->x1);
84: MatMult(B,lsc->x1,lsc->x0);
85: if (lsc->scale) {
86: VecPointwiseMult(lsc->x0,lsc->x0,lsc->scale);
87: }
88: MatMult(A,lsc->x0,lsc->y0);
89: if (lsc->scale) {
90: VecPointwiseMult(lsc->y0,lsc->y0,lsc->scale);
91: }
92: MatMult(C,lsc->y0,lsc->x1);
93: KSPSolve(lsc->kspL,lsc->x1,y);
94: return(0);
95: }
99: static PetscErrorCode PCReset_LSC(PC pc)
100: {
101: PC_LSC *lsc = (PC_LSC*)pc->data;
105: VecDestroy(&lsc->x0);
106: VecDestroy(&lsc->y0);
107: VecDestroy(&lsc->x1);
108: VecDestroy(&lsc->scale);
109: KSPDestroy(&lsc->kspL);
110: MatDestroy(&lsc->L);
111: return(0);
112: }
116: static PetscErrorCode PCDestroy_LSC(PC pc)
117: {
121: PCReset_LSC(pc);
122: PetscFree(pc->data);
123: return(0);
124: }
128: static PetscErrorCode PCSetFromOptions_LSC(PetscOptions *PetscOptionsObject,PC pc)
129: {
130: PC_LSC *lsc = (PC_LSC*)pc->data;
134: PetscOptionsHead(PetscOptionsObject,"LSC options");
135: {
136: PetscOptionsBool("-pc_lsc_scale_diag","Use diagonal of velocity block (A) for scaling","None",lsc->scalediag,&lsc->scalediag,NULL);
137: }
138: PetscOptionsTail();
139: return(0);
140: }
144: static PetscErrorCode PCView_LSC(PC pc,PetscViewer viewer)
145: {
146: PC_LSC *jac = (PC_LSC*)pc->data;
148: PetscBool iascii;
151: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
152: if (iascii) {
153: PetscViewerASCIIPushTab(viewer);
154: KSPView(jac->kspL,viewer);
155: PetscViewerASCIIPopTab(viewer);
156: }
157: return(0);
158: }
160: /*MC
161: PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators
163: Options Database Key:
164: . -pc_lsc_scale_diag - Use the diagonal of A for scaling
166: Level: intermediate
168: Notes:
169: This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but
170: it can be used for any Schur complement system. Consider the Schur complement
172: .vb
173: S = A11 - A10 inv(A00) A01
174: .ve
176: PCLSC currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to
177: inv(S) is given by
179: .vb
180: inv(A10 A01) A10 A00 A01 inv(A10 A01)
181: .ve
183: The product A10 A01 can be computed for you, but you can provide it (this is
184: usually more efficient anyway). In the case of incompressible flow, A10 A10 is a Laplacian, call it L. The current
185: interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix.
187: If you had called KSPSetOperators(ksp,S,Sp), S should have type MATSCHURCOMPLEMENT and Sp can be any type you
188: like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
189: For example, you might have setup code like this
191: .vb
192: PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
193: PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
194: .ve
196: And then your Jacobian assembly would look like
198: .vb
199: PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
200: PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
201: if (L) { assembly L }
202: if (Lp) { assemble Lp }
203: .ve
205: With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L
207: .vb
208: -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
209: .ve
211: Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.
213: References:
214: + Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006.
215: - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow, 2001.
217: Concepts: physics based preconditioners, block preconditioners
219: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCFIELDSPLIT,
220: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
221: MatCreateSchurComplement()
222: M*/
226: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
227: {
228: PC_LSC *lsc;
232: PetscNewLog(pc,&lsc);
233: pc->data = (void*)lsc;
235: pc->ops->apply = PCApply_LSC;
236: pc->ops->applytranspose = 0;
237: pc->ops->setup = PCSetUp_LSC;
238: pc->ops->reset = PCReset_LSC;
239: pc->ops->destroy = PCDestroy_LSC;
240: pc->ops->setfromoptions = PCSetFromOptions_LSC;
241: pc->ops->view = PCView_LSC;
242: pc->ops->applyrichardson = 0;
243: return(0);
244: }