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