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:
26! *  Email:
27! *  Web:     http://elmerice.elmerfem.org
28! *
29! *  Original Date:
30! *
31! *****************************************************************************
32! Compute the integrated gradient of the Cost function for the
33! Adjoint Inverse Problem
34!
35! Serial/Parallel(not halo?)   and 2D/3D
36!
37! Need:
38! - Navier-stokes and Adjoint Problems Solutions
39! - Name of Beta variable and Grad Variable
40! - Power formulation if 'Slip Coef'=10^Beta and optimization on Beta
41! - Beta2 formulation if 'Slip Coef'=Beta^2  and optimization on Beta
42! - Lambda : the regularization coefficient (Jr=0.5 Lambda (dBeta/dx)^2)
43!
44! If DJDBeta should be set to zero for floating ice shelves the following is
45! to be used (default is not to do this):
46! - FreeSlipShelves (logical, default false)
47! - mask name (string, default GroundedMask)
48! Note that if FreeSlipShelves is true not only is DJDBeta set to zero for
49! floating ice, but also the regularisation term NodalRegb.
50!
51! *****************************************************************************
52SUBROUTINE DJDBeta_Adjoint( Model,Solver,dt,TransientSimulation )
53!------------------------------------------------------------------------------
54!******************************************************************************
55  USE DefUtils
56  IMPLICIT NONE
57!------------------------------------------------------------------------------
58  TYPE(Solver_t) :: Solver
59  TYPE(Model_t) :: Model
60
61  REAL(KIND=dp) :: dt
62  LOGICAL :: TransientSimulation
63!
64  TYPE(ValueList_t), POINTER :: BC
65
66  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName,NeumannSolName,AdjointSolName
67  CHARACTER(LEN=MAX_NAME_LEN) :: VarSolName,GradSolName
68  TYPE(Element_t),POINTER ::  Element
69  TYPE(Variable_t), POINTER :: PointerToVariable, BetaVariable, VeloSolN,VeloSolD
70  TYPE(ValueList_t), POINTER :: SolverParams
71  TYPE(Nodes_t) :: ElementNodes
72  TYPE(GaussIntegrationPoints_t) :: IntegStuff
73
74  REAL(KIND=dp), POINTER ::  VariableValues(:),VelocityN(:),VelocityD(:),BetaValues(:)
75  INTEGER, POINTER :: Permutation(:), VeloNPerm(:),VeloDPerm(:),BetaPerm(:),NodeIndexes(:)
76
77  real(kind=dp),allocatable :: VisitedNode(:),db(:),Basis(:),dBasisdx(:,:)
78  real(kind=dp),allocatable :: nodalbetab(:),NodalRegb(:)
79  real(kind=dp) :: betab
80  real(kind=dp) :: u,v,w,SqrtElementMetric,s
81  real(kind=dp) :: Lambda
82  REAL(KIND=dp) :: Normal(3),Tangent(3),Tangent2(3),Vect(3)
83
84  integer :: i,j,k,e,t,n,NMAX,NActiveNodes,DIM
85  integer :: p,q
86
87  logical :: PowerFormulation,Beta2Formulation
88  Logical ::  Firsttime=.true.,Found,stat,UnFoundFatal=.TRUE.
89  Logical :: NormalTangential1,NormalTangential2
90
91  ! Variables for setting DJDBeta to zero for ice shelves.
92  TYPE(Variable_t), POINTER   :: PointerToMask => NULL()
93  REAL(KIND=dp), POINTER      :: MaskValues(:) => NULL()
94  INTEGER, POINTER            :: MaskPerm(:) => NULL()
95  LOGICAL                     :: FreeSlipShelves
96  CHARACTER(LEN=MAX_NAME_LEN) :: MaskName
97
98  save FreeSlipShelves,MaskName
99  save SolverName,NeumannSolName,AdjointSolName,VarSolName,GradSolName
100  save VisitedNode,db,Basis,dBasisdx,nodalbetab,NodalRegb
101  save Firsttime,DIM,Lambda
102  save ElementNodes
103  save PowerFormulation,Beta2Formulation
104
105  If (Firsttime) then
106
107     DIM = CoordinateSystemDimension()
108     WRITE(SolverName, '(A)') 'DJDBeta_Adjoint'
109
110   CALL Info(SolverName,'***********************',level=0)
111   CALL Info(SolverName,' This solver has been replaced by:',level=0)
112   CALL Info(SolverName,'   AdjointStokes_GradientBetaSolver  ',level=0)
113   CALL Info(SolverName,' See documentation under:   ',level=0)
114   CALL Info(SolverName,'   elmerice/Solvers/Documentation   ',level=0)
115   CALL Info(SolverName,'***********************',level=0)
116   CALL FATAL(SolverName,' Use new solver !!')
117
118     NMAX=Solver % Mesh % NumberOfNodes
119     allocate(VisitedNode(NMAX),db(NMAX),  &
120          nodalbetab(Model %  MaxElementNodes),&
121          NodalRegb(Model %  MaxElementNodes),&
122          Basis(Model % MaxElementNodes),  &
123          dBasisdx(Model % MaxElementNodes,3))
124
125!!!!!!!!!!! get Solver Variables
126     SolverParams => GetSolverParams()
127
128     NeumannSolName =  GetString( SolverParams,'Flow Solution Name', Found)
129     IF(.NOT.Found) THEN
130        CALL WARN(SolverName,'Keyword >Neumann Solution Name< not found in section >Solver<')
131        CALL WARN(SolverName,'Taking default value >Flow Solution<')
132        WRITE(NeumannSolName,'(A)') 'Flow Solution'
133     END IF
134     AdjointSolName =  GetString( SolverParams,'Adjoint Solution Name', Found)
135     IF(.NOT.Found) THEN
136        CALL WARN(SolverName,'Keyword >Adjoint Solution Name< not found in section >Solver<')
137        CALL WARN(SolverName,'Taking default value >Adjoint<')
138        WRITE(AdjointSolName,'(A)') 'Adjoint'
139     END IF
140     VarSolName =  GetString( SolverParams,'Optimized Variable Name', Found)
141     IF(.NOT.Found) THEN
142        CALL WARN(SolverName,'Keyword >Optimized Variable Name< not found  in section >Solver<')
143        CALL WARN(SolverName,'Taking default value >Beta<')
144        WRITE(VarSolName,'(A)') 'Beta'
145     END IF
146     GradSolName =  GetString( SolverParams,'Gradient Variable Name', Found)
147     IF(.NOT.Found) THEN
148        CALL WARN(SolverName,'Keyword >Gradient Variable Name< not found  in section >Solver<')
149        CALL WARN(SolverName,'Taking default value >DJDB<')
150        WRITE(GradSolName,'(A)') 'DJDB'
151     END IF
152     PowerFormulation=GetLogical( SolverParams, 'PowerFormulation', Found)
153     IF(.NOT.Found) THEN
154        CALL WARN(SolverName,'Keyword >PowerFormulation< not found  in section >Equation<')
155        CALL WARN(SolverName,'Taking default value >FALSE<')
156        PowerFormulation=.FALSE.
157     END IF
158     Beta2Formulation=GetLogical( SolverParams, 'Beta2Formulation', Found)
159     IF(.NOT.Found) THEN
160        CALL WARN(SolverName,'Keyword >Beta2Formulation< not found  in section >Equation<')
161        CALL WARN(SolverName,'Taking default value >FALSE<')
162        Beta2Formulation=.FALSE.
163     END IF
164     IF (PowerFormulation.and.Beta2Formulation) then
165        WRITE(Message,'(A)') 'Can t be PowerFormulation and Beta2Formulation in the same time'
166        CALL FATAL(SolverName,Message)
167     End if
168     Lambda =  GetConstReal( SolverParams,'Lambda', Found)
169     IF(.NOT.Found) THEN
170        CALL WARN(SolverName,'Keyword >Lambda< not found  in section  >Solver<')
171        CALL WARN(SolverName,'Taking default value Lambda=0.0')
172        Lambda = 0.0
173     End if
174
175     FreeSlipShelves=GetLogical( SolverParams, 'FreeSlipShelves', Found)
176     IF(.NOT.Found) THEN
177        CALL WARN(SolverName,'Keyword >FreeSlipShelves< not found in solver params')
178        CALL WARN(SolverName,'Taking default value >FALSE<')
179        FreeSlipShelves=.FALSE.
180     END IF
181     IF (FreeSlipShelves) THEN
182        MaskName =  GetString( SolverParams,'mask name', Found)
183        IF(.NOT.Found) THEN
184           CALL WARN(SolverName,'Keyword >mask name< not found in solver section')
185           CALL WARN(SolverName,'Taking default value >GroundedMask<')
186           WRITE(MaskName,'(A)') 'GroundedMask'
187        END IF
188     END IF
189
190!!! End of First visit
191     Firsttime=.false.
192  Endif
193
194  PointerToVariable => VariableGet( Model % Mesh % Variables, GradSolName,UnFoundFatal=UnFoundFatal)
195  VariableValues => PointerToVariable % Values
196  Permutation => PointerToVariable % Perm
197  VariableValues=0._dp
198
199  BetaVariable => VariableGet( Model % Mesh % Variables, VarSolName,UnFoundFatal=UnFoundFatal)
200  BetaValues => BetaVariable % Values
201  BetaPerm => BetaVariable % Perm
202
203  VeloSolN => VariableGet( Model % Mesh % Variables, NeumannSolName,UnFoundFatal=UnFoundFatal)
204  VelocityN => VeloSolN % Values
205  VeloNPerm => VeloSolN % Perm
206
207  VeloSolD => VariableGet( Model % Mesh % Variables, AdjointSolName,UnFoundFatal=UnFoundFatal)
208  VelocityD => VeloSolD % Values
209  VeloDPerm => VeloSolD % Perm
210
211  IF (FreeSlipShelves) THEN
212     PointerToMask => VariableGet( Model % Variables, MaskName, UnFoundFatal=.TRUE.)
213     MaskValues => PointerToMask % Values
214     MaskPerm => PointerToMask % Perm
215  END IF
216
217  VisitedNode=0.0_dp
218  db=0.0_dp
219
220  DO e=1,Solver % NumberOfActiveElements
221     Element => GetActiveElement(e)
222     CALL GetElementNodes( ElementNodes )
223     n = GetElementNOFNodes()
224     NodeIndexes => Element % NodeIndexes
225
226     ! Compute Nodal Value of DJDBeta
227     BC => GetBC()
228     if (.NOT.ASSOCIATED(BC)) CYCLE
229
230     NormalTangential1 = GetLogical( BC, &
231          'Normal-Tangential Velocity', Found )
232     IF (.NOT.Found) then
233        NormalTangential1 = GetLogical( BC, &
234             'Normal-Tangential '//trim(NeumannSolName), Found)
235     END IF
236     NormalTangential2 = GetLogical( BC, &
237          'Normal-Tangential '//trim(AdjointSolName), Found)
238     IF (NormalTangential1.NEQV.NormalTangential2) then
239        WRITE(Message,'(A,I1,A,I1)') &
240             'NormalTangential Velocity is : ',NormalTangential1, &
241             'But NormalTangential Adjoint is : ',NormalTangential2
242        CALL FATAL(SolverName,Message)
243     ENDIF
244     IF (.NOT.NormalTangential1) then
245        WRITE(Message,'(A)') &
246             'ALWAYS USE Normal-Tangential COORDINATES with SlipCoef 2=SlipCoef 3'
247        CALL FATAL(SolverName,Message)
248     ENDIF
249
250     VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp
251
252     ! Compute Integrated Nodal Value of DJDBeta
253     nodalbetab=0.0_dp
254     NodalRegb=0.0_dp
255
256
257     IntegStuff = GaussPoints( Element )
258     DO t=1,IntegStuff % n
259        U = IntegStuff % u(t)
260        V = IntegStuff % v(t)
261        W = IntegStuff % w(t)
262        stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, &
263             Basis,dBasisdx )
264
265        s = SqrtElementMetric * IntegStuff % s(t)
266
267        ! compute gradient from Stokes and adjoint computation
268        ! follow the computation of the stiffMatrix as done in the NS solver
269        Normal = NormalVector( Element, ElementNodes, u,v,.TRUE. )
270        SELECT CASE( Element % TYPE % DIMENSION )
271        CASE(1)
272           Tangent(1) =  Normal(2)
273           Tangent(2) = -Normal(1)
274           Tangent(3) =  0.0_dp
275           Tangent2   =  0.0_dp
276        CASE(2)
277           CALL TangentDirections( Normal, Tangent, Tangent2 )
278        END SELECT
279
280
281        betab=0.0_dp
282        Do p=1,n
283           Do q=1,n
284
285              Do i=2,dim
286                 SELECT CASE(i)
287                 CASE(2)
288                    Vect = Tangent
289                 CASE(3)
290                    Vect = Tangent2
291                 END SELECT
292
293                 Do j=1,DIM
294                    Do k=1,DIM
295                       betab = betab + s *  Basis(q) * Basis(p) * Vect(j) * Vect(k) * &
296                            (- VelocityN((DIM+1)*(VeloNPerm(NodeIndexes(q))-1)+k) * &
297                            VelocityD((DIM+1)*(VeloDPerm(NodeIndexes(p))-1)+j))
298                    End Do !on k
299                 End Do !on j
300              End Do !on i
301           End Do !on q
302        End Do !on p
303
304        nodalbetab(1:n)=nodalbetab(1:n)+betab*Basis(1:n)
305
306        If (Lambda /= 0.0) then
307           NodalRegb(1:n)=NodalRegb(1:n)+&
308                s*Lambda*(SUM(dBasisdx(1:n,1)*BetaValues(BetaPerm(NodeIndexes(1:n))))*dBasisdx(1:n,1))
309           IF (DIM.eq.3) then
310              NodalRegb(1:n)=NodalRegb(1:n)+&
311                   s*Lambda*(SUM(dBasisdx(1:n,2)*BetaValues(BetaPerm(NodeIndexes(1:n))))*dBasisdx(1:n,2))
312           End if
313        End if
314     End DO ! on IPs
315
316     IF (PowerFormulation) then
317        nodalbetab(1:n)=nodalbetab(1:n)*(10**(BetaValues(BetaPerm(NodeIndexes(1:n)))))*log(10.0)
318     ENDIF
319     IF (Beta2Formulation) then
320        nodalbetab(1:n)=nodalbetab(1:n)*2.0_dp*BetaValues(BetaPerm(NodeIndexes(1:n)))
321     END IF
322
323     ! Set regularisation to zero for floating points
324     IF (FreeSlipShelves) THEN
325        DO t=1,n
326           IF ( MaskValues(MaskPerm(NodeIndexes(t))).LT.0.0_dp ) THEN
327              NodalRegb(t) = 0.0_dp
328           END IF
329        END DO
330     END IF
331
332     db(NodeIndexes(1:n)) = db(NodeIndexes(1:n)) + nodalbetab(1:n) + NodalRegb(1:n)
333  End do ! on elements
334
335  Do t=1,Solver % Mesh % NumberOfNodes
336     if (VisitedNode(t).lt.1.0_dp) cycle
337     VariableValues(Permutation(t)) = db(t)
338     IF (FreeSlipShelves) THEN
339        IF ( MaskValues(MaskPerm(t)).LT.0.0_dp ) THEN
340           VariableValues(Permutation(t)) = 0.0_dp
341        END IF
342     END IF
343  End do
344
345  IF (FreeSlipShelves) THEN
346     NULLIFY(PointerToMask)
347     NULLIFY(MaskValues)
348     NULLIFY(MaskPerm)
349  END IF
350
351  Return
352
353CONTAINS
354
355  function calcNorm(v) result(v2)
356    implicit none
357    real(kind=dp) :: v(3),v2
358
359    v2=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
360  end function calcNorm
361
362  function scalar(v1,v2) result(vr)
363    implicit none
364    real(kind=dp) :: v2(3),v1(3),vr
365
366    vr=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
367  end function scalar
368
369!------------------------------------------------------------------------------
370END SUBROUTINE DJDBeta_Adjoint
371!------------------------------------------------------------------------------
372
373
374