Actual source code: snesob.c

      1: #include <petsc/private/snesimpl.h>

      3: /*@C
      4:   SNESSetObjective - Sets the objective function minimized by some of the `SNES` linesearch methods, used instead of the 2-norm of the residual in the line search

      6:   Logically Collective

      8:   Input Parameters:
      9: + snes - the `SNES` context
     10: . obj  - objective evaluation routine; see `SNESObjectiveFn` for the calling sequence
     11: - ctx  - [optional] user-defined context for private data for the objective evaluation routine (may be `NULL`)

     13:   Level: intermediate

     15:   Note:
     16:   Some of the `SNESLineSearch` methods attempt to minimize a given objective provided by this function to determine a step length.

     18:   If not provided then this defaults to the two-norm of the function evaluation (set with `SNESSetFunction()`)

     20:   This is not used in the `SNESLINESEARCHCP` line search.

     22: .seealso: [](ch_snes), `SNES`, `SNESLineSearch()`, `SNESGetObjective()`, `SNESComputeObjective()`, `SNESSetFunction()`, `SNESSetJacobian()`,
     23:           `SNESObjectiveFn`, `SNESSetObjectiveDomainError()`
     24: @*/
     25: PetscErrorCode SNESSetObjective(SNES snes, SNESObjectiveFn *obj, PetscCtx ctx)
     26: {
     27:   DM dm;

     29:   PetscFunctionBegin;
     31:   PetscCall(SNESGetDM(snes, &dm));
     32:   PetscCall(DMSNESSetObjective(dm, obj, ctx));
     33:   PetscFunctionReturn(PETSC_SUCCESS);
     34: }

     36: /*@C
     37:   SNESGetObjective - Returns the objective function set with `SNESSetObjective()`

     39:   Not Collective

     41:   Input Parameter:
     42: . snes - the `SNES` context

     44:   Output Parameters:
     45: + obj - objective evaluation routine (or `NULL`); see `SNESObjectiveFn` for the calling sequence
     46: - ctx - the function context (or `NULL`)

     48:   Level: advanced

     50: .seealso: [](ch_snes), `SNES`, `SNESSetObjective()`, `SNESGetSolution()`, `SNESObjectiveFn`
     51: @*/
     52: PetscErrorCode SNESGetObjective(SNES snes, SNESObjectiveFn **obj, PetscCtxRt ctx)
     53: {
     54:   DM dm;

     56:   PetscFunctionBegin;
     58:   PetscCall(SNESGetDM(snes, &dm));
     59:   PetscCall(DMSNESGetObjective(dm, obj, ctx));
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: /*@
     64:   SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()`

     66:   Collective

     68:   Input Parameters:
     69: + snes - the `SNES` context
     70: - X    - the state vector

     72:   Output Parameter:
     73: . ob - the objective value

     75:   Level: developer

     77:   Notes:
     78:   `SNESComputeObjective()` is typically used within line-search routines,
     79:   so users would not generally call this routine themselves.

     81:   When solving for $F(x) = b$, this routine computes $objective(x) - x^T b$ where $objective(x)$ is the function provided with `SNESSetObjective()`

     83: .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()`
     84: @*/
     85: PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob)
     86: {
     87:   DM     dm;
     88:   DMSNES sdm;

     90:   PetscFunctionBegin;
     93:   PetscAssertPointer(ob, 3);
     94:   PetscCall(SNESGetDM(snes, &dm));
     95:   PetscCall(DMGetDMSNES(dm, &sdm));
     96:   PetscCheck(sdm->ops->computeobjective, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
     97:   PetscCall(PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0));
     98:   PetscCall(sdm->ops->computeobjective(snes, X, ob, sdm->objectivectx));
     99:   PetscCall(PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0));
    100:   if (snes->vec_rhs) {
    101:     PetscScalar dot;

    103:     PetscCall(VecDot(snes->vec_rhs, X, &dot));
    104:     *ob -= PetscRealPart(dot);
    105:   }
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: /*@C
    110:   SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function

    112:   Collective

    114:   Input Parameters:
    115: + snes - the `SNES` context
    116: . X    - the state vector
    117: - ctx  - the (ignored) function context

    119:   Output Parameter:
    120: . F - the function value

    122:   Options Database Keys:
    123: + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6
    124: - -snes_fd_function     - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference

    126:   Level: advanced

    128:   Notes:
    129:   This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function
    130:   `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`.
    131:   Therefore, it should be used for debugging purposes only.  Using it in conjunction with
    132:   `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite
    133:   noisy.  This is often necessary, but should be done with care, even when debugging
    134:   small problems.

    136:   This uses quadratic interpolation of the objective to form each value in the function.

    138: .seealso: [](ch_snes), `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()`, `SNESObjectiveFn`
    139: @*/
    140: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, PetscCtx ctx)
    141: {
    142:   Vec         Xh;
    143:   PetscInt    i, N, start, end;
    144:   PetscReal   ob, ob1, ob2, ob3, fob, dx, eps = 1e-6;
    145:   PetscScalar fv, xv;

    147:   PetscFunctionBegin;
    148:   PetscCall(VecDuplicate(X, &Xh));
    149:   PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES");
    150:   PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL));
    151:   PetscOptionsEnd();
    152:   PetscCall(VecSet(F, 0.));

    154:   PetscCall(VecNorm(X, NORM_2, &fob));

    156:   PetscCall(VecGetSize(X, &N));
    157:   PetscCall(VecGetOwnershipRange(X, &start, &end));
    158:   PetscCall(SNESComputeObjective(snes, X, &ob));

    160:   if (fob > 0.) dx = 1e-6 * fob;
    161:   else dx = 1e-6;

    163:   for (i = 0; i < N; i++) {
    164:     /* compute the 1st value */
    165:     PetscCall(VecCopy(X, Xh));
    166:     if (i >= start && i < end) {
    167:       xv = dx;
    168:       PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
    169:     }
    170:     PetscCall(VecAssemblyBegin(Xh));
    171:     PetscCall(VecAssemblyEnd(Xh));
    172:     PetscCall(SNESComputeObjective(snes, Xh, &ob1));

    174:     /* compute the 2nd value */
    175:     PetscCall(VecCopy(X, Xh));
    176:     if (i >= start && i < end) {
    177:       xv = 2. * dx;
    178:       PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
    179:     }
    180:     PetscCall(VecAssemblyBegin(Xh));
    181:     PetscCall(VecAssemblyEnd(Xh));
    182:     PetscCall(SNESComputeObjective(snes, Xh, &ob2));

    184:     /* compute the 3rd value */
    185:     PetscCall(VecCopy(X, Xh));
    186:     if (i >= start && i < end) {
    187:       xv = -dx;
    188:       PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES));
    189:     }
    190:     PetscCall(VecAssemblyBegin(Xh));
    191:     PetscCall(VecAssemblyEnd(Xh));
    192:     PetscCall(SNESComputeObjective(snes, Xh, &ob3));

    194:     if (i >= start && i < end) {
    195:       /* set this entry to be the gradient of the objective */
    196:       fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx);
    197:       if (PetscAbsScalar(fv) > eps) {
    198:         PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
    199:       } else {
    200:         fv = 0.;
    201:         PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES));
    202:       }
    203:     }
    204:   }
    205:   PetscCall(VecDestroy(&Xh));
    206:   PetscCall(VecAssemblyBegin(F));
    207:   PetscCall(VecAssemblyEnd(F));
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }