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