1 #include <petsc/private/snesimpl.h> 2 3 /*MC 4 SNESObjectiveFunction - functional form used to convey an objective function to the nonlinear solver 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. 23 24 Logically Collective on snes 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(0); 52 } 53 54 /*@C 55 SNESGetObjective - Returns the objective function. 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(0); 79 } 80 81 /*@C 82 SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()` 83 84 Collective on snes 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 if (sdm->ops->computeobjective) { 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 } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 113 PetscFunctionReturn(0); 114 } 115 116 /*@C 117 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function 118 119 Collective on snes 120 121 Input Parameters: 122 + snes - the `SNES` context 123 . X - the state vector 124 - ctx - the (ignored) function context 125 126 Output Parameter: 127 . F - the function value 128 129 Options Database Keys: 130 + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6 131 - -snes_fd_function - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference 132 133 Notes: 134 This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function 135 `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`. 136 Therefore, it should be used for debugging purposes only. Using it in conjunction with 137 `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite 138 noisy. This is often necessary, but should be done with care, even when debugging 139 small problems. 140 141 Note that this uses quadratic interpolation of the objective to form each value in the function. 142 143 Level: advanced 144 145 .seealso: `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()` 146 @*/ 147 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *ctx) 148 { 149 Vec Xh; 150 PetscInt i, N, start, end; 151 PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6; 152 PetscScalar fv, xv; 153 154 PetscFunctionBegin; 155 PetscCall(VecDuplicate(X, &Xh)); 156 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES"); 157 PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL)); 158 PetscOptionsEnd(); 159 PetscCall(VecSet(F, 0.)); 160 161 PetscCall(VecNorm(X, NORM_2, &fob)); 162 163 PetscCall(VecGetSize(X, &N)); 164 PetscCall(VecGetOwnershipRange(X, &start, &end)); 165 PetscCall(SNESComputeObjective(snes, X, &ob)); 166 167 if (fob > 0.) dx = 1e-6 * fob; 168 else dx = 1e-6; 169 170 for (i = 0; i < N; i++) { 171 /* compute the 1st value */ 172 PetscCall(VecCopy(X, Xh)); 173 if (i >= start && i < end) { 174 xv = dx; 175 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 176 } 177 PetscCall(VecAssemblyBegin(Xh)); 178 PetscCall(VecAssemblyEnd(Xh)); 179 PetscCall(SNESComputeObjective(snes, Xh, &ob1)); 180 181 /* compute the 2nd value */ 182 PetscCall(VecCopy(X, Xh)); 183 if (i >= start && i < end) { 184 xv = 2. * dx; 185 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 186 } 187 PetscCall(VecAssemblyBegin(Xh)); 188 PetscCall(VecAssemblyEnd(Xh)); 189 PetscCall(SNESComputeObjective(snes, Xh, &ob2)); 190 191 /* compute the 3rd value */ 192 PetscCall(VecCopy(X, Xh)); 193 if (i >= start && i < end) { 194 xv = -dx; 195 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 196 } 197 PetscCall(VecAssemblyBegin(Xh)); 198 PetscCall(VecAssemblyEnd(Xh)); 199 PetscCall(SNESComputeObjective(snes, Xh, &ob3)); 200 201 if (i >= start && i < end) { 202 /* set this entry to be the gradient of the objective */ 203 fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx); 204 if (PetscAbsScalar(fv) > eps) { 205 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 206 } else { 207 fv = 0.; 208 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 209 } 210 } 211 } 212 PetscCall(VecDestroy(&Xh)); 213 214 PetscCall(VecAssemblyBegin(F)); 215 PetscCall(VecAssemblyEnd(F)); 216 PetscFunctionReturn(0); 217 } 218