Actual source code: snesob.c

petsc-master 2019-10-16
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: .      F - current function/gradient
 14: .      obj - real to hold the objective value
 15: -      ctx - optional user-defined objective context

 17:    Level: advanced

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


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

 26:    Logically Collective on SNES

 28:    Input Parameters:
 29: +  snes - the SNES context
 30: .  obj - objective evaluation routine; see SNESObjectiveFunction for details
 31: -  ctx - [optional] user-defined context for private data for the
 32:          function evaluation routine (may be NULL)

 34:    Level: intermediate

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

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

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

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

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

 57:    Not Collective

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

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

 66:    Level: advanced

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

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

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

 85:    Collective on SNES

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

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

 94:    Level: advanced

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

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


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

122:    Collective on SNES

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

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

132:    Options Database Key:
133: +  -snes_fd_function_eps - The differencing parameter
134: -  -snes_fd_function - Compute function from user provided objective with finite difference

136:    Notes:
137:    SNESObjectiveComputeFunctionDefaultFD is similar in character to SNESComputeJacobianDefault.
138:    Therefore, it should be used for debugging purposes only.  Using it in conjunction with
139:    SNESComputeJacobianDefault is excessively costly and produces a Jacobian that is quite
140:    noisy.  This is often necessary, but should be done with a grain of salt, even when debugging
141:    small problems.

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

145:    Level: advanced

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

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

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

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

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

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

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

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

204:     if (i >= start && i<end) {
205:       /* set this entry to be the gradient of the objective */
206:       fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx);
207:       if (PetscAbsScalar(fv) > eps) {
208:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
209:       } else {
210:         fv   = 0.;
211:         VecSetValues(F,1,&i,&fv,INSERT_VALUES);
212:       }
213:     }
214:   }

216:   VecDestroy(&Xh);

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