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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

160:   Level: advanced

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

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

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

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

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

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

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

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