Actual source code: wp.c
petsc-3.9.4 2018-09-11
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>
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;
68: if (!(ctx->count % ctx->recomputeperiod)) {
69: if (hctx->computenormU || !ctx->ncurrenth) {
70: VecNorm(U,NORM_2,&normU);
71: hctx->normUfact = PetscSqrtReal(1.0+normU);
72: }
73: VecNorm(a,NORM_2,&norma);
74: if (norma == 0.0) {
75: *zeroa = PETSC_TRUE;
76: return(0);
77: }
78: *zeroa = PETSC_FALSE;
79: *h = ctx->error_rel*hctx->normUfact/norma;
80: } else {
81: *h = ctx->currenth;
82: }
83: return(0);
84: }
86: /*
87: MatMFFDView_WP - Prints information about this particular
88: method for computing h. Note that this does not print the general
89: information about the matrix free, that is printed by the calling
90: routine.
92: Input Parameters:
93: + ctx - the matrix free context
94: - viewer - the PETSc viewer
96: */
97: static PetscErrorCode MatMFFDView_WP(MatMFFD ctx,PetscViewer viewer)
98: {
99: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
101: PetscBool iascii;
104: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
105: if (iascii) {
106: if (hctx->computenormU) {
107: PetscViewerASCIIPrintf(viewer," Computes normU\n");
108: } else {
109: PetscViewerASCIIPrintf(viewer," Does not compute normU\n");
110: }
111: }
112: return(0);
113: }
115: /*
116: MatMFFDSetFromOptions_WP - Looks in the options database for
117: any options appropriate for this method
119: Input Parameter:
120: . ctx - the matrix free context
122: */
123: static PetscErrorCode MatMFFDSetFromOptions_WP(PetscOptionItems *PetscOptionsObject,MatMFFD ctx)
124: {
126: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
129: PetscOptionsHead(PetscOptionsObject,"Walker-Pernice options");
130: PetscOptionsBool("-mat_mffd_compute_normu","Compute the norm of u","MatMFFDWPSetComputeNormU", hctx->computenormU,&hctx->computenormU,NULL);
131: PetscOptionsTail();
132: return(0);
133: }
135: /*
136: MatMFFDDestroy_WP - Frees the space allocated by
137: MatCreateMFFD_WP().
139: Input Parameter:
140: . ctx - the matrix free context
142: Notes: does not free the ctx, that is handled by the calling routine
144: */
145: static PetscErrorCode MatMFFDDestroy_WP(MatMFFD ctx)
146: {
150: PetscFree(ctx->hctx);
151: return(0);
152: }
154: PetscErrorCode MatMFFDWPSetComputeNormU_P(Mat mat,PetscBool flag)
155: {
156: MatMFFD ctx = (MatMFFD)mat->data;
157: MatMFFD_WP *hctx = (MatMFFD_WP*)ctx->hctx;
160: hctx->computenormU = flag;
161: return(0);
162: }
164: /*@
165: MatMFFDWPSetComputeNormU - Sets whether it computes the ||U|| used by the WP
166: PETSc routine for computing h. With any Krylov solver this need only
167: be computed during the first iteration and kept for later.
169: Input Parameters:
170: + A - the matrix created with MatCreateSNESMF()
171: - flag - PETSC_TRUE causes it to compute ||U||, PETSC_FALSE uses the previous value
173: Options Database Key:
174: . -mat_mffd_compute_normu <true,false> - true by default, false can save calculations but you
175: must be sure that ||U|| has not changed in the mean time.
177: Level: advanced
179: Notes:
180: See the manual page for MATMFFD_WP for a complete description of the
181: algorithm used to compute h.
183: .seealso: MatMFFDSetFunctionError(), MatCreateSNESMF()
185: @*/
186: PetscErrorCode MatMFFDWPSetComputeNormU(Mat A,PetscBool flag)
187: {
192: PetscTryMethod(A,"MatMFFDWPSetComputeNormU_C",(Mat,PetscBool),(A,flag));
193: return(0);
194: }
196: /*
197: MatCreateMFFD_WP - Standard PETSc code for
198: computing h with matrix-free finite differences.
200: Input Parameter:
201: . ctx - the matrix free context created by MatCreateMFFD()
203: */
204: PETSC_EXTERN PetscErrorCode MatCreateMFFD_WP(MatMFFD ctx)
205: {
207: MatMFFD_WP *hctx;
210: /* allocate my own private data structure */
211: PetscNewLog(ctx,&hctx);
212: ctx->hctx = (void*)hctx;
213: hctx->computenormU = PETSC_FALSE;
215: /* set the functions I am providing */
216: ctx->ops->compute = MatMFFDCompute_WP;
217: ctx->ops->destroy = MatMFFDDestroy_WP;
218: ctx->ops->view = MatMFFDView_WP;
219: ctx->ops->setfromoptions = MatMFFDSetFromOptions_WP;
221: PetscObjectComposeFunction((PetscObject)ctx->mat,"MatMFFDWPSetComputeNormU_C",MatMFFDWPSetComputeNormU_P);
222: return(0);
223: }