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: if (!(ctx->count % ctx->recomputeperiod)) {
64: if (hctx->computenormU || !ctx->ncurrenth) {
65: VecNorm(U, NORM_2, &normU);
66: hctx->normUfact = PetscSqrtReal(1.0 + normU);
67: }
68: VecNorm(a, NORM_2, &norma);
69: if (norma == 0.0) {
70: *zeroa = PETSC_TRUE;
71: return 0;
72: }
73: *zeroa = PETSC_FALSE;
74: *h = ctx->error_rel * hctx->normUfact / norma;
75: } else {
76: *h = ctx->currenth;
77: }
78: return 0;
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: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
98: if (iascii) {
99: if (hctx->computenormU) {
100: PetscViewerASCIIPrintf(viewer, " Computes normU\n");
101: } else {
102: PetscViewerASCIIPrintf(viewer, " Does not compute normU\n");
103: }
104: }
105: return 0;
106: }
108: /*
109: MatMFFDSetFromOptions_WP - Looks in the options database for
110: any options appropriate for this method
112: Input Parameter:
113: . ctx - the matrix free context
115: */
116: static PetscErrorCode MatMFFDSetFromOptions_WP(MatMFFD ctx, PetscOptionItems *PetscOptionsObject)
117: {
118: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
120: PetscOptionsHeadBegin(PetscOptionsObject, "Walker-Pernice options");
121: PetscOptionsBool("-mat_mffd_compute_normu", "Compute the norm of u", "MatMFFDWPSetComputeNormU", hctx->computenormU, &hctx->computenormU, NULL);
122: PetscOptionsHeadEnd();
123: return 0;
124: }
126: /*
127: MatMFFDDestroy_WP - Frees the space allocated by
128: MatCreateMFFD_WP().
130: Input Parameter:
131: . ctx - the matrix free context
133: Notes:
134: does not free the ctx, that is handled by the calling routine
136: */
137: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
138: {
139: PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", NULL);
140: PetscFree(ctx->hctx);
141: return 0;
142: }
144: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat, PetscBool flag)
145: {
146: MatMFFD ctx = (MatMFFD)mat->data;
147: MatMFFD_WP *hctx = (MatMFFD_WP *)ctx->hctx;
149: hctx->computenormU = flag;
150: return 0;
151: }
153: /*@
154: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the Walker-Pernice
155: PETSc routine for computing h. With any Krylov solver this need only
156: be computed during the first iteration and kept for later.
158: Input Parameters:
159: + A - the `MATMFFD` matrix
160: - flag - `PETSC_TRUE` causes it to compute ||U||, `PETSC_FALSE` uses the previous value
162: Options Database Key:
163: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
164: must be sure that ||U|| has not changed in the mean time.
166: Level: advanced
168: Note:
169: See the manual page for `MATMFFD_WP` for a complete description of the
170: algorithm used to compute h.
172: .seealso: `MATMFFD_WP`, `MATMFFD`, `MatMFFDSetFunctionError()`, `MatCreateSNESMF()`
173: @*/
174: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A, PetscBool flag)
175: {
177: PetscTryMethod(A, "MatMFFDWPSetComputeNormU_C", (Mat, PetscBool), (A, flag));
178: return 0;
179: }
181: /*
182: MatCreateMFFD_WP - Standard PETSc code for
183: computing h with matrix-free finite differences.
185: Input Parameter:
186: . ctx - the matrix free context created by MatCreateMFFD()
188: */
189: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
190: {
191: MatMFFD_WP *hctx;
193: /* allocate my own private data structure */
194: PetscNew(&hctx);
195: ctx->hctx = (void *)hctx;
196: hctx->computenormU = PETSC_FALSE;
198: /* set the functions I am providing */
199: ctx->ops->compute = MatMFFDCompute_WP;
200: ctx->ops->destroy = MatMFFDDestroy_WP;
201: ctx->ops->view = MatMFFDView_WP;
202: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
204: PetscObjectComposeFunction((PetscObject)ctx->mat, "MatMFFDWPSetComputeNormU_C", MatMFFDWPSetComputeNormU_P);
205: return 0;
206: }