Actual source code: wp.c


  2: /*MC
  3:      MATMFFD_WP - Implements an approach for computing the differencing parameter
  4:         h used with the finite difference based matrix-free Jacobian.

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

  8:    Options Database Key:
  9: .   -mat_mffd_compute_normu -Compute the norm of u every time see `MatMFFDWPSetComputeNormU()`

 11:    Level: intermediate

 13:    Notes:
 14:    || U || does not change between linear iterations so is reused

 16:    In `KSPGMRES` || a || == 1 and so does not need to ever be computed except at restart
 17:     when it is recomputed.  Thus equires no global collectives when used with `KSPGMRES`

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

 22:    Reference:
 23: .  * -  M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative
 24:       Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998,
 25:       vol 19, pp. 302--318.

 27: .seealso: `MATMFFD`, `MATMFFD_DS`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_DS`
 28: M*/

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

 35:    See snesmfjdef.c for  a full set of comments on the routines below.
 36: */
 37: #include <petsc/private/matimpl.h>
 38: #include <../src/mat/impls/mffd/mffdimpl.h>

 40: typedef struct {
 41:   PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
 42:   PetscBool computenormU;
 43: } MatMFFD_WP;

 45: /*
 46:      MatMFFDCompute_WP - code for
 47:    computing h with matrix-free finite differences.

 49:   Input Parameters:
 50: +   ctx - the matrix free context
 51: .   U - the location at which you want the Jacobian
 52: -   a - the direction you want the derivative

 54:   Output Parameter:
 55: .   h - the scale computed

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

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

 81: /*
 82:    MatMFFDView_WP - Prints information about this particular
 83:      method for computing h. Note that this does not print the general
 84:      information about the matrix free, that is printed by the calling
 85:      routine.

 87:   Input Parameters:
 88: +   ctx - the matrix free context
 89: -   viewer - the PETSc viewer

 91: */
 92: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer)
 93: {
 94:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
 95:   PetscBool   iascii;

 97:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 98:   if (iascii) {
 99:     if (hctx->computenormU) {
100:       PetscViewerASCIIPrintf(viewer, "    Computes normU\n");
101:     } else {
102:       PetscViewerASCIIPrintf(viewer, "    Does not compute normU\n");
103:     }
104:   }
105:   return 0;
106: }

108: /*
109:    MatMFFDSetFromOptions_WP - Looks in the options database for
110:      any options appropriate for this method

112:   Input Parameter:
113: .  ctx - the matrix free context

115: */
116: static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
117: {
118:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;

120:   PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options");
121:   PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL);
122:   PetscOptionsHeadEnd();
123:   return 0;
124: }

126: /*
127:    MatMFFDDestroy_WP - Frees the space allocated by
128:        MatCreateMFFD_WP().

130:   Input Parameter:
131: .  ctx - the matrix free context

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

136: */
137: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
138: {
139:   PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL);
140:   PetscFree(ctx->hctx);
141:   return 0;
142: }

144: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
145: {
146:   MatMFFD     ctx  = (MatMFFD)mat->data;
147:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;

149:   hctx->computenormU = flag;
150:   return 0;
151: }

153: /*@
154:     MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice
155:              PETSc routine for computing h. With any Krylov solver this need only
156:              be computed during the first iteration and kept for later.

158:   Input Parameters:
159: +   A - the `MATMFFD` matrix
160: -   flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value

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

166:   Level: advanced

168:   Note:
169:    See the manual page for `MATMFFD_WP` for a complete description of the
170:    algorithm used to compute h.

172: .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
173: @*/
174: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
175: {
177:   PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
178:   return 0;
179: }

181: /*
182:      MatCreateMFFD_WP - Standard PETSc code for
183:    computing h with matrix-free finite differences.

185:    Input Parameter:
186: .  ctx - the matrix free context created by MatCreateMFFD()

188: */
189: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
190: {
191:   MatMFFD_WP *hctx;

193:   /* allocate my own private data structure */
194:   PetscNew(&hctx);
195:   ctx->hctx          = (void *)hctx;
196:   hctx->computenormU = PETSC_FALSE;

198:   /* set the functions I am providing */
199:   ctx->ops->compute        = MatMFFDCompute_WP;
200:   ctx->ops->destroy        = MatMFFDDestroy_WP;
201:   ctx->ops->view           = MatMFFDView_WP;
202:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;

204:   PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P);
205:   return 0;
206: }