Actual source code: lsc.c

petsc-master 2015-07-30
Report Typos and Errors
  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: }