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: [](ch_snes), `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: [](ch_snes), `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 `SNESObjectiveFunction` for details 64 - ctx - the function context (or `NULL`) 65 66 Level: advanced 67 68 .seealso: [](ch_snes), `SNES`, `SNESSetObjective()`, `SNESObjectiveFunction`, `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 Notes: 96 `SNESComputeObjective()` is typically used within line-search routines, 97 so users would not generally call this routine themselves. 98 99 When solving for $F(x) = b$, this routine computes $objective(x) - x^T b$ where $objective(x)$ is the function provided with `SNESSetObjective()` 100 101 .seealso: [](ch_snes), `SNESLineSearch`, `SNES`, `SNESSetObjective()`, `SNESGetSolution()` 102 @*/ 103 PetscErrorCode SNESComputeObjective(SNES snes, Vec X, PetscReal *ob) 104 { 105 DM dm; 106 DMSNES sdm; 107 108 PetscFunctionBegin; 109 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 110 PetscValidHeaderSpecific(X, VEC_CLASSID, 2); 111 PetscAssertPointer(ob, 3); 112 PetscCall(SNESGetDM(snes, &dm)); 113 PetscCall(DMGetDMSNES(dm, &sdm)); 114 PetscCheck(sdm->ops->computeobjective, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 115 PetscCall(PetscLogEventBegin(SNES_ObjectiveEval, snes, X, 0, 0)); 116 PetscCall((sdm->ops->computeobjective)(snes, X, ob, sdm->objectivectx)); 117 PetscCall(PetscLogEventEnd(SNES_ObjectiveEval, snes, X, 0, 0)); 118 if (snes->vec_rhs) { 119 PetscScalar dot; 120 121 PetscCall(VecDot(snes->vec_rhs, X, &dot)); 122 *ob -= PetscRealPart(dot); 123 } 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 /*@C 128 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective function 129 130 Collective 131 132 Input Parameters: 133 + snes - the `SNES` context 134 . X - the state vector 135 - ctx - the (ignored) function context 136 137 Output Parameter: 138 . F - the function value 139 140 Options Database Keys: 141 + -snes_fd_function_eps - Tolerance for including non-zero entries into the gradient, default is 1.e-6 142 - -snes_fd_function - Computes function from user provided objective function (set with `SNESSetObjective()`) with finite difference 143 144 Level: advanced 145 146 Notes: 147 This function can be used with `SNESSetFunction()` to have the nonlinear function solved for with `SNES` defined by the gradient of an objective function 148 `SNESObjectiveComputeFunctionDefaultFD()` is similar in character to `SNESComputeJacobianDefault()`. 149 Therefore, it should be used for debugging purposes only. Using it in conjunction with 150 `SNESComputeJacobianDefault()` is excessively costly and produces a Jacobian that is quite 151 noisy. This is often necessary, but should be done with care, even when debugging 152 small problems. 153 154 This uses quadratic interpolation of the objective to form each value in the function. 155 156 .seealso: [](ch_snes), `SNESSetObjective()`, `SNESSetFunction()`, `SNESComputeObjective()`, `SNESComputeJacobianDefault()` 157 @*/ 158 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes, Vec X, Vec F, void *ctx) 159 { 160 Vec Xh; 161 PetscInt i, N, start, end; 162 PetscReal ob, ob1, ob2, ob3, fob, dx, eps = 1e-6; 163 PetscScalar fv, xv; 164 165 PetscFunctionBegin; 166 PetscCall(VecDuplicate(X, &Xh)); 167 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing parameters", "SNES"); 168 PetscCall(PetscOptionsReal("-snes_fd_function_eps", "Tolerance for nonzero entries in fd function", "None", eps, &eps, NULL)); 169 PetscOptionsEnd(); 170 PetscCall(VecSet(F, 0.)); 171 172 PetscCall(VecNorm(X, NORM_2, &fob)); 173 174 PetscCall(VecGetSize(X, &N)); 175 PetscCall(VecGetOwnershipRange(X, &start, &end)); 176 PetscCall(SNESComputeObjective(snes, X, &ob)); 177 178 if (fob > 0.) dx = 1e-6 * fob; 179 else dx = 1e-6; 180 181 for (i = 0; i < N; i++) { 182 /* compute the 1st value */ 183 PetscCall(VecCopy(X, Xh)); 184 if (i >= start && i < end) { 185 xv = dx; 186 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 187 } 188 PetscCall(VecAssemblyBegin(Xh)); 189 PetscCall(VecAssemblyEnd(Xh)); 190 PetscCall(SNESComputeObjective(snes, Xh, &ob1)); 191 192 /* compute the 2nd value */ 193 PetscCall(VecCopy(X, Xh)); 194 if (i >= start && i < end) { 195 xv = 2. * dx; 196 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 197 } 198 PetscCall(VecAssemblyBegin(Xh)); 199 PetscCall(VecAssemblyEnd(Xh)); 200 PetscCall(SNESComputeObjective(snes, Xh, &ob2)); 201 202 /* compute the 3rd value */ 203 PetscCall(VecCopy(X, Xh)); 204 if (i >= start && i < end) { 205 xv = -dx; 206 PetscCall(VecSetValues(Xh, 1, &i, &xv, ADD_VALUES)); 207 } 208 PetscCall(VecAssemblyBegin(Xh)); 209 PetscCall(VecAssemblyEnd(Xh)); 210 PetscCall(SNESComputeObjective(snes, Xh, &ob3)); 211 212 if (i >= start && i < end) { 213 /* set this entry to be the gradient of the objective */ 214 fv = (-ob2 + 6. * ob1 - 3. * ob - 2. * ob3) / (6. * dx); 215 if (PetscAbsScalar(fv) > eps) { 216 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 217 } else { 218 fv = 0.; 219 PetscCall(VecSetValues(F, 1, &i, &fv, INSERT_VALUES)); 220 } 221 } 222 } 223 PetscCall(VecDestroy(&Xh)); 224 PetscCall(VecAssemblyBegin(F)); 225 PetscCall(VecAssemblyEnd(F)); 226 PetscFunctionReturn(PETSC_SUCCESS); 227 } 228