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, that will be used instead of the 2-norm of the residual
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, used instead of the 2-norm of the residual
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;
47: PetscFunctionBegin;
49: PetscCall(SNESGetDM(snes, &dm));
50: PetscCall(DMSNESSetObjective(dm, obj, ctx));
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@C
55: SNESGetObjective - Returns the objective function set with `SNESSetObjective()`
57: Not Collective
59: Input Parameter:
60: . snes - the `SNES` context
62: Output Parameters:
63: + obj - objective evaluation routine (or `NULL`); see `SNESObjectFunction` for details
64: - ctx - the function context (or `NULL`)
66: Level: advanced
68: .seealso: `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
69: @*/
70: PetscErrorCode SNESGetObjective(SNES snes, PetscErrorCode (**obj)(SNES, Vec, PetscReal *, void *), void **ctx)
71: {
72: DM dm;
74: PetscFunctionBegin;
76: PetscCall(SNESGetDM(snes, &dm));
77: PetscCall(DMSNESGetObjective(dm, obj, ctx));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*@C
82: SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()`
84: Collective
86: Input Parameters:
87: + snes - the `SNES` context
88: - X - the state vector
90: Output Parameter:
91: . ob - the objective value
93: Level: developer
95: .seealso: `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
96: @*/
97: PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob)
98: {
99: DM dm;
100: DMSNES sdm;
102: PetscFunctionBegin;
106: PetscCall(SNESGetDM(snes, &dm));
107: PetscCall(DMGetDMSNES(dm, &sdm));
108: PetscCheck(sdm->ops->computeobjective, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
109: PetscCall(PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0));
110: PetscCall((sdm->ops->computeobjective)(snes, X, ob, sdm->objectivectx));
111: PetscCall(PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: /*@C
116: SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function
118: Collective
120: Input Parameters:
121: + snes - the `SNES` context
122: . X - the state vector
123: - ctx - the (ignored) function context
125: Output Parameter:
126: . F - the function value
128: Options Database Keys:
129: + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6
130: - -snes_fd_function - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference
132: Level: advanced
134: Notes:
135: This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function
136: `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`.
137: Therefore, it should be used for debugging purposes only. Using it in conjunction with
138: `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite
139: noisy. This is often necessary, but should be done with care, even when debugging
140: small problems.
142: Note that this uses quadratic interpolation of the objective to form each value in the function.
144: .seealso: `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()`
145: @*/
146: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *ctx)
147: {
148: Vec Xh;
149: PetscInt i, N, start, end;
150: PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6;
151: PetscScalar fv, xv;
153: PetscFunctionBegin;
154: PetscCall(VecDuplicate(X, &Xh));
155: PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES");
156: PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL));
157: PetscOptionsEnd();
158: PetscCall(VecSet(F, 0.));
160: PetscCall(VecNorm(X, NORM_2, &fob));
162: PetscCall(VecGetSize(X, &N));
163: PetscCall(VecGetOwnershipRange(X, &start, &end));
164: PetscCall(SNESComputeObjective(snes, X, &ob));
166: if (fob > 0.) dx = 1e-6 * fob;
167: else dx = 1e-6;
169: for (i = 0; i < N; i++) {
170: /* compute the 1st value */
171: PetscCall(VecCopy(X, Xh));
172: if (i >= start && i < end) {
173: xv = dx;
174: PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
175: }
176: PetscCall(VecAssemblyBegin(Xh));
177: PetscCall(VecAssemblyEnd(Xh));
178: PetscCall(SNESComputeObjective(snes, Xh, &ob1));
180: /* compute the 2nd value */
181: PetscCall(VecCopy(X, Xh));
182: if (i >= start && i < end) {
183: xv = 2. * dx;
184: PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
185: }
186: PetscCall(VecAssemblyBegin(Xh));
187: PetscCall(VecAssemblyEnd(Xh));
188: PetscCall(SNESComputeObjective(snes, Xh, &ob2));
190: /* compute the 3rd value */
191: PetscCall(VecCopy(X, Xh));
192: if (i >= start && i < end) {
193: xv = -dx;
194: PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
195: }
196: PetscCall(VecAssemblyBegin(Xh));
197: PetscCall(VecAssemblyEnd(Xh));
198: PetscCall(SNESComputeObjective(snes, Xh, &ob3));
200: if (i >= start && i < end) {
201: /* set this entry to be the gradient of the objective */
202: fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx);
203: if (PetscAbsScalar(fv) > eps) {
204: PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
205: } else {
206: fv = 0.;
207: PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
208: }
209: }
210: }
211: PetscCall(VecDestroy(&Xh));
213: PetscCall(VecAssemblyBegin(F));
214: PetscCall(VecAssemblyEnd(F));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }