Actual source code: snesob.c
1: #include <petsc/private/snesimpl.h>
3: /*MC
4: SNESObjectiveFunction - functional form used to convey an objective function to the nonlinear solver
6: Synopsis:
7: #include <petscsnes.h>
8: SNESObjectiveFunction(SNES snes,Vec x,PetscReal *obj,void *ctx);
10: Input Parameters:
11: + snes - the SNES context
12: . X - solution
13: . obj - real to hold the objective value
14: - ctx - optional user-defined objective context
16: Level: advanced
18: .seealso: `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESSetObjective()`, `SNESGetObjective()`, `SNESJacobianFunction`, `SNESFunction`
19: M*/
21: /*@C
22: SNESSetObjective - Sets the objective function minimized by some of the `SNES` linesearch methods.
24: Logically Collective
26: Input Parameters:
27: + snes - the `SNES` context
28: . obj - objective evaluation routine; see `SNESObjectiveFunction` for details
29: - ctx - [optional] user-defined context for private data for the
30: function evaluation routine (may be NULL)
32: Level: intermediate
34: Note:
35: Some of the `SNESLineSearch` methods attempt to minimize a given objective provided by this function to determine a step length.
37: If not provided then this defaults to the two norm of the function evaluation (set with `SNESSetFunction()`)
39: This is not used in the `SNESLINESEARCHCP` line search.
41: .seealso: `SNES`, `SNESLineSearch()`, `SNESGetObjective()`, `SNESComputeObjective()`, `SNESSetFunction()`, `SNESSetJacobian()`, `SNESObjectiveFunction`
42: @*/
43: PetscErrorCode SNESSetObjective(SNES snes, PetscErrorCode (*obj)(SNES, Vec, PetscReal *, void *), void *ctx)
44: {
45: DM dm;
48: SNESGetDM(snes, &dm);
49: DMSNESSetObjective(dm, obj, ctx);
50: return 0;
51: }
53: /*@C
54: SNESGetObjective - Returns the objective function.
56: Not Collective
58: Input Parameter:
59: . snes - the `SNES` context
61: Output Parameters:
62: + obj - objective evaluation routine (or NULL); see `SNESObjectFunction` for details
63: - ctx - the function context (or NULL)
65: Level: advanced
67: .seealso: `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
68: @*/
69: PetscErrorCode SNESGetObjective(SNES snes, PetscErrorCode (**obj)(SNES, Vec, PetscReal *, void *), void **ctx)
70: {
71: DM dm;
74: SNESGetDM(snes, &dm);
75: DMSNESGetObjective(dm, obj, ctx);
76: return 0;
77: }
79: /*@C
80: SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()`
82: Collective
84: Input Parameters:
85: + snes - the `SNES` context
86: - X - the state vector
88: Output Parameter:
89: . ob - the objective value
91: Level: developer
93: .seealso: `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
94: @*/
95: PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob)
96: {
97: DM dm;
98: DMSNES sdm;
103: SNESGetDM(snes, &dm);
104: DMGetDMSNES(dm, &sdm);
105: if (sdm->ops->computeobjective) {
106: PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0);
107: (sdm->ops->computeobjective)(snes, X, ob, sdm->objectivectx);
108: PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0);
109: } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
110: return 0;
111: }
113: /*@C
114: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function
116: Collective
118: Input Parameters:
119: + snes - the `SNES` context
120: . X - the state vector
121: - ctx - the (ignored) function context
123: Output Parameter:
124: . F - the function value
126: Options Database Keys:
127: + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6
128: - -snes_fd_function - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference
130: Notes:
131: This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function
132: `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`.
133: Therefore, it should be used for debugging purposes only. Using it in conjunction with
134: `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite
135: noisy. This is often necessary, but should be done with care, even when debugging
136: small problems.
138: Note that this uses quadratic interpolation of the objective to form each value in the function.
140: Level: advanced
142: .seealso: `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()`
143: @*/
144: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *ctx)
145: {
146: Vec Xh;
147: PetscInt i, N, start, end;
148: PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6;
149: PetscScalar fv, xv;
151: VecDuplicate(X, &Xh);
152: PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES");
153: PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL);
154: PetscOptionsEnd();
155: VecSet(F, 0.);
157: VecNorm(X, NORM_2, &fob);
159: VecGetSize(X, &N);
160: VecGetOwnershipRange(X, &start, &end);
161: SNESComputeObjective(snes, X, &ob);
163: if (fob > 0.) dx = 1e-6 * fob;
164: else dx = 1e-6;
166: for (i = 0; i < N; i++) {
167: /* compute the 1st value */
168: VecCopy(X, Xh);
169: if (i >= start && i < end) {
170: xv = dx;
171: VecSetValues(Xh, 1, &i, &xv, ADD_VALUES);
172: }
173: VecAssemblyBegin(Xh);
174: VecAssemblyEnd(Xh);
175: SNESComputeObjective(snes, Xh, &ob1);
177: /* compute the 2nd value */
178: VecCopy(X, Xh);
179: if (i >= start && i < end) {
180: xv = 2. * dx;
181: VecSetValues(Xh, 1, &i, &xv, ADD_VALUES);
182: }
183: VecAssemblyBegin(Xh);
184: VecAssemblyEnd(Xh);
185: SNESComputeObjective(snes, Xh, &ob2);
187: /* compute the 3rd value */
188: VecCopy(X, Xh);
189: if (i >= start && i < end) {
190: xv = -dx;
191: VecSetValues(Xh, 1, &i, &xv, ADD_VALUES);
192: }
193: VecAssemblyBegin(Xh);
194: VecAssemblyEnd(Xh);
195: SNESComputeObjective(snes, Xh, &ob3);
197: if (i >= start && i < end) {
198: /* set this entry to be the gradient of the objective */
199: fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx);
200: if (PetscAbsScalar(fv) > eps) {
201: VecSetValues(F, 1, &i, &fv, INSERT_VALUES);
202: } else {
203: fv = 0.;
204: VecSetValues(F, 1, &i, &fv, INSERT_VALUES);
205: }
206: }
207: }
208: VecDestroy(&Xh);
210: VecAssemblyBegin(F);
211: VecAssemblyEnd(F);
212: return 0;
213: }