Actual source code: wp.c

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

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

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

 10:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

 62:   PetscFunctionBegin;
 63:   if (!(ctx->count % ctx->recomputeperiod)) {
 64:     if (hctx->computenormU || !ctx->ncurrenth) {
 65:       PetscCall(VecNorm(U, NORM_2, &normU));
 66:       hctx->normUfact = PetscSqrtReal(1.0 + normU);
 67:     }
 68:     PetscCall(VecNorm(a, NORM_2, &norma));
 69:     if (norma == 0.0) {
 70:       *zeroa = PETSC_TRUE;
 71:       PetscFunctionReturn(PETSC_SUCCESS);
 72:     }
 73:     *zeroa = PETSC_FALSE;
 74:     *h     = ctx->error_rel * hctx->normUfact / norma;
 75:   } else {
 76:     *h = ctx->currenth;
 77:   }
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 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:   PetscFunctionBegin;
 98:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 99:   if (iascii) {
100:     if (hctx->computenormU) {
101:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Computes normU\n"));
102:     } else {
103:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Does not compute normU\n"));
104:     }
105:   }
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

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

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

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

121:   PetscFunctionBegin;
122:   PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options");
123:   PetscCall(PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL));
124:   PetscOptionsHeadEnd();
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
129: {
130:   PetscFunctionBegin;
131:   PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL));
132:   PetscCall(PetscFree(ctx->hctx));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: static PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
137: {
138:   MatMFFD     ctx  = (MatMFFD)mat->data;
139:   MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;

141:   PetscFunctionBegin;
142:   hctx->computenormU = flag;
143:   PetscFunctionReturn(PETSC_SUCCESS);
144: }

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

151:   Input Parameters:
152: + A    - the `MATMFFD` matrix
153: - flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value

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

159:   Level: advanced

161:   Note:
162:   See the manual page for `MATMFFD_WP` for a complete description of the
163:   algorithm used to compute h.

165: .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
166: @*/
167: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
168: {
169:   PetscFunctionBegin;
171:   PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: /*
176:      MatCreateMFFD_WP - Standard PETSc code for
177:    computing h with matrix-free finite differences.

179:    Input Parameter:
180: .  ctx - the matrix-free context created by MatCreateMFFD()

182: */
183: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
184: {
185:   MatMFFD_WP *hctx;

187:   PetscFunctionBegin;
188:   /* allocate my own private data structure */
189:   PetscCall(PetscNew(&hctx));
190:   ctx->hctx          = (void *)hctx;
191:   hctx->computenormU = PETSC_FALSE;

193:   /* set the functions I am providing */
194:   ctx->ops->compute        = MatMFFDCompute_WP;
195:   ctx->ops->destroy        = MatMFFDDestroy_WP;
196:   ctx->ops->view           = MatMFFDView_WP;
197:   ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;

199:   PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P));
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }