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