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. From Walker-Pernice {cite}`pw98`
5: $$
6: h = error_rel * sqrt(1 + ||u||) / ||a||
7: $$
9: Options Database Key:
10: . -mat_mffd_compute_normu -Compute the norm of u every time see `MatMFFDWPSetComputeNormU()`
12: Level: intermediate
14: Notes:
15: $ || U || $ does not change between linear iterations so is reused
17: In `KSPGMRES` $ || a || == 1 $ and so does not need to ever be computed except at restart
18: when it is recomputed. Thus requires no global collectives when used with `KSPGMRES`
20: .seealso: `MATMFFD`, `MATMFFD_DS`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_DS`
21: M*/
23: /*
24: This include file defines the data structure MatMFFD that
25: includes information about the computation of h. It is shared by
26: all implementations that people provide.
28: See snesmfjdef.c for a full set of comments on the routines below.
29: */
30: #include <petsc/private/matimpl.h>
31: #include <../src/mat/impls/mffd/mffdimpl.h>
33: typedef struct {
34: PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
35: PetscBool computenormU;
36: } MatMFFD_WP;
38: /*
39: MatMFFDCompute_WP - code for
40: computing h with matrix-free finite differences.
42: Input Parameters:
43: + ctx - the matrix-free context
44: . U - the location at which you want the Jacobian
45: - a - the direction you want the derivative
47: Output Parameter:
48: . h - the scale computed
49: */
50: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx, Vec U, Vec a, PetscScalar *h, PetscBool *zeroa)
51: {
52: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
53: PetscReal normU, norma;
55: PetscFunctionBegin;
56: if (!(ctx->count % ctx->recomputeperiod)) {
57: if (hctx->computenormU || !ctx->ncurrenth) {
58: PetscCall(VecNorm(U, NORM_2, &normU));
59: hctx->normUfact = PetscSqrtReal(1.0 + normU);
60: }
61: PetscCall(VecNorm(a, NORM_2, &norma));
62: if (norma == 0.0) {
63: *zeroa = PETSC_TRUE;
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
66: *zeroa = PETSC_FALSE;
67: *h = ctx->error_rel * hctx->normUfact / norma;
68: } else {
69: *h = ctx->currenth;
70: }
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*
75: MatMFFDView_WP - Prints information about this particular
76: method for computing h. Note that this does not print the general
77: information about the matrix-free, that is printed by the calling
78: routine.
80: Input Parameters:
81: + ctx - the matrix-free context
82: - viewer - the PETSc viewer
84: */
85: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer)
86: {
87: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
88: PetscBool iascii;
90: PetscFunctionBegin;
91: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
92: if (iascii) {
93: if (hctx->computenormU) {
94: PetscCall(PetscViewerASCIIPrintf(viewer, " Computes normU\n"));
95: } else {
96: PetscCall(PetscViewerASCIIPrintf(viewer, " Does not compute normU\n"));
97: }
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*
103: MatMFFDSetFromOptions_WP - Looks in the options database for
104: any options appropriate for this method
106: Input Parameter:
107: . ctx - the matrix-free context
109: */
110: static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
111: {
112: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
114: PetscFunctionBegin;
115: PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options");
116: PetscCall(PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL));
117: PetscOptionsHeadEnd();
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
122: {
123: PetscFunctionBegin;
124: PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL));
125: PetscCall(PetscFree(ctx->hctx));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
130: {
131: MatMFFD ctx = (MatMFFD)mat->data;
132: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
134: PetscFunctionBegin;
135: hctx->computenormU = flag;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: /*@
140: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice {cite}`pw98`
141: PETSc routine for computing h. With any Krylov solver this need only
142: be computed during the first iteration and kept for later.
144: Input Parameters:
145: + A - the `MATMFFD` matrix
146: - flag - `PETSC_TRUE` causes it to compute $||U||$, `PETSC_FALSE` uses the previous value
148: Options Database Key:
149: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
150: must be sure that $||U||$ has not changed in the mean time.
152: Level: advanced
154: Note:
155: See the manual page for `MATMFFD_WP` for a complete description of the
156: algorithm used to compute h.
158: .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
159: @*/
160: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
161: {
162: PetscFunctionBegin;
164: PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*
169: MatCreateMFFD_WP - Standard PETSc code for
170: computing h with matrix-free finite differences.
172: Input Parameter:
173: . ctx - the matrix-free context created by MatCreateMFFD()
175: */
176: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
177: {
178: MatMFFD_WP *hctx;
180: PetscFunctionBegin;
181: /* allocate my own private data structure */
182: PetscCall(PetscNew(&hctx));
183: ctx->hctx = (void *)hctx;
184: hctx->computenormU = PETSC_FALSE;
186: /* set the functions I am providing */
187: ctx->ops->compute = MatMFFDCompute_WP;
188: ctx->ops->destroy = MatMFFDDestroy_WP;
189: ctx->ops->view = MatMFFDView_WP;
190: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
192: PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }