1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Authors: f. Gillet-Chaulet (IGE, Grenoble,France) 26! * Email: 27! * Web: http://elmerice.elmerfem.org 28! * 29! * Original Date: April 2020 30! * 31! ***************************************************************************** 32SUBROUTINE Adjoint_GradientValidation_init0(Model,Solver,dt,TransientSimulation ) 33!------------------------------------------------------------------------------ 34 USE DefUtils 35 IMPLICIT NONE 36!------------------------------------------------------------------------------ 37 TYPE(Solver_t), TARGET :: Solver 38 TYPE(Model_t) :: Model 39 REAL(KIND=dp) :: dt 40 LOGICAL :: TransientSimulation 41!------------------------------------------------------------------------------ 42!------------------------------------------------------------------------------ 43! Local variables 44!------------------------------------------------------------------------------ 45 CHARACTER(LEN=MAX_NAME_LEN) :: Name 46 47 Name = ListGetString( Solver % Values, 'Equation',UnFoundFatal=.TRUE.) 48 CALL ListAddNewString( Solver % Values,'Variable',& 49 '-nooutput '//TRIM(Name)//'_var') 50 CALL ListAddLogical(Solver % Values, 'Optimize Bandwidth',.FALSE.) 51 CALL ListAddInteger(Solver % Values, 'Nonlinear System Norm Degree',0) 52END SUBROUTINE Adjoint_GradientValidation_init0 53! ***************************************************************************** 54SUBROUTINE Adjoint_GradientValidation ( Model,Solver,dt,TransientSimulation ) 55! ***************************************************************************** 56! Compare the total derivative of the cost function computed as: 57! (1) dJ=P.G where P is a perturbation vector of the variable of interest 58! G is the gradient of the cost function computed by an inverse method 59! (2) [J(V+hP)-J(V)]/h : forward finite difference computation of the derivative 60! V is the variable of interest 61! h is the step size 62! 63! 64! Compute (1) from at the first iteration and update V=Vini+hP, h=1 65! Compute (2) for all the other iteration with h^i+1=h^i/2 66! 67! Serial/parallel 2D/3D 68! 69! Keyword in Solver section of the .sif: 70! Cost Variable Name 71! Optimized Variable Name 72! Perturbed Variable Name !optional: take -g if not given 73! Gradient Variable Name 74! Result File 75! 76! Output: in result File: h , abs((1)-(2))/(1) , (1), (2) 77! 78! 79!------------------------------------------------------------------------------ 80!****************************************************************************** 81 USE DefUtils 82 IMPLICIT NONE 83!------------------------------------------------------------------------------ 84 TYPE(Solver_t) :: Solver 85 TYPE(Model_t) :: Model 86 87 REAL(KIND=dp) :: dt 88 LOGICAL :: TransientSimulation 89! 90 CHARACTER(LEN=MAX_NAME_LEN) :: SolverName 91 TYPE(Element_t),POINTER :: Element 92 TYPE(ValueList_t), POINTER :: SolverParams 93 TYPE(Variable_t), POINTER :: Var,PVar,CostVar,GradVar 94 REAL(KIND=dp), POINTER :: Values(:),PValues(:),CostValues(:),GradValues(:) 95 INTEGER, POINTER :: Perm(:),PPerm(:),GradPerm(:) 96 INTEGER, POINTER :: NodeIndexes(:) 97 98 REAL(KIND=dp),allocatable :: x(:),xp(:),g(:) 99 REAL(KIND=dp) :: J,J0,h,dJ,dJd 100 101 integer :: i,c,t,n,NMAX,NActiveNodes 102 integer :: ierr 103 integer,allocatable :: ActiveNodes(:) 104 integer,allocatable :: NewNode(:) 105 integer,parameter :: io=20 106 integer :: MyPe=-1 107 integer,SAVE :: DOFs 108 109 Logical :: FirstVisit=.true.,Found,UnFoundFatal=.TRUE. 110 Logical :: Parallel 111 Logical :: haveP 112 logical,allocatable :: VisitedNode(:) 113 114 CHARACTER(LEN=MAX_NAME_LEN) :: CostSolName,VarSolName,GradSolName,PSolName,ResultFile 115 CHARACTER*10 :: date,temps 116 117! 118 save FirstVisit 119 save SolverName 120 save CostSolName,VarSolName,GradSolName,ResultFile 121 save ActiveNodes,NActiveNodes 122 save x,xp 123 save J0,h,dJ 124 save MyPe 125 save Parallel 126 127 128! Read Constant from sif solver section 129 IF(FirstVisit) Then 130 131 WRITE(SolverName, '(A)') 'GradientValidation' 132 133 ! Check we have a parallel run 134 Parallel = .FALSE. 135 IF(ASSOCIATED(Solver % Matrix % ParMatrix)) Then 136 IF ( Solver % Matrix % ParMatrix % ParEnv % PEs > 1 ) Parallel =.True. 137 MyPe=ParEnv % MyPe 138 End if 139 140 141 SolverParams => GetSolverParams() 142 143 CostSolName = ListGetString( SolverParams,'Cost Variable Name', UnFoundFatal=.TRUE.) 144 CostVar => VariableGet( Solver % Mesh % Variables, CostSolName,UnFoundFatal=UnFoundFatal) 145 CostValues => CostVar % Values 146 147 VarSolName = ListGetString( SolverParams,'Optimized Variable Name', UnFoundFatal=.TRUE.) 148 Var => VariableGet( Solver % Mesh % Variables, VarSolName,UnFoundFatal=UnFoundFatal) 149 Values => Var % Values 150 Perm => Var % Perm 151 DOFs = Var % DOFs 152 153 GradSolName = ListGetString( SolverParams,'Gradient Variable Name', UnFoundFatal=.TRUE.) 154 GradVar => VariableGet( Solver % Mesh % Variables, GradSolName,UnFoundFatal=UnFoundFatal) 155 GradValues => GradVar % Values 156 GradPerm => GradVar % Perm 157 IF (GradVar%DOFs.NE.DOFs) & 158 CALL FATAL(SolverName,'DOFs not corresponding for Gradient Variable') 159 160 PSolName = ListGetString( SolverParams,'Perturbation Variable Name', Found=HaveP) 161 IF (HAVEP) THEN 162 PVar => VariableGet( Solver % Mesh % Variables, PSolName,UnFoundFatal=UnFoundFatal) 163 PValues => PVar % Values 164 PPerm => PVar % Perm 165 IF (PVar%DOFs.NE.DOFs) & 166 CALL FATAL(SolverName,'DOFs not corresponding for Perturbation Variable') 167 ENDIF 168 169!!!!!!!!!!!!find active nodes 170 NMAX=Solver % Mesh % NumberOfNodes 171 allocate(VisitedNode(NMAX),NewNode(NMAX)) 172 VisitedNode=.false. 173 NewNode=-1 174 175 NActiveNodes=0 176 DO t=1,Solver % NumberOfActiveElements 177 Element => GetActiveElement(t) 178 n = GetElementNOFNodes() 179 NodeIndexes => Element % NodeIndexes 180 Do i=1,n 181 if (VisitedNode(NodeIndexes(i))) then 182 cycle 183 else 184 VisitedNode(NodeIndexes(i))=.true. 185 NActiveNodes=NActiveNodes+1 186 NewNode(NActiveNodes)=NodeIndexes(i) 187 endif 188 End do 189 End do 190 191 if (NActiveNodes.eq.0) THEN 192 WRITE(Message,'(A)') 'NActiveNodes = 0 !!!' 193 CALL FATAL(SolverName,Message) 194 End if 195 196 allocate(ActiveNodes(NActiveNodes),x(DOFS*NActiveNodes),xp(DOFs*NActiveNodes),g(DOFs*NActiveNodes)) 197 ActiveNodes(1:NActiveNodes)=NewNode(1:NActiveNodes) 198 199 deallocate(VisitedNode,NewNode) 200 201!!!!!!! Solver Params 202 203 204 ResultFile=GetString( SolverParams,'Result File',Found) 205 IF(Found) Then 206 IF ((Parallel.AND.(MyPe.EQ.0)).OR.(.NOT.Parallel)) Then 207 open(io,file=trim(ResultFile)) 208 CALL DATE_AND_TIME(date,temps) 209 write(io,'(a1,a2,a1,a2,a1,a4,5x,a2,a1,a2,a1,a2)')'#',date(5:6),'/',date(7:8),'/',date(1:4), & 210 temps(1:2),':',temps(3:4),':',temps(5:6) 211 write(io,'(A)') '# step size, relative error, Adjoint total der., FD total der.' 212 close(io) 213 ENDIF 214 ELSE 215 CALL FATAL(SolverName,'Keyword <Result File> Not Found') 216 ENDIF 217 218!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 219 DO i=1,NActiveNodes 220 DO c=1,DOFs 221 x(DOFs*(i-1)+c)=Values(DOFs*(Perm(ActiveNodes(i))-1)+c) 222 g(DOFs*(i-1)+c)=GradValues(DOFs *(GradPerm(ActiveNodes(i))-1)+c) 223 IF (HaveP) THEN 224 xp(DOFs*(i-1)+c)=PValues(DOFs*(PPerm(ActiveNodes(i))-1)+c) 225 ELSE 226 xp(DOFs*(i-1)+c)=-g(DOFs*(i-1)+c) 227 ENDIF 228 END DO 229 END DO 230 231 !!!!!! total derivative from gradient 232 dJ=0._dp 233 dJ=SUM(xp(:)*g(:)) 234 deallocate(g) 235 236 IF (Parallel) THEN 237 CALL MPI_ALLREDUCE(dJ,dJ,1,& 238 MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr) 239 ENDIF 240 241 !!!!! Store cost value at first iteration 242 J0=CostValues(1) 243 244!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 245 !!! new values for x=x+h*xp 246 h=1.0_dp 247 Do i=1,NActiveNodes 248 Do c=1,DOFs 249 Values(DOFs*(Perm(ActiveNodes(i))-1)+c)=x(DOFs*(i-1)+c)+h*xp(DOFs*(i-1)+c) 250 End do 251 END DO 252 253 254 FirstVisit=.FALSE. 255 Return 256 End if 257!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 258 259 CostVar => VariableGet( Solver % Mesh % Variables, CostSolName,UnFoundFatal=UnFoundFatal) 260 CostValues => CostVar % Values 261 262 Var => VariableGet( Solver % Mesh % Variables, VarSolName,UnFoundFatal=UnFoundFatal) 263 Values => Var % Values 264 Perm => Var % Perm 265 266 267 J=CostValues(1) 268 269 dJd=(J-J0)/h 270 271 IF (Parallel) MyPe=ParEnv % MyPe 272 IF ((Parallel.AND.(MyPe.EQ.0)).OR.(.NOT.Parallel)) Then 273 open(io,file=trim(ResultFile),position='append') 274 write(io,'(4(e15.8,2x))') h,abs(dJ-dJd)/abs(dJ),dJ,dJd 275 close(io) 276 ENDIF 277 278 h=h/2.0_dp 279 Do i=1,NActiveNodes 280 Do c=1,DOFs 281 Values(DOFs*(Perm(ActiveNodes(i))-1)+c)=x(DOFs*(i-1)+c)+h*xp(DOFs*(i-1)+c) 282 End do 283 END DO 284 285 Return 286!------------------------------------------------------------------------------ 287END SUBROUTINE Adjoint_GradientValidation 288!------------------------------------------------------------------------------ 289 290 291