Actual source code: snesob.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <petsc/private/snesimpl.h>

  3: /*MC
  4:     SNESObjectiveFunction - functional form used to convey the objective function to the nonlinear solver

  6:      Synopsis:
  7:  #include <petscsnes.h>
  8:        SNESObjectiveFunction(SNES snes,Vec x,PetscReal *obj,void *ctx);

 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

 16:    Level: advanced

 18: .seealso:   SNESSetFunction(), SNESGetFunction(), SNESSetObjective(), SNESGetObjective(), SNESJacobianFunction, SNESFunction
 19: M*/


 22: /*@C
 23:    SNESSetObjective - Sets the objective function minimized by some of the SNES linesearch methods.

 25:    Logically Collective on SNES

 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)

 33:    Level: intermediate

 35:    Note: This is not used in the SNESLINESEARCHCP line search.

 37:          If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction())

 39: .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction
 40: @*/
 41: PetscErrorCode  SNESSetObjective(SNES snes,PetscErrorCode (*obj)(SNES,Vec,PetscReal*,void*),void *ctx)
 42: {
 44:   DM             dm;

 48:   SNESGetDM(snes,&dm);
 49:   DMSNESSetObjective(dm,obj,ctx);
 50:   return(0);
 51: }

 53: /*@C
 54:    SNESGetObjective - Returns the objective function.

 56:    Not Collective

 58:    Input Parameter:
 59: .  snes - the SNES context

 61:    Output Parameter:
 62: +  obj - objective evaluation routine (or NULL); see SNESObjectFunction for details
 63: -  ctx - the function context (or NULL)

 65:    Level: advanced

 67: .seealso: SNESSetObjective(), SNESGetSolution()
 68: @*/
 69: PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**obj)(SNES,Vec,PetscReal*,void*),void **ctx)
 70: {
 72:   DM             dm;

 76:   SNESGetDM(snes,&dm);
 77:   DMSNESGetObjective(dm,obj,ctx);
 78:   return(0);
 79: }

 81: /*@C
 82:    SNESComputeObjective - Computes the objective.

 84:    Collective on SNES

 86:    Input Parameter:
 87: +  snes - the SNES context
 88: -  X    - the state vector

 90:    Output Parameter:
 91: .  ob   - the objective value

 93:    Level: advanced

 95: .seealso: SNESSetObjective(), SNESGetSolution()
 96: @*/
 97: PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob)
 98: {
100:   DM             dm;
101:   DMSNES         sdm;

107:   SNESGetDM(snes,&dm);
108:   DMGetDMSNES(dm,&sdm);
109:   if (sdm->ops->computeobjective) {
110:     PetscLogEventBegin(SNES_ObjectiveEval,snes,X,0,0);
111:     (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);
112:     PetscLogEventEnd(SNES_ObjectiveEval,snes,X,0,0);
113:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective().");
114:   return(0);
115: }


118: /*@C
119:    SNESObjectiveComputeFunctionDefaultFD - Computes the gradient of a user provided objective

121:    Collective on SNES

123:    Input Parameter:
124: +  snes - the SNES context
125: .  X    - the state vector
126: -  ctx  - the (ignored) function context

128:    Output Parameter:
129: .  F   - the function value

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

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.

142:    Note that this uses quadratic interpolation of the objective to form each value in the function.

144:    Level: advanced

146: .seealso: SNESSetFunction(), SNESComputeObjective(), SNESComputeJacobianDefault()
147: @*/
148: PetscErrorCode SNESObjectiveComputeFunctionDefaultFD(SNES snes,Vec X,Vec F,void *ctx)
149: {
150:   Vec            Xh;
152:   PetscInt       i,N,start,end;
153:   PetscReal      ob,ob1,ob2,ob3,fob,dx,eps=1e-6;
154:   PetscScalar    fv,xv;

157:   VecDuplicate(X,&Xh);
158:   PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing parameters","SNES");
159:   PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,NULL);
160:   PetscOptionsEnd();
161:   VecSet(F,0.);

163:   VecNorm(X,NORM_2,&fob);

165:   VecGetSize(X,&N);
166:   VecGetOwnershipRange(X,&start,&end);
167:   SNESComputeObjective(snes,X,&ob);

169:   if (fob > 0.) dx =1e-6*fob;
170:   else dx = 1e-6;

172:   for (i=0; i<N; i++) {
173:     /* compute the 1st value */
174:     VecCopy(X,Xh);
175:     if (i>= start && i<end) {
176:       xv   = dx;
177:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
178:     }
179:     VecAssemblyBegin(Xh);
180:     VecAssemblyEnd(Xh);
181:     SNESComputeObjective(snes,Xh,&ob1);

183:     /* compute the 2nd value */
184:     VecCopy(X,Xh);
185:     if (i>= start && i<end) {
186:       xv   = 2.*dx;
187:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
188:     }
189:     VecAssemblyBegin(Xh);
190:     VecAssemblyEnd(Xh);
191:     SNESComputeObjective(snes,Xh,&ob2);

193:     /* compute the 3rd value */
194:     VecCopy(X,Xh);
195:     if (i>= start && i<end) {
196:       xv   = -dx;
197:       VecSetValues(Xh,1,&i,&xv,ADD_VALUES);
198:     }
199:     VecAssemblyBegin(Xh);
200:     VecAssemblyEnd(Xh);
201:     SNESComputeObjective(snes,Xh,&ob3);

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:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
208:       } else {
209:         fv   = 0.;
210:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
211:       }
212:     }
213:   }

215:   VecDestroy(&Xh);

217:   VecAssemblyBegin(F);
218:   VecAssemblyEnd(F);
219:   return(0);
220: }