1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: mlInterface.c,v 1.1 1999/06/09 20:24:12 knepley Exp $";
  3: #endif
  4: /*
  5:    Defines the interface for the  multilevel preconditioner due to
  6:    Vivek Sarin for AIJ matrices.
  7: */
 8:  #include src/sles/pc/pcimpl.h
 9:  #include petscsnes.h
 10:  #include ml.h

 12: #undef  __FUNCT__
 14: /*@C PCMultiLevelSetFields
 15:         This function sets the shape and test function fields for the
 16:   gradient operator. It must be called prior to PCSetUp().

 18:   Collective on PC

 20:   Input Parameters:
 21: + pc     - The preconditioner context
 22: . sField - The shape function field
 23: - tField - The test function field

 25:   Level: intermediate

 27: .keywords multilevel, fields
 28: .seealso PCMultiLevelSetGradientOperator, PCMultiLevelSetNonlinearIterate
 29: @*/
 30: int PCMultiLevelSetFields(PC pc, int sField, int tField)
 31: {
 32:   PC_Multilevel *ml;

 36:   ml = (PC_Multilevel *) pc->data;
 37:   ml->sField = sField;
 38:   ml->tField = tField;
 39:   return(0);
 40: }

 42: #undef  __FUNCT__
 44: /*@C PCMultiLevelSetNonlinearIterate
 45:         This function sets the solution vector used in a Newton iteration. If
 46:   the inner solve is preconditioned with ML, the current iterate is
 47:   necessary to perform the reduction of the system.

 49:   Collective on PC

 51:   Input Parameters:
 52: + pc - The preconditioner context
 53: - u  - The current iterate in a Newton solve

 55:   Level: intermediate

 57: .keywords multilevel, nonlinear
 58: .seealso PCMultiLevelSetFields, PCMultiLevelSetGradientOperator
 59: @*/
 60: int PCMultiLevelSetNonlinearIterate(PC pc, GVec u)
 61: {
 62:   PC_Multilevel *ml;

 67:   ml = (PC_Multilevel *) pc->data;
 68:   ml->nonlinearIterate = u;
 69:   return(0);
 70: }

 72: #undef  __FUNCT__
 74: /*@C PCMultiLevelSetGradientOperator
 75:         This function sets the operator to use as the gradient and its transpose
 76:   in the multilevel scheme. It must be called prior to PCSetUp().

 78:   Collective on PC

 80:   Input Parameters:
 81: + pc    - The preconditioner context
 82: . grad  - The gradient operator
 83: . div   - The divergence operator
 84: - alpha - The scalar multiplier for the gradient operator

 86:   Level: intermediate

 88: .keywords multilevel, operator
 89: .seealso PCMultiLevelSetFields, PCMultiLevelSetNonlinearIterate
 90: @*/
 91: int PCMultiLevelSetGradientOperator(PC pc, int Grad, int Div, PetscScalar alpha)
 92: {
 93:   PC_Multilevel *ml;

 97:   ml = (PC_Multilevel *) pc->data;
 98:   ml->gradOp    = Grad;
 99:   ml->divOp     = Div;
100:   ml->gradAlpha = alpha;
101:   return(0);
102: }

104: #undef  __FUNCT__
106: /*@C PCMultiLevelBuildSolution
107:         This function recovers the solution of the global problem
108:   from the solution of the reduced system.

110:   Collective on GVec

112:   Input Parameters:
113: + pc - The preconditioner context
114: - u  - The solution of the reduced system

116:   Output Parameter:
117: . u  - The global solution

119:   Level: beginner

121: .keywords multilevel, solution
122: .seealso PCMultiLevelApplyGradientGetMultiplier
123: @*/
124: int PCMultiLevelBuildSolution(PC pc, GVec u)
125: {
126:   PC_Multilevel *ml;
127:   Grid           grid;
128:   PetscScalar    alpha;
129:   int            ierr;

134:   if (pc->setupcalled < 2) {
135:     PCSetUp(pc);
136:   }

138:   ml = (PC_Multilevel *) pc->data;
139:   if (ml->reduceSol != PETSC_NULL) {
140:     ierr  = GVecGetGrid(u, &grid);
141:     ierr  = GridGetBCMultiplier(grid, &alpha);
142:     alpha = -alpha;
143:     ierr  = VecAXPY(&alpha, ml->reduceSol, u);
144:   }
145:   return(0);
146: }

148: #undef  __FUNCT__
150: /*@C PCMultiLevelGetMultiplier
151:         This function recovers the multiplier of the global problem
152:   from the solution of the reduced system.

154:   Collective on GVec

156:   Input Parameters:
157: + pc - The preconditioner context
158: - u  - The solution of the reduced system

160:   Output Parameter:
161: . p  - The global multiplier

163:   Level: beginner

165: .keywords multilevel, solution, multiplier
166: .seealso PCMultiLevelBuildSolution
167: @*/
168: int PCMultiLevelGetMultiplier(PC pc, GVec u, GVec p)
169: {
170:   GVec        tempVec;
171:   PetscScalar minusOne = -1.0;
172:   PetscScalar zero     =  0.0;
173:   int         ierr;

179:   if ((pc->vec == PETSC_NULL) || (pc->mat == PETSC_NULL)) {
180:     VecSet(&zero, p);
181:     return(0);
182:   }
183:   if (pc->setupcalled < 2) {
184:     PCSetUp(pc);
185:   }

187:   /* Multiply P_2 (P^{-1} u)_2 by A */
188:   GVecDuplicate(u, &tempVec);
189:   MatMult(pc->mat, u, tempVec);

191:   /* Let A P_2 (P^{-1} u)_2 = f - A P_1 D^{-1} Z^T g - A P_2 (P^{-1} u)_2 --
192:        This is really unstaisfactory, but how do I get the Rhs after solution of
193:        the system. Perhaps I should have a flag for storing the multiplier.
194:   */
195:   VecAYPX(&minusOne, pc->vec, tempVec);

197:   /* Apply P^T_1 */
198:   PCMultiLevelApplyP1Trans(pc, tempVec, p);

200:   /* Apply Z D^{-1} */
201:   PCMultiLevelApplyDInv(pc, p, p);
202:   PCMultiLevelApplyV(pc, p, p);

204:   /* Cleanup */
205:   VecDestroy(tempVec);
206:   return(0);
207: }