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`
4
5 $$
6 h = error_rel * sqrt(1 + ||u||) / ||a||
7 $$
8
9 Options Database Key:
10 . -mat_mffd_compute_normu -Compute the norm of u every time see `MatMFFDWPSetComputeNormU()`
11
12 Level: intermediate
13
14 Notes:
15 $ || U || $ does not change between linear iterations so is reused
16
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`
19
20 .seealso: `MATMFFD`, `MATMFFD_DS`, `MatCreateMFFD()`, `MatCreateSNESMF()`, `MATMFFD_DS`
21 M*/
22
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.
27
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> /*I "petscmat.h" I*/
32
33 typedef struct {
34 PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
35 PetscBool computenormU;
36 } MatMFFD_WP;
37
38 /*
39 MatMFFDCompute_WP - code for
40 computing h with matrix-free finite differences.
41
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
46
47 Output Parameter:
48 . h - the scale computed
49 */
MatMFFDCompute_WP(MatMFFD ctx,Vec U,Vec a,PetscScalar * h,PetscBool * zeroa)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;
54
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 }
73
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.
79
80 Input Parameters:
81 + ctx - the matrix-free context
82 - viewer - the PETSc viewer
83
84 */
MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)85 static PetscErrorCode MatMFFDView_WP(MatMFFD ctx, PetscViewer viewer)
86 {
87 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
88 PetscBool isascii;
89
90 PetscFunctionBegin;
91 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
92 if (isascii) {
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 }
101
102 /*
103 MatMFFDSetFromOptions_WP - Looks in the options database for
104 any options appropriate for this method
105
106 Input Parameter:
107 . ctx - the matrix-free context
108
109 */
MatMFFDSetFromOptions_WP(MatMFFD ctx,PetscOptionItems PetscOptionsObject)110 static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems PetscOptionsObject)
111 {
112 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
113
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 }
120
MatMFFDDestroy_WP(MatMFFD ctx)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 }
128
MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)129 static PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
130 {
131 MatMFFD ctx = (MatMFFD)mat->data;
132 MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
133
134 PetscFunctionBegin;
135 hctx->computenormU = flag;
136 PetscFunctionReturn(PETSC_SUCCESS);
137 }
138
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.
143
144 Input Parameters:
145 + A - the `MATMFFD` matrix
146 - flag - `PETSC_TRUE` causes it to compute $||U||$, `PETSC_FALSE` uses the previous value
147
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.
151
152 Level: advanced
153
154 Note:
155 See the manual page for `MATMFFD_WP` for a complete description of the
156 algorithm used to compute h.
157
158 .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
159 @*/
MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)160 PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
161 {
162 PetscFunctionBegin;
163 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
164 PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
165 PetscFunctionReturn(PETSC_SUCCESS);
166 }
167
168 /*
169 MatCreateMFFD_WP - Standard PETSc code for
170 computing h with matrix-free finite differences.
171
172 Input Parameter:
173 . ctx - the matrix-free context created by MatCreateMFFD()
174
175 */
MatCreateMFFD_WP(MatMFFD ctx)176 PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
177 {
178 MatMFFD_WP *hctx;
179
180 PetscFunctionBegin;
181 /* allocate my own private data structure */
182 PetscCall(PetscNew(&hctx));
183 ctx->hctx = (void *)hctx;
184 hctx->computenormU = PETSC_FALSE;
185
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;
191
192 PetscCall(PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P));
193 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195