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 . F - current function/gradient 14 . obj - real to hold the objective value 15 - ctx - optional user-defined objective context 16 17 .seealso: SNESSetFunction(), SNESGetFunction(), SNESSetObjectiveFunction(), SNESGetObjectiveFunction() 18 M*/ 19 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "SNESSetObjective" 23 /*@C 24 SNESSetObjective - Sets the objective function minimized by the SNES methods. 25 26 Logically Collective on SNES 27 28 Input Parameters: 29 + snes - the SNES context 30 . SNESObjectiveFunction - objective evaluation routine 31 - ctx - [optional] user-defined context for private data for the 32 function evaluation routine (may be PETSC_NULL) 33 34 Level: beginner 35 36 Note: If not provided then this defaults to the two norm of the function evaluation (set with SNESSetFunction()) 37 38 .keywords: SNES, nonlinear, set, objective 39 40 .seealso: SNESGetObjective(), SNESComputeObjective(), SNESSetFunction(), SNESSetJacobian(), SNESObjectiveFunction 41 @*/ 42 PetscErrorCode SNESSetObjective(SNES snes,PetscErrorCode (*SNESObjectiveFunction)(SNES,Vec,PetscReal*,void*),void *ctx) 43 { 44 PetscErrorCode ierr; 45 DM dm; 46 47 PetscFunctionBegin; 48 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 49 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 50 ierr = DMSNESSetObjective(dm,SNESObjectiveFunction,ctx);CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESGetObjective" 56 /*@C 57 SNESGetObjective - Returns the objective function. 58 59 Not Collective 60 61 Input Parameter: 62 . snes - the SNES context 63 64 Output Parameter: 65 + func - the function (or PETSC_NULL) 66 - ctx - the function context (or PETSC_NULL) 67 68 Level: advanced 69 70 .keywords: SNES, nonlinear, get, objective 71 72 .seealso: SNESSetObjective(), SNESGetSolution() 73 @*/ 74 PetscErrorCode SNESGetObjective(SNES snes,PetscErrorCode (**SNESObjectiveFunction)(SNES,Vec,PetscReal*,void*),void **ctx) 75 { 76 PetscErrorCode ierr; 77 DM dm; 78 79 PetscFunctionBegin; 80 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 81 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 82 ierr = DMSNESGetObjective(dm,SNESObjectiveFunction,ctx);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "SNESComputeObjective" 88 /*@C 89 SNESComputeObjective - Computes the objective. 90 91 Collective on SNES 92 93 Input Parameter: 94 + snes - the SNES context 95 - X - the state vector 96 97 Output Parameter: 98 . ob - the objective value 99 100 Level: advanced 101 102 .keywords: SNES, nonlinear, compute, objective 103 104 .seealso: SNESSetObjective(), SNESGetSolution() 105 @*/ 106 PetscErrorCode SNESComputeObjective(SNES snes,Vec X,PetscReal *ob) 107 { 108 PetscErrorCode ierr; 109 DM dm; 110 DMSNES sdm; 111 112 PetscFunctionBegin; 113 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 114 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 115 PetscValidPointer(ob,3); 116 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 117 ierr = DMGetDMSNES(dm,&sdm);CHKERRQ(ierr); 118 if (sdm->ops->computeobjective) { 119 ierr = (sdm->ops->computeobjective)(snes,X,ob,sdm->objectivectx);CHKERRQ(ierr); 120 } else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetObjective() before SNESComputeObjective()."); 121 PetscFunctionReturn(0); 122 } 123 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "SNESDefaultObjectiveComputeFunctionFD" 127 PetscErrorCode SNESDefaultObjectiveComputeFunctionFD(SNES snes,Vec X,Vec F,void *ctx) 128 { 129 /* Quadratically interpolates the change in the objective based upon a change in a particular direction */ 130 131 Vec Xh; 132 PetscErrorCode ierr; 133 PetscInt i,N,start,end; 134 PetscReal ob,ob1,ob2,ob3,fob,dx,eps=1e-6; 135 PetscScalar fv,xv; 136 137 PetscFunctionBegin; 138 ierr = VecDuplicate(X,&Xh);CHKERRQ(ierr); 139 ierr = PetscOptionsReal("-snes_fd_function_eps","Tolerance for nonzero entries in fd function","None",eps,&eps,PETSC_NULL);CHKERRQ(ierr); 140 ierr = VecSet(F,0.);CHKERRQ(ierr); 141 142 ierr = VecNorm(X,NORM_2,&fob);CHKERRQ(ierr); 143 144 ierr = VecGetSize(X,&N);CHKERRQ(ierr); 145 ierr = VecGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 146 ierr = SNESComputeObjective(snes,X,&ob);CHKERRQ(ierr); 147 148 if (fob > 0.) dx =1e-6*fob; 149 else dx = 1e-6; 150 151 for (i=0; i<N; i++) { 152 /* compute the 1st value */ 153 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 154 if (i>= start && i<end) { 155 xv = dx; 156 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 157 } 158 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 159 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 160 ierr = SNESComputeObjective(snes,Xh,&ob1);CHKERRQ(ierr); 161 162 /* compute the 2nd value */ 163 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 164 if (i>= start && i<end) { 165 xv = 2.*dx; 166 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 167 } 168 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 169 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 170 ierr = SNESComputeObjective(snes,Xh,&ob2);CHKERRQ(ierr); 171 172 /* compute the 3rd value */ 173 ierr = VecCopy(X,Xh);CHKERRQ(ierr); 174 if (i>= start && i<end) { 175 xv = -dx; 176 ierr = VecSetValues(Xh,1,&i,&xv,ADD_VALUES);CHKERRQ(ierr); 177 } 178 ierr = VecAssemblyBegin(Xh);CHKERRQ(ierr); 179 ierr = VecAssemblyEnd(Xh);CHKERRQ(ierr); 180 ierr = SNESComputeObjective(snes,Xh,&ob3);CHKERRQ(ierr); 181 182 if (i >= start && i<end) { 183 /* set this entry to be the gradient of the objective */ 184 fv = (-ob2 + 6.*ob1 - 3.*ob -2.*ob3) / (6.*dx); 185 if (PetscAbsScalar(fv) > eps) { 186 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 187 } else { 188 fv = 0.; 189 ierr = VecSetValues(F,1,&i,&fv,INSERT_VALUES);CHKERRQ(ierr); 190 } 191 } 192 } 193 194 ierr = VecDestroy(&Xh);CHKERRQ(ierr); 195 196 ierr = VecAssemblyBegin(F);CHKERRQ(ierr); 197 ierr = VecAssemblyEnd(F);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200