1 #include <petsc/private/snesimpl.h> 2 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 5 6 Synopsis: 7 #include <petscsnes.h> 8 SNESObjectiveFunction(SNES snes,Vec x,PetscReal *obj,void *ctx); 9 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 15 16 Level: advanced 17 18 .seealso: `SNES`, `SNESSetFunction()`, `SNESGetFunction()`, `SNESSetObjective()`, `SNESGetObjective()`, `SNESJacobianFunction`, `SNESFunction` 19 M*/ 20 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 23 24 Logically Collective 25 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`) 31 32 Level: intermediate 33 34 Note: 35 Some of the `SNESLineSearch` methods attempt to minimize a given objective provided by this function to determine a step length. 36 37 If not provided then this defaults to the two norm of the function evaluation (set with `SNESSetFunction()`) 38 39 This is not used in the `SNESLINESEARCHCP` line search. 40 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; 46 47 PetscFunctionBegin; 48 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 49 PetscCall(SNESGetDM(snes, &dm)); 50 PetscCall(DMSNESSetObjective(dm, obj, ctx)); 51 PetscFunctionReturn(PETSC_SUCCESS); 52 } 53 54 /*@C 55 SNESGetObjective - Returns the objective function set with `SNESSetObjective()` 56 57 Not Collective 58 59 Input Parameter: 60 . snes - the `SNES` context 61 62 Output Parameters: 63 + obj - objective evaluation routine (or `NULL`); see `SNESObjectFunction` for details 64 - ctx - the function context (or `NULL`) 65 66 Level: advanced 67 68 .seealso: `SNES`, `SNESSetObjective()`, `SNESGetSolution()` 69 @*/ 70 PetscErrorCode SNESGetObjective(SNES snes, PetscErrorCode (**obj)(SNES, Vec, PetscReal *, void *), void **ctx) 71 { 72 DM dm; 73 74 PetscFunctionBegin; 75 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 76 PetscCall(SNESGetDM(snes, &dm)); 77 PetscCall(DMSNESGetObjective(dm, obj, ctx)); 78 PetscFunctionReturn(PETSC_SUCCESS); 79 } 80 81 /*@C 82 SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()` 83 84 Collective 85 86 Input Parameters: 87 + snes - the `SNES` context 88 - X - the state vector 89 90 Output Parameter: 91 . ob - the objective value 92 93 Level: developer 94 95 .seealso: `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()` 96 @*/ 97 PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob) 98 { 99 DM dm; 100 DMSNES sdm; 101 102 PetscFunctionBegin; 103 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 104 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 105 PetscValidRealPointer(ob, 3); 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 } 114 115 /*@C 116 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function 117 118 Collective 119 120 Input Parameters: 121 + snes - the `SNES` context 122 . X - the state vector 123 - ctx - the (ignored) function context 124 125 Output Parameter: 126 . F - the function value 127 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 131 132 Level: advanced 133 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. 141 142 Note that this uses quadratic interpolation of the objective to form each value in the function. 143 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; 152 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.)); 159 160 PetscCall(VecNorm(X, NORM_2, &fob)); 161 162 PetscCall(VecGetSize(X, &N)); 163 PetscCall(VecGetOwnershipRange(X, &start, &end)); 164 PetscCall(SNESComputeObjective(snes, X, &ob)); 165 166 if (fob > 0.) dx = 1e-6 * fob; 167 else dx = 1e-6; 168 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)); 179 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)); 189 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)); 199 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)); 212 213 PetscCall(VecAssemblyBegin(F)); 214 PetscCall(VecAssemblyEnd(F)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217