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