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