Actual source code: wp.c
petsc-3.13.6 2020-09-29
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 everytime see MatMFFDWPSetComputeNormU()
22: Level: intermediate
24: Notes:
25: Requires no global collectives when used with GMRES
27: Formula used:
28: F'(u)*a = [F(u+h*a) - F(u)]/h where
30: .seealso: MATMFFD, MatCreateMFFD(), MatCreateSNESMF(), MATMFFD_DS
32: M*/
34: /*
35: This include file defines the data structure MatMFFD that
36: includes information about the computation of h. It is shared by
37: all implementations that people provide.
39: See snesmfjdef.c for a full set of comments on the routines below.
40: */
41: #include <petsc/private/matimpl.h>
42: #include <../src/mat/impls/mffd/mffdimpl.h>
44: typedef struct {
45: PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
46: PetscBool computenormU;
47: } MatMFFD_WP;
49: /*
50: MatMFFDCompute_WP - Standard PETSc code for
51: computing h with matrix-free finite differences.
53: Input Parameters:
54: + ctx - the matrix free context
55: . U - the location at which you want the Jacobian
56: - a - the direction you want the derivative
58: Output Parameter:
59: . h - the scale computed
61: */
62: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa)
63: {
64: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
65: PetscReal normU,norma;
69: if (!(ctx->count % ctx->recomputeperiod)) {
70: if (hctx->computenormU || !ctx->ncurrenth) {
71: VecNorm(U,NORM_2,&normU);
72: hctx->normUfact = PetscSqrtReal(1.0+normU);
73: }
74: VecNorm(a,NORM_2,&norma);
75: if (norma == 0.0) {
76: *zeroa = PETSC_TRUE;
77: return(0);
78: }
79: *zeroa = PETSC_FALSE;
80: *h = ctx->error_rel*hctx->normUfact/norma;
81: } else {
82: *h = ctx->currenth;
83: }
84: return(0);
85: }
87: /*
88: MatMFFDView_WP - Prints information about this particular
89: method for computing h. Note that this does not print the general
90: information about the matrix free, that is printed by the calling
91: routine.
93: Input Parameters:
94: + ctx - the matrix free context
95: - viewer - the PETSc viewer
97: */
98: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
99: {
100: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
102: PetscBool iascii;
105: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
106: if (iascii) {
107: if (hctx->computenormU) {
108: PetscViewerASCIIPrintf(viewer," Computes normU\n");
109: } else {
110: PetscViewerASCIIPrintf(viewer," Does not compute normU\n");
111: }
112: }
113: return(0);
114: }
116: /*
117: MatMFFDSetFromOptions_WP - Looks in the options database for
118: any options appropriate for this method
120: Input Parameter:
121: . ctx - the matrix free context
123: */
124: static PetscErrorCode MatMFFDSetFromOptions_WP(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
125: {
127: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
130: PetscOptionsHead(PetscOptionsObject,"Walker-Pernice options");
131: PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,NULL);
132: PetscOptionsTail();
133: return(0);
134: }
136: /*
137: MatMFFDDestroy_WP - Frees the space allocated by
138: MatCreateMFFD_WP().
140: Input Parameter:
141: . ctx - the matrix free context
143: Notes:
144: does not free the ctx, that is handled by the calling routine
146: */
147: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
148: {
152: PetscFree(ctx->hctx);
153: return(0);
154: }
156: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)
157: {
158: MatMFFD ctx = (MatMFFD)mat->data;
159: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
162: hctx->computenormU = flag;
163: return(0);
164: }
166: /*@
167: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
168: PETSc routine for computing h. With any Krylov solver this need only
169: be computed during the first iteration and kept for later.
171: Input Parameters:
172: + A - the matrix created with MatCreateSNESMF()
173: - flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value
175: Options Database Key:
176: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
177: must be sure that ||U|| has not changed in the mean time.
179: Level: advanced
181: Notes:
182: See the manual page for MATMFFD_WP for a complete description of the
183: algorithm used to compute h.
185: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
187: @*/
188: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)
189: {
194: PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
195: return(0);
196: }
198: /*
199: MatCreateMFFD_WP - Standard PETSc code for
200: computing h with matrix-free finite differences.
202: Input Parameter:
203: . ctx - the matrix free context created by MatCreateMFFD()
205: */
206: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
207: {
209: MatMFFD_WP *hctx;
212: /* allocate my own private data structure */
213: PetscNewLog(ctx,&hctx);
214: ctx->hctx = (void*)hctx;
215: hctx->computenormU = PETSC_FALSE;
217: /* set the functions I am providing */
218: ctx->ops->compute = MatMFFDCompute_WP;
219: ctx->ops->destroy = MatMFFDDestroy_WP;
220: ctx->ops->view = MatMFFDView_WP;
221: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
223: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C",MatMFFDWPSetComputeNormU_P);
224: return(0);
225: }