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: }