Actual source code: lsc.c
petsc-3.12.5 2020-03-29
2: #include <petsc/private/pcimpl.h>
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;
13: static PetscErrorCode PCLSCAllocate_Private(PC pc)
14: {
15: PC_LSC *lsc = (PC_LSC*)pc->data;
16: Mat A;
20: if (lsc->allocated) return(0);
21: KSPCreate(PetscObjectComm((PetscObject)pc),&lsc->kspL);
22: KSPSetErrorIfNotConverged(lsc->kspL,pc->erroriffailure);
23: PetscObjectIncrementTabLevel((PetscObject)lsc->kspL,(PetscObject)pc,1);
24: KSPSetType(lsc->kspL,KSPPREONLY);
25: KSPSetOptionsPrefix(lsc->kspL,((PetscObject)pc)->prefix);
26: KSPAppendOptionsPrefix(lsc->kspL,"lsc_");
27: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,NULL,NULL,NULL);
28: MatCreateVecs(A,&lsc->x0,&lsc->y0);
29: MatCreateVecs(pc->pmat,&lsc->x1,NULL);
30: if (lsc->scalediag) {
31: VecDuplicate(lsc->x0,&lsc->scale);
32: }
33: lsc->allocated = PETSC_TRUE;
34: return(0);
35: }
37: static PetscErrorCode PCSetUp_LSC(PC pc)
38: {
39: PC_LSC *lsc = (PC_LSC*)pc->data;
40: Mat L,Lp,B,C;
44: PCLSCAllocate_Private(pc);
45: PetscObjectQuery((PetscObject)pc->mat,"LSC_L",(PetscObject*)&L);
46: if (!L) {PetscObjectQuery((PetscObject)pc->pmat,"LSC_L",(PetscObject*)&L);}
47: PetscObjectQuery((PetscObject)pc->pmat,"LSC_Lp",(PetscObject*)&Lp);
48: if (!Lp) {PetscObjectQuery((PetscObject)pc->mat,"LSC_Lp",(PetscObject*)&Lp);}
49: if (!L) {
50: MatSchurComplementGetSubMatrices(pc->mat,NULL,NULL,&B,&C,NULL);
51: if (!lsc->L) {
52: MatMatMult(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lsc->L);
53: } else {
54: MatMatMult(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&lsc->L);
55: }
56: Lp = L = lsc->L;
57: }
58: if (lsc->scale) {
59: Mat Ap;
60: MatSchurComplementGetSubMatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);
61: MatGetDiagonal(Ap,lsc->scale); /* Should be the mass matrix, but we don't have plumbing for that yet */
62: VecReciprocal(lsc->scale);
63: }
64: KSPSetOperators(lsc->kspL,L,Lp);
65: KSPSetFromOptions(lsc->kspL);
66: return(0);
67: }
69: static PetscErrorCode PCApply_LSC(PC pc,Vec x,Vec y)
70: {
71: PC_LSC *lsc = (PC_LSC*)pc->data;
72: Mat A,B,C;
76: MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,&B,&C,NULL);
77: KSPSolve(lsc->kspL,x,lsc->x1);
78: KSPCheckSolve(lsc->kspL,pc,lsc->x1);
79: MatMult(B,lsc->x1,lsc->x0);
80: if (lsc->scale) {
81: VecPointwiseMult(lsc->x0,lsc->x0,lsc->scale);
82: }
83: MatMult(A,lsc->x0,lsc->y0);
84: if (lsc->scale) {
85: VecPointwiseMult(lsc->y0,lsc->y0,lsc->scale);
86: }
87: MatMult(C,lsc->y0,lsc->x1);
88: KSPSolve(lsc->kspL,lsc->x1,y);
89: KSPCheckSolve(lsc->kspL,pc,y);
90: return(0);
91: }
93: static PetscErrorCode PCReset_LSC(PC pc)
94: {
95: PC_LSC *lsc = (PC_LSC*)pc->data;
99: VecDestroy(&lsc->x0);
100: VecDestroy(&lsc->y0);
101: VecDestroy(&lsc->x1);
102: VecDestroy(&lsc->scale);
103: KSPDestroy(&lsc->kspL);
104: MatDestroy(&lsc->L);
105: return(0);
106: }
108: static PetscErrorCode PCDestroy_LSC(PC pc)
109: {
113: PCReset_LSC(pc);
114: PetscFree(pc->data);
115: return(0);
116: }
118: static PetscErrorCode PCSetFromOptions_LSC(PetscOptionItems *PetscOptionsObject,PC pc)
119: {
120: PC_LSC *lsc = (PC_LSC*)pc->data;
124: PetscOptionsHead(PetscOptionsObject,"LSC options");
125: {
126: PetscOptionsBool("-pc_lsc_scale_diag","Use diagonal of velocity block (A) for scaling","None",lsc->scalediag,&lsc->scalediag,NULL);
127: }
128: PetscOptionsTail();
129: return(0);
130: }
132: static PetscErrorCode PCView_LSC(PC pc,PetscViewer viewer)
133: {
134: PC_LSC *jac = (PC_LSC*)pc->data;
136: PetscBool iascii;
139: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
140: if (iascii) {
141: PetscViewerASCIIPushTab(viewer);
142: KSPView(jac->kspL,viewer);
143: PetscViewerASCIIPopTab(viewer);
144: }
145: return(0);
146: }
148: /*MC
149: PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators
151: Options Database Key:
152: . -pc_lsc_scale_diag - Use the diagonal of A for scaling
154: Level: intermediate
156: Notes:
157: This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but
158: it can be used for any Schur complement system. Consider the Schur complement
160: .vb
161: S = A11 - A10 inv(A00) A01
162: .ve
164: PCLSC currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to
165: inv(S) is given by
167: .vb
168: inv(A10 A01) A10 A00 A01 inv(A10 A01)
169: .ve
171: The product A10 A01 can be computed for you, but you can provide it (this is
172: usually more efficient anyway). In the case of incompressible flow, A10 A10 is a Laplacian, call it L. The current
173: interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix.
175: If you had called KSPSetOperators(ksp,S,Sp), S should have type MATSCHURCOMPLEMENT and Sp can be any type you
176: like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
177: For example, you might have setup code like this
179: .vb
180: PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
181: PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
182: .ve
184: And then your Jacobian assembly would look like
186: .vb
187: PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
188: PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
189: if (L) { assembly L }
190: if (Lp) { assemble Lp }
191: .ve
193: With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L
195: .vb
196: -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
197: .ve
199: Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.
201: References:
202: + 1. - Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006.
203: - 2. - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier Stokes equations for incompressible flow, 2001.
205: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCFIELDSPLIT,
206: PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(),
207: MatCreateSchurComplement()
208: M*/
210: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
211: {
212: PC_LSC *lsc;
216: PetscNewLog(pc,&lsc);
217: pc->data = (void*)lsc;
219: pc->ops->apply = PCApply_LSC;
220: pc->ops->applytranspose = 0;
221: pc->ops->setup = PCSetUp_LSC;
222: pc->ops->reset = PCReset_LSC;
223: pc->ops->destroy = PCDestroy_LSC;
224: pc->ops->setfromoptions = PCSetFromOptions_LSC;
225: pc->ops->view = PCView_LSC;
226: pc->ops->applyrichardson = 0;
227: return(0);
228: }