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: }