2: /*MC
  3:      MATMFFD_WP - Implements an alternative approach for computing the differencing parameter
  4:         h used with the finite difference based matrix-free Jacobian.  This code
  5:         implements the strategy of M. Pernice and H. Walker:

  7:       h = error_rel * sqrt(1 + ||U||) / ||a||

  9:       Notes:
 10:         1) || U || does not change between linear iterations so is reused
 11:         2) In GMRES || a || == 1 and so does not need to ever be computed except at restart
 12:            when it is recomputed.

 14:       Reference:  M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative 
 15:       Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998, 
 16:       vol 19, pp. 302--318.

 18:    Options Database Keys:
 19: .   -mat_mffd_compute_normu -Compute the norm of u everytime see MatMFFDWPSetComputeNormU()


 22:    Level: intermediate

 24:    Notes: Requires no global collectives when used with GMRES

 26:    Formula used:
 27:      F'(u)*a = [F(u+h*a) - F(u)]/h where

 29: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_DS

 31: M*/

 33: /*
 34:     This include file defines the data structure  MatMFFD that 
 35:    includes information about the computation of h. It is shared by 
 36:    all implementations that people provide.

 38:    See snesmfjdef.c for  a full set of comments on the routines below.
 39: */
 40: #include <petsc-private/matimpl.h>
 41: #include <../src/mat/impls/mffd/mffdimpl.h>   /*I  "petscmat.h"   I*/

 43: typedef struct {
 44:   PetscReal  normUfact;                   /* previous sqrt(1.0 + || U ||) */
 45:   PetscBool  computenormU;
 46: } MatMFFD_WP;

 50: /*
 51:      MatMFFDCompute_WP - Standard PETSc code for 
 52:    computing h with matrix-free finite differences.

 54:   Input Parameters:
 55: +   ctx - the matrix free context
 56: .   U - the location at which you want the Jacobian
 57: -   a - the direction you want the derivative

 59:   Output Parameter:
 60: .   h - the scale computed

 62: */
 63: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool  *zeroa)
 64: {
 65:   MatMFFD_WP    *hctx = (MatMFFD_WP*)ctx->hctx;
 66:   PetscReal      normU,norma;

 70:   if (!(ctx->count % ctx->recomputeperiod)) {
 71:     if (hctx->computenormU || !ctx->ncurrenth) {
 72:       VecNorm(U,NORM_2,&normU);
 73:       hctx->normUfact = PetscSqrtReal(1.0+normU);
 74:     }
 75:     VecNorm(a,NORM_2,&norma);
 76:     if (norma == 0.0) {
 77:       *zeroa = PETSC_TRUE;
 78:       return(0);
 79:     }
 80:     *zeroa = PETSC_FALSE;
 81:     *h = ctx->error_rel*hctx->normUfact/norma;
 82:   } else {
 83:     *h = ctx->currenth;
 84:   }
 85:   return(0);
 86: }

 90: /*
 91:    MatMFFDView_WP - Prints information about this particular 
 92:      method for computing h. Note that this does not print the general
 93:      information about the matrix free, that is printed by the calling
 94:      routine.

 96:   Input Parameters:
 97: +   ctx - the matrix free context
 98: -   viewer - the PETSc viewer

100: */
101: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
102: {
103:   MatMFFD_WP     *hctx = (MatMFFD_WP *)ctx->hctx;
105:   PetscBool      iascii;

108:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
109:   if (iascii) {
110:     if (hctx->computenormU){ PetscViewerASCIIPrintf(viewer,"    Computes normU\n");}
111:     else                   { PetscViewerASCIIPrintf(viewer,"    Does not compute normU\n");}
112:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES matrix-free WP",((PetscObject)viewer)->type_name);
113:   return(0);
114: }

118: /*
119:    MatMFFDSetFromOptions_WP - Looks in the options database for 
120:      any options appropriate for this method

122:   Input Parameter:
123: .  ctx - the matrix free context

125: */
126: static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx)
127: {
129:   MatMFFD_WP     *hctx = (MatMFFD_WP*)ctx->hctx;

132:   PetscOptionsHead("Walker-Pernice options");
133:     PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,0);
134:   PetscOptionsTail();
135:   return(0);
136: }

140: /*
141:    MatMFFDDestroy_WP - Frees the space allocated by 
142:        MatCreateMFFD_WP(). 

144:   Input Parameter:
145: .  ctx - the matrix free context

147:    Notes: does not free the ctx, that is handled by the calling routine

149: */
150: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
151: {
154:   PetscFree(ctx->hctx);
155:   return(0);
156: }

158: EXTERN_C_BEGIN
161: PetscErrorCode  MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool  flag)
162: {
163:   MatMFFD     ctx = (MatMFFD)mat->data;
164:   MatMFFD_WP  *hctx = (MatMFFD_WP*)ctx->hctx;

167:   hctx->computenormU = flag;
168:   return(0);
169: }
170: EXTERN_C_END

174: /*@
175:     MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
176:              PETSc routine for computing h. With any Krylov solver this need only 
177:              be computed during the first iteration and kept for later.

179:   Input Parameters:
180: +   A - the matrix created with MatCreateSNESMF()
181: -   flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value

183:   Options Database Key:
184: .   -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you 
185:               must be sure that ||U|| has not changed in the mean time.

187:   Level: advanced

189:   Notes:
190:    See the manual page for MATMFFD_WP for a complete description of the
191:    algorithm used to compute h.

193: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()

195: @*/
196: PetscErrorCode  MatMFFDWPSetComputeNormU(Mat A,PetscBool  flag)
197: {

202:   PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
203:   return(0);
204: }

206: EXTERN_C_BEGIN
209: /*
210:      MatCreateMFFD_WP - Standard PETSc code for 
211:    computing h with matrix-free finite differences.

213:    Input Parameter:
214: .  ctx - the matrix free context created by MatCreateMFFD()

216: */
217: PetscErrorCode  MatCreateMFFD_WP(MatMFFD ctx)
218: {
220:   MatMFFD_WP     *hctx;


224:   /* allocate my own private data structure */
225:   PetscNewLog(ctx,MatMFFD_WP,&hctx);
226:   ctx->hctx          = (void*)hctx;
227:   hctx->computenormU = PETSC_FALSE;

229:   /* set the functions I am providing */
230:   ctx->ops->compute        = MatMFFDCompute_WP;
231:   ctx->ops->destroy        = MatMFFDDestroy_WP;
232:   ctx->ops->view           = MatMFFDView_WP;
233:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;

235:   PetscObjectComposeFunctionDynamic((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C","MatMFFDWPSetComputeNormU_P",MatMFFDWPSetComputeNormU_P);

237:   return(0);
238: }
239: EXTERN_C_END