Actual source code: wp.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

  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>

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

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

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

 57:   Output Parameter:
 58: .   h - the scale computed

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

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

 86: /*
 87:    MatMFFDView_WP - Prints information about this particular
 88:      method for computing h. Note that this does not print the general
 89:      information about the matrix free, that is printed by the calling
 90:      routine.

 92:   Input Parameters:
 93: +   ctx - the matrix free context
 94: -   viewer - the PETSc viewer

 96: */
 97: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
 98: {
 99:   MatMFFD_WP     *hctx = (MatMFFD_WP*)ctx->hctx;
101:   PetscBool      iascii;

104:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
105:   if (iascii) {
106:     if (hctx->computenormU) {
107:        PetscViewerASCIIPrintf(viewer,"    Computes normU\n");
108:     } else {
109:        PetscViewerASCIIPrintf(viewer,"    Does not compute normU\n");
110:     }
111:   }
112:   return(0);
113: }

115: /*
116:    MatMFFDSetFromOptions_WP - Looks in the options database for
117:      any options appropriate for this method

119:   Input Parameter:
120: .  ctx - the matrix free context

122: */
123: static PetscErrorCode MatMFFDSetFromOptions_WP(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
124: {
126:   MatMFFD_WP     *hctx = (MatMFFD_WP*)ctx->hctx;

129:   PetscOptionsHead(PetscOptionsObject,"Walker-Pernice options");
130:   PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,NULL);
131:   PetscOptionsTail();
132:   return(0);
133: }

135: /*
136:    MatMFFDDestroy_WP - Frees the space allocated by
137:        MatCreateMFFD_WP().

139:   Input Parameter:
140: .  ctx - the matrix free context

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

144: */
145: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
146: {

150:   PetscFree(ctx->hctx);
151:   return(0);
152: }

154: PetscErrorCode  MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)
155: {
156:   MatMFFD    ctx   = (MatMFFD)mat->data;
157:   MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;

160:   hctx->computenormU = flag;
161:   return(0);
162: }

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

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

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

177:   Level: advanced

179:   Notes:
180:    See the manual page for MATMFFD_WP for a complete description of the
181:    algorithm used to compute h.

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

185: @*/
186: PetscErrorCode  MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)
187: {

192:   PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
193:   return(0);
194: }

196: /*
197:      MatCreateMFFD_WP - Standard PETSc code for
198:    computing h with matrix-free finite differences.

200:    Input Parameter:
201: .  ctx - the matrix free context created by MatCreateMFFD()

203: */
204: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
205: {
207:   MatMFFD_WP     *hctx;

210:   /* allocate my own private data structure */
211:   PetscNewLog(ctx,&hctx);
212:   ctx->hctx          = (void*)hctx;
213:   hctx->computenormU = PETSC_FALSE;

215:   /* set the functions I am providing */
216:   ctx->ops->compute        = MatMFFDCompute_WP;
217:   ctx->ops->destroy        = MatMFFDDestroy_WP;
218:   ctx->ops->view           = MatMFFDView_WP;
219:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;

221:   PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C",MatMFFDWPSetComputeNormU_P);
222:   return(0);
223: }