Actual source code: wp.c
2: /*MC
3: MATMFFD_WP - Implements an alternative approach for computing the differencing parameter
4: h used with the finite difference based matrix-free Jacobian. This code
5: implements the strategy of M. Pernice and H. Walker:
7: h = error_rel * sqrt(1 + ||U||) / ||a||
9: Notes:
10: 1) || U || does not change between linear iterations so is reused
11: 2) In GMRES || a || == 1 and so does not need to ever be computed except at restart
12: when it is recomputed.
14: Reference: M. Pernice and H. F. Walker, "NITSOL: A Newton Iterative
15: Solver for Nonlinear Systems", SIAM J. Sci. Stat. Comput.", 1998,
16: vol 19, pp. 302--318.
18: Options Database Keys:
19: . -mat_mffd_compute_normu -Compute the norm of u every time see MatMFFDWPSetComputeNormU()
21: Level: intermediate
23: Notes:
24: Requires no global collectives when used with GMRES
26: Formula used:
27: F'(u)*a = [F(u+h*a) - F(u)]/h where
29: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_DS
31: M*/
33: /*
34: This include file defines the data structure MatMFFD that
35: includes information about the computation of h. It is shared by
36: all implementations that people provide.
38: See snesmfjdef.c for a full set of comments on the routines below.
39: */
40: #include <petsc/private/matimpl.h>
41: #include <../src/mat/impls/mffd/mffdimpl.h>
43: typedef struct {
44: PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
45: PetscBool computenormU;
46: } MatMFFD_WP;
48: /*
49: MatMFFDCompute_WP - Standard PETSc code for
50: computing h with matrix-free finite differences.
52: Input Parameters:
53: + ctx - the matrix free context
54: . U - the location at which you want the Jacobian
55: - a - the direction you want the derivative
57: Output Parameter:
58: . h - the scale computed
60: */
61: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa)
62: {
63: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
64: PetscReal normU,norma;
66: if (!(ctx->count % ctx->recomputeperiod)) {
67: if (hctx->computenormU || !ctx->ncurrenth) {
68: VecNorm(U,NORM_2,&normU);
69: hctx->normUfact = PetscSqrtReal(1.0+normU);
70: }
71: VecNorm(a,NORM_2,&norma);
72: if (norma == 0.0) {
73: *zeroa = PETSC_TRUE;
74: return 0;
75: }
76: *zeroa = PETSC_FALSE;
77: *h = ctx->error_rel*hctx->normUfact/norma;
78: } else {
79: *h = ctx->currenth;
80: }
81: return 0;
82: }
84: /*
85: MatMFFDView_WP - Prints information about this particular
86: method for computing h. Note that this does not print the general
87: information about the matrix free, that is printed by the calling
88: routine.
90: Input Parameters:
91: + ctx - the matrix free context
92: - viewer - the PETSc viewer
94: */
95: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
96: {
97: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
98: PetscBool iascii;
100: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
101: if (iascii) {
102: if (hctx->computenormU) {
103: PetscViewerASCIIPrintf(viewer," Computes normU\n");
104: } else {
105: PetscViewerASCIIPrintf(viewer," Does not compute normU\n");
106: }
107: }
108: return 0;
109: }
111: /*
112: MatMFFDSetFromOptions_WP - Looks in the options database for
113: any options appropriate for this method
115: Input Parameter:
116: . ctx - the matrix free context
118: */
119: static PetscErrorCode MatMFFDSetFromOptions_WP(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
120: {
121: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
123: PetscOptionsHead(PetscOptionsObject,"Walker-Pernice options");
124: PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,NULL);
125: PetscOptionsTail();
126: return 0;
127: }
129: /*
130: MatMFFDDestroy_WP - Frees the space allocated by
131: MatCreateMFFD_WP().
133: Input Parameter:
134: . ctx - the matrix free context
136: Notes:
137: does not free the ctx, that is handled by the calling routine
139: */
140: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
141: {
142: PetscFree(ctx->hctx);
143: return 0;
144: }
146: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)
147: {
148: MatMFFD ctx = (MatMFFD)mat->data;
149: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
151: hctx->computenormU = flag;
152: return 0;
153: }
155: /*@
156: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
157: PETSc routine for computing h. With any Krylov solver this need only
158: be computed during the first iteration and kept for later.
160: Input Parameters:
161: + A - the matrix created with MatCreateSNESMF()
162: - flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value
164: Options Database Key:
165: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
166: must be sure that ||U|| has not changed in the mean time.
168: Level: advanced
170: Notes:
171: See the manual page for MATMFFD_WP for a complete description of the
172: algorithm used to compute h.
174: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
176: @*/
177: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)
178: {
180: PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
181: return 0;
182: }
184: /*
185: MatCreateMFFD_WP - Standard PETSc code for
186: computing h with matrix-free finite differences.
188: Input Parameter:
189: . ctx - the matrix free context created by MatCreateMFFD()
191: */
192: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
193: {
194: MatMFFD_WP *hctx;
196: /* allocate my own private data structure */
197: PetscNewLog(ctx,&hctx);
198: ctx->hctx = (void*)hctx;
199: hctx->computenormU = PETSC_FALSE;
201: /* set the functions I am providing */
202: ctx->ops->compute = MatMFFDCompute_WP;
203: ctx->ops->destroy = MatMFFDDestroy_WP;
204: ctx->ops->view = MatMFFDView_WP;
205: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
207: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C",MatMFFDWPSetComputeNormU_P);
208: return 0;
209: }