1 #include <petsc/private/snesimpl.h>
2
3 /*MC
4 SNESObjectiveFunction - functional form used to convey the 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: SNESSetFunction(), SNESGetFunction(), SNESSetObjective(), SNESGetObjective(), SNESJacobianFunction, SNESFunction
19 M*/
20
21
22 /*@C
23 SNESSetObjective - Sets the objective function minimized by some of the SNES linesearch methods.
24
25 Logically Collective on SNES
26
27 Input Parameters:
28 + snes - the SNES context
29 . obj - objective evaluation routine; see SNESObjectiveFunction for details
30 - ctx - [optional] user-defined context for private data for the
31 function evaluation routine (may be NULL)
32
33 Level: intermediate
34
35 Note: This is not used in the SNESLINESEARCHCP line search.
36
37 If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction())
38
39 .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction
40 @*/
SNESSetObjective(SNES snes,PetscErrorCode (* obj)(SNES,Vec,PetscReal *,void *),void * ctx)41 PetscErrorCode SNESSetObjective(SNES snes,PetscErrorCode (*obj)(SNES,Vec,PetscReal*,void*),void *ctx)
42 {
43 PetscErrorCode ierr;
44 DM dm;
45
46 PetscFunctionBegin;
47 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
48 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
49 ierr = DMSNESSetObjective(dm,obj,ctx);CHKERRQ(ierr);
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 Parameter:
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: SNESSetObjective(), SNESGetSolution()
68 @*/
SNESGetObjective(SNES snes,PetscErrorCode (** obj)(SNES,Vec,PetscReal *,void *),void ** ctx)69 PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**obj)(SNES,Vec,PetscReal*,void*),void **ctx)
70 {
71 PetscErrorCode ierr;
72 DM dm;
73
74 PetscFunctionBegin;
75 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
76 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
77 ierr = DMSNESGetObjective(dm,obj,ctx);CHKERRQ(ierr);
78 PetscFunctionReturn(0);
79 }
80
81 /*@C
82 SNESComputeObjective - Computes the objective.
83
84 Collective on SNES
85
86 Input Parameter:
87 + snes - the SNES context
88 - X - the state vector
89
90 Output Parameter:
91 . ob - the objective value
92
93 Level: advanced
94
95 .seealso: SNESSetObjective(), SNESGetSolution()
96 @*/
SNESComputeObjective(SNES snes,Vec X,PetscReal * ob)97 PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
98 {
99 PetscErrorCode ierr;
100 DM dm;
101 DMSNES sdm;
102
103 PetscFunctionBegin;
104 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
105 PetscValidHeaderSpecific(X,VEC_CLASSID,2);
106 PetscValidRealPointer(ob,3);
107 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
108 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr);
109 if (sdm->ops->computeobjective) {
110 ierr = PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0);CHKERRQ(ierr);
111 ierr = (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);CHKERRQ(ierr);
112 ierr = PetscLogEventEnd(SNES_ObjectiveEval,snes,X,0,0);CHKERRQ(ierr);
113 } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
114 PetscFunctionReturn(0);
115 }
116
117
118 /*@C
119 SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective
120
121 Collective on SNES
122
123 Input Parameter:
124 + snes - the SNES context
125 . X - the state vector
126 - ctx - the (ignored) function context
127
128 Output Parameter:
129 . F - the function value
130
131 Options Database Key:
132 + -snes_fd_function_eps - The differencing parameter
133 - -snes_fd_function - Compute function from user provided objective with finite difference
134
135 Notes:
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 a grain of salt, 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 Level: advanced
145
146 .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
147 @*/
SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void * ctx)148 PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
149 {
150 Vec Xh;
151 PetscErrorCode ierr;
152 PetscInt i,N,start,end;
153 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
154 PetscScalar fv,xv;
155
156 PetscFunctionBegin;
157 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr);
158 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");CHKERRQ(ierr);
159 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);CHKERRQ(ierr);
160 ierr = PetscOptionsEnd();CHKERRQ(ierr);
161 ierr = VecSet(F,0.);CHKERRQ(ierr);
162
163 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr);
164
165 ierr = VecGetSize(X,&N);CHKERRQ(ierr);
166 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr);
167 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr);
168
169 if (fob > 0.) dx =1e-6*fob;
170 else dx = 1e-6;
171
172 for (i=0; i<N; i++) {
173 /* compute the 1st value */
174 ierr = VecCopy(X,Xh);CHKERRQ(ierr);
175 if (i>= start && i<end) {
176 xv = dx;
177 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr);
178 }
179 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr);
180 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr);
181 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr);
182
183 /* compute the 2nd value */
184 ierr = VecCopy(X,Xh);CHKERRQ(ierr);
185 if (i>= start && i<end) {
186 xv = 2.*dx;
187 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr);
188 }
189 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr);
190 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr);
191 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr);
192
193 /* compute the 3rd value */
194 ierr = VecCopy(X,Xh);CHKERRQ(ierr);
195 if (i>= start && i<end) {
196 xv = -dx;
197 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr);
198 }
199 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr);
200 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr);
201 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr);
202
203 if (i >= start && i<end) {
204 /* set this entry to be the gradient of the objective */
205 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
206 if (PetscAbsScalar(fv) > eps) {
207 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr);
208 } else {
209 fv = 0.;
210 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr);
211 }
212 }
213 }
214
215 ierr = VecDestroy(&Xh);CHKERRQ(ierr);
216
217 ierr = VecAssemblyBegin(F);CHKERRQ(ierr);
218 ierr = VecAssemblyEnd(F);CHKERRQ(ierr);
219 PetscFunctionReturn(0);
220 }
221