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 DM dm; 45 46 PetscFunctionBegin; 47 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 48 PetscCall(SNESGetDM(snes, &dm)); 49 PetscCall(DMSNESSetObjective(dm, obj, ctx)); 50 PetscFunctionReturn(0); 51 } 52 53 /*@C 54 SNESGetObjective - Returns the objective function. 55 56 Not Collective 57 58 Input Parameter: 59 . snes - the `SNES` context 60 61 Output Parameters: 62 + obj - objective evaluation routine (or NULL); see `SNESObjectFunction` for details 63 - ctx - the function context (or NULL) 64 65 Level: advanced 66 67 .seealso: `SNES`, `SNESSetObjective()`, `SNESGetSolution()` 68 @*/ 69 PetscErrorCode SNESGetObjective(SNES snes, PetscErrorCode (**obj)(SNES, Vec, PetscReal *, void *), void **ctx) { 70 DM dm; 71 72 PetscFunctionBegin; 73 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 74 PetscCall(SNESGetDM(snes, &dm)); 75 PetscCall(DMSNESGetObjective(dm, obj, ctx)); 76 PetscFunctionReturn(0); 77 } 78 79 /*@C 80 SNESComputeObjective - Computes the objective function that has been provided by `SNESSetObjective()` 81 82 Collective on snes 83 84 Input Parameters: 85 + snes - the `SNES` context 86 - X - the state vector 87 88 Output Parameter: 89 . ob - the objective value 90 91 Level: developer 92 93 .seealso: `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()` 94 @*/ 95 PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob) { 96 DM dm; 97 DMSNES sdm; 98 99 PetscFunctionBegin; 100 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 101 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 102 PetscValidRealPointer(ob, 3); 103 PetscCall(SNESGetDM(snes, &dm)); 104 PetscCall(DMGetDMSNES(dm, &sdm)); 105 if (sdm->ops->computeobjective) { 106 PetscCall(PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0)); 107 PetscCall((sdm->ops->computeobjective)(snes, X, ob, sdm->objectivectx)); 108 PetscCall(PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0)); 109 } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 110 PetscFunctionReturn(0); 111 } 112 113 /*@C 114 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function 115 116 Collective on snes 117 118 Input Parameters: 119 + snes - the `SNES` context 120 . X - the state vector 121 - ctx - the (ignored) function context 122 123 Output Parameter: 124 . F - the function value 125 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 129 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. 137 138 Note that this uses quadratic interpolation of the objective to form each value in the function. 139 140 Level: advanced 141 142 .seealso: `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()` 143 @*/ 144 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *ctx) { 145 Vec Xh; 146 PetscInt i, N, start, end; 147 PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6; 148 PetscScalar fv, xv; 149 150 PetscFunctionBegin; 151 PetscCall(VecDuplicate(X, &Xh)); 152 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES"); 153 PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL)); 154 PetscOptionsEnd(); 155 PetscCall(VecSet(F, 0.)); 156 157 PetscCall(VecNorm(X, NORM_2, &fob)); 158 159 PetscCall(VecGetSize(X, &N)); 160 PetscCall(VecGetOwnershipRange(X, &start, &end)); 161 PetscCall(SNESComputeObjective(snes, X, &ob)); 162 163 if (fob > 0.) dx = 1e-6 * fob; 164 else dx = 1e-6; 165 166 for (i = 0; i < N; i++) { 167 /* compute the 1st value */ 168 PetscCall(VecCopy(X, Xh)); 169 if (i >= start && i < end) { 170 xv = dx; 171 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 172 } 173 PetscCall(VecAssemblyBegin(Xh)); 174 PetscCall(VecAssemblyEnd(Xh)); 175 PetscCall(SNESComputeObjective(snes, Xh, &ob1)); 176 177 /* compute the 2nd value */ 178 PetscCall(VecCopy(X, Xh)); 179 if (i >= start && i < end) { 180 xv = 2. * dx; 181 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 182 } 183 PetscCall(VecAssemblyBegin(Xh)); 184 PetscCall(VecAssemblyEnd(Xh)); 185 PetscCall(SNESComputeObjective(snes, Xh, &ob2)); 186 187 /* compute the 3rd value */ 188 PetscCall(VecCopy(X, Xh)); 189 if (i >= start && i < end) { 190 xv = -dx; 191 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 192 } 193 PetscCall(VecAssemblyBegin(Xh)); 194 PetscCall(VecAssemblyEnd(Xh)); 195 PetscCall(SNESComputeObjective(snes, Xh, &ob3)); 196 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 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 202 } else { 203 fv = 0.; 204 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 205 } 206 } 207 } 208 PetscCall(VecDestroy(&Xh)); 209 210 PetscCall(VecAssemblyBegin(F)); 211 PetscCall(VecAssemblyEnd(F)); 212 PetscFunctionReturn(0); 213 } 214