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: }