Actual source code: wp.c
petsc-3.7.5 2017-01-01
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: 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> /*I "petscmat.h" I*/
43: typedef struct {
44: PetscReal normUfact; /* previous sqrt(1.0 + || U ||) */
45: PetscBool computenormU;
46: } MatMFFD_WP;
50: /*
51: MatMFFDCompute_WP - Standard PETSc code for
52: computing h with matrix-free finite differences.
54: Input Parameters:
55: + ctx - the matrix free context
56: . U - the location at which you want the Jacobian
57: - a - the direction you want the derivative
59: Output Parameter:
60: . h - the scale computed
62: */
63: static PetscErrorCode MatMFFDCompute_WP(MatMFFD ctx,Vec U,Vec a,PetscScalar *h,PetscBool *zeroa)
64: {
65: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
66: PetscReal normU,norma;
70: if (!(ctx->count % ctx->recomputeperiod)) {
71: if (hctx->computenormU || !ctx->ncurrenth) {
72: VecNorm(U,NORM_2,&normU);
73: hctx->normUfact = PetscSqrtReal(1.0+normU);
74: }
75: VecNorm(a,NORM_2,&norma);
76: if (norma == 0.0) {
77: *zeroa = PETSC_TRUE;
78: return(0);
79: }
80: *zeroa = PETSC_FALSE;
81: *h = ctx->error_rel*hctx->normUfact/norma;
82: } else {
83: *h = ctx->currenth;
84: }
85: return(0);
86: }
90: /*
91: MatMFFDView_WP - Prints information about this particular
92: method for computing h. Note that this does not print the general
93: information about the matrix free, that is printed by the calling
94: routine.
96: Input Parameters:
97: + ctx - the matrix free context
98: - viewer - the PETSc viewer
100: */
101: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
102: {
103: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
105: PetscBool iascii;
108: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
109: if (iascii) {
110: if (hctx->computenormU) {
111: PetscViewerASCIIPrintf(viewer," Computes normU\n");
112: } else {
113: PetscViewerASCIIPrintf(viewer," Does not compute normU\n");
114: }
115: }
116: return(0);
117: }
121: /*
122: MatMFFDSetFromOptions_WP - Looks in the options database for
123: any options appropriate for this method
125: Input Parameter:
126: . ctx - the matrix free context
128: */
129: static PetscErrorCode MatMFFDSetFromOptions_WP(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
130: {
132: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
135: PetscOptionsHead(PetscOptionsObject,"Walker-Pernice options");
136: PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,NULL);
137: PetscOptionsTail();
138: return(0);
139: }
143: /*
144: MatMFFDDestroy_WP - Frees the space allocated by
145: MatCreateMFFD_WP().
147: Input Parameter:
148: . ctx - the matrix free context
150: Notes: does not free the ctx, that is handled by the calling routine
152: */
153: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
154: {
158: PetscFree(ctx->hctx);
159: return(0);
160: }
164: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)
165: {
166: MatMFFD ctx = (MatMFFD)mat->data;
167: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
170: hctx->computenormU = flag;
171: return(0);
172: }
176: /*@
177: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
178: PETSc routine for computing h. With any Krylov solver this need only
179: be computed during the first iteration and kept for later.
181: Input Parameters:
182: + A - the matrix created with MatCreateSNESMF()
183: - flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value
185: Options Database Key:
186: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
187: must be sure that ||U|| has not changed in the mean time.
189: Level: advanced
191: Notes:
192: See the manual page for MATMFFD_WP for a complete description of the
193: algorithm used to compute h.
195: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
197: @*/
198: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)
199: {
204: PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
205: return(0);
206: }
210: /*
211: MatCreateMFFD_WP - Standard PETSc code for
212: computing h with matrix-free finite differences.
214: Input Parameter:
215: . ctx - the matrix free context created by MatCreateMFFD()
217: */
218: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
219: {
221: MatMFFD_WP *hctx;
224: /* allocate my own private data structure */
225: PetscNewLog(ctx,&hctx);
226: ctx->hctx = (void*)hctx;
227: hctx->computenormU = PETSC_FALSE;
229: /* set the functions I am providing */
230: ctx->ops->compute = MatMFFDCompute_WP;
231: ctx->ops->destroy = MatMFFDDestroy_WP;
232: ctx->ops->view = MatMFFDView_WP;
233: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
235: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C",MatMFFDWPSetComputeNormU_P);
236: return(0);
237: }