CoCalc -- AdjointSSA_CostTaubSolver.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/Solvers/AdjointSSA/AdjointSSA_CostTaubSolver.F90
³²⁰⁶ views
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
SUBROUTINE AdjointSSA_CostTaubSolver_init0(Model,Solver,dt,TransientSimulation )
24
!------------------------------------------------------------------------------
25
  USE DefUtils
26
  IMPLICIT NONE
27
!------------------------------------------------------------------------------
28
  TYPE(Solver_t), TARGET :: Solver
29
  TYPE(Model_t) :: Model
30
  REAL(KIND=dp) :: dt
31
  LOGICAL :: TransientSimulation
32
!------------------------------------------------------------------------------
33
! Local variables
34
!------------------------------------------------------------------------------
35
  CHARACTER(LEN=MAX_NAME_LEN) :: Name
36

37
  Name = ListGetString( Solver % Values, 'Equation',UnFoundFatal=.TRUE.)
38
  CALL ListAddNewString( Solver % Values,'Variable',&
39
          '-nooutput '//TRIM(Name)//'_var')
40
  CALL ListAddLogical(Solver % Values, 'Optimize Bandwidth',.FALSE.)
41

42
END SUBROUTINE AdjointSSA_CostTaubSolver_init0
43
! ******************************************************************************
44
! *****************************************************************************
45
SUBROUTINE AdjointSSA_CostTaubSolver( Model,Solver,dt,TransientSimulation )
46
! *****************************************************************************
47
!------------------------------------------------------------------------------
48
!
49
!  Compute a cost function that penalises 1rst derivative of the basal shear stress
50
!     J=int_{Pb dimension} 0.5 * ( dTau_b/dx )**2
51
!
52
!     The basal friction is computed at the nodes following:
53
!       Tau_b=beta Velocity_nom**fm
54
!          with fm the friction exponent
55
!
56
!   and provides the derivatives with respect to beta en velocity
57
!
58
!   Be Careful, by default, this solver 
59
!      - reset derivatives wrt beta to 0 (Set Reset DJDBeta = Logical False,
60
!                   if values have been computed in a previous cost function)
61
!      - do not reset Cost and derivatives wrt velocity to 0 (Set Reset Cost value = Logical True,
62
!                    if values have been computed in a previous cost function)
63
!
64
!     OUTPUT are : J ; DJDBeta ; Velocityb
65
!
66
!     INPUT PARAMETERS are:
67
!
68
!      In solver section:
69
!               Reset Cost Value = Logical (default = .FALSE.)
70
!               Reset DJDBeta = Logical (default = .True.)
71
!               Cost Filename = File (default = 'CostTaub.dat')
72
!               Cost Variable Name = String (default= 'CostValue')
73
!               DJDBeta Name = String (default= 'DJDBeta')
74
!               Lambda = Real (default = 1.0)
75
!
76
!
77
!      Variables
78
!                SSAVelocity (solution of the SSA pb)
79
!                Velocityb (forcing for the adjoint pb)
80
!
81
!     In Material:
82
!       Keywords related to SSA Friction law (only linear and Weertman)
83
!------------------------------------------------------------------------------
84
!******************************************************************************
85
  USE DefUtils
86
  IMPLICIT NONE
87
!------------------------------------------------------------------------------
88
  TYPE(Solver_t) :: Solver
89
  TYPE(Model_t) :: Model
90

91
  REAL(KIND=dp) :: dt
92
  LOGICAL :: TransientSimulation
93
!  
94
  TYPE(Variable_t), POINTER :: CostVar
95
  TYPE(Variable_t), POINTER :: TimeVar
96
  TYPE(Variable_t), POINTER :: DJDVariable !derivée/beta
97
  TYPE(Variable_t), POINTER :: VVar ! ssavelocity
98
  TYPE(Variable_t), POINTER :: VbVar ! velocityb
99
  REAL(KIND=dp),POINTER :: DJDValues(:),VValues(:),VbValues(:)
100
  INTEGER,POINTER :: DJDPerm(:),VPerm(:),VbPerm(:)
101
  INTEGER :: DOFs 
102
  TYPE(Nodes_t),SAVE :: ElementNodes
103
  INTEGER :: N
104
  TYPE(ValueList_t), POINTER :: SolverParams
105
  TYPE(ValueList_t), POINTER :: Material
106
  TYPE(Element_t), POINTER :: Element
107
  INTEGER, POINTER :: NodeIndexes(:)
108
  TYPE(GaussIntegrationPoints_t) :: IntegStuff
109
  REAL(KIND=dp) :: U,V,W,SqrtElementMetric
110
  REAL(KIND=dp),ALLOCATABLE,SAVE :: Basis(:), dBasisdx(:,:)
111
  REAL(KIND=dp),ALLOCATABLE,SAVE :: Beta(:),Vnode(:,:),Vn(:),Taub(:),NodalBetaDer(:)
112
  REAL(KIND=dp),ALLOCATABLE,SAVE :: Betab(:),Vnodeb(:,:),Vnb(:),Taubb(:)
113
  REAL(KIND=dp) :: Cost,Cost_S
114
  REAL(KIND=dp) :: coeff,coeffb,s,Lambda
115
  REAL(KIND=dp) :: fm
116

117
  CHARACTER(LEN=MAX_NAME_LEN),SAVE :: CostSolName
118
  CHARACTER(LEN=MAX_NAME_LEN),SAVE :: CostFile
119
  CHARACTER(LEN=MAX_NAME_LEN) :: DefaultCostFile="CostTaub.dat"
120
  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName="CostTaub"
121
  CHARACTER(LEN=MAX_NAME_LEN) :: SName
122
  CHARACTER(LEN=MAX_NAME_LEN) :: Friction
123
  CHARACTER*10 :: date,temps
124

125
  LOGICAL :: Reset 
126
  LOGICAL :: ResetCost
127
  LOGICAL :: Found
128
  LOGICAL :: stat
129
  LOGICAL :: HaveBetaDer
130
  LOGICAL, SAVE :: Parallel
131
  LOGICAL, SAVE :: Firsttime=.TRUE.
132

133
  INTEGER :: i,t,p,j
134
  INTEGER :: ierr
135

136

137
  SolverParams => GetSolverParams()
138

139
  Lambda =  GetConstReal( SolverParams,'Lambda', Found)
140
  IF(.NOT.Found) THEN
141
     CALL WARN(SolverName,'Keyword >Lambda< not found  in section >Solver<')
142
     CALL WARN(SolverName,'Taking default value Lambda=1.0')
143
     Lambda = 1.0
144
  End if
145

146
  ResetCost =  GetLogical( SolverParams,'Reset Cost Value', Found)
147
  IF(.NOT.Found) ResetCost=.FALSE.
148

149
!!! SOME INITIALISATION AT FIRST TIME
150
  If (Firsttime) then
151
    N = model % MaxElementNodes
152
    allocate(ElementNodes % x(N), ElementNodes % y(N), ElementNodes % z(N))
153
    allocate( Basis(N),dBasisdx(N,3))
154
    allocate( Beta(N),Vnode(N,3),Vn(N),Taub(N),NodalBetaDer(N))
155
    allocate( Betab(N),Vnodeb(N,3),Vnb(N),Taubb(N))
156

157
!!!!!!! Check for parallel run 
158
    Parallel = .FALSE.
159
    IF ( ASSOCIATED( Solver % Matrix % ParMatrix ) ) THEN
160
            IF ( Solver %  Matrix % ParMatrix % ParEnv % PEs > 1 )  THEN
161
                    Parallel = .TRUE.
162
            END IF
163
    END IF
164

165
    !! check some kws
166
    SName =  GetString( SolverParams,'DJDBeta Name', Found)
167
    IF(.NOT.Found) THEN
168
          CALL WARN(SolverName,'Keyword >DJDBeta Name< not found in section >Solver<')
169
          CALL WARN(SolverName,'Taking default value >DJDBeta<')
170
          WRITE(SName,'(A)') 'DJDBeta'
171
          CALL ListAddString(  SolverParams, 'DJDBeta Name', TRIM(SName))
172
    END IF
173
    !!
174
    CostSolName =  GetString( SolverParams,'Cost Variable Name', Found)
175
    IF(.NOT.Found) THEN
176
       CALL WARN(SolverName,'Keyword >Cost Variable Name< not found  in section >Solver<')
177
       CALL WARN(SolverName,'Taking default value >CostValue<')
178
       WRITE(CostSolName,'(A)') 'CostValue'
179
    END IF
180

181

182
    CostFile = ListGetString(Solver % Values,'Cost Filename',Found )
183
    IF (.NOT. Found) CostFile = DefaultCostFile
184
    CALL DATE_AND_TIME(date,temps)
185
    If (Parallel) then
186
        if (ParEnv % MyPe.EQ.0) then
187
           OPEN (12, FILE=CostFile)
188
                   write(12,1000) date(5:6),date(7:8),date(1:4),temps(1:2),temps(3:4),temps(5:6)
189
                   write(12,1001) Lambda
190
                   write(12,'(A)') '# iter, Jreg'
191
           CLOSE(12)
192
         End if
193
    Else
194
           OPEN (12, FILE=CostFile)
195
                   write(12,1000) date(5:6),date(7:8),date(1:4),temps(1:2),temps(3:4),temps(5:6)
196
                   write(12,1001) Lambda
197
                   write(12,'(A)') '# iter, Jreg'
198
           CLOSE(12)
199
    End if
200
   
201
  !!! End of First visit
202
    Firsttime=.false.
203
  Endif
204
!!!! INITIALISATION DONE
205

206
    SName =  ListGetString( SolverParams,'DJDBeta Name', UnFoundFatal=.TRUE.)
207
    DJDVariable => VariableGet( Solver % Mesh % Variables,TRIM(SName),UnFoundFatal=.TRUE.)
208
    DJDValues => DJDVariable % Values
209
    DJDPerm => DJDVariable % Perm
210
    Reset =  ListGetLogical( SolverParams,'Reset DJDBeta', Found)
211
    if (Reset.OR.(.NOT.Found)) DJDValues=0.0_dp ! reset to zero herz!
212
    
213

214
    VVar => VariableGet( Solver % Mesh % Variables, 'ssavelocity',UnFoundFatal=.TRUE.)
215
    VValues => VVar % Values
216
    VPerm => VVar % Perm
217
    DOFs = VVar % DOFs
218

219
    VbVar => VariableGet( Solver % Mesh % Variables, 'velocityb',UnFoundFatal=.TRUE.)
220
    VbValues => VbVar % Values
221
    VbPerm => VbVar % Perm
222
    IF (ResetCost) VbValues = 0._dp
223
    IF (VbVar%DOFs.NE.DOFs) CALL FATAL(SolverName,'Dimension error')
224

225
    Cost=0._dp
226

227
    DO t=1,Solver % NumberOfActiveElements
228
       Element => GetActiveElement(t)
229
       IF (CheckPassiveElement(Element)) THEN
230
          CYCLE
231
       END IF
232
       n = GetElementNOFNodes(Element)
233

234
       NodeIndexes => Element % NodeIndexes
235

236
! set coords of highest occurring dimension to zero (to get correct path element)
237
       !-------------------------------------------------------------------------------
238
       ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes(1:n))
239
       IF (DOFs.EQ.2 ) THEN
240
         ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes(1:n))
241
       ELSE
242
         ElementNodes % y(1:n) = 0.0_dp
243
       ENDIF
244
       ElementNodes % z(1:n) = 0.0_dp
245

246
      ! Get Friction coefficient
247
       Material => GetMaterial(Element)
248
       Friction = ListGetString(Material, 'SSA Friction Law', Found,UnFoundFatal=.TRUE.)
249

250
       SELECT CASE(Friction)
251
         CASE('linear')
252
          fm = 1.0_dp
253
         CASE('weertman')
254
          fm = ListGetConstReal( Material, 'SSA Friction Exponent', UnFoundFatal=.TRUE.)
255
         CASE DEFAULT
256
           CALL FATAL(SolverName,'Friction should be linear or Weertman')
257
      END SELECT
258
      
259
       Beta(1:n) = ListGetReal( Material, 'SSA Friction Parameter', n, NodeIndexes,UnFoundFatal=.TRUE.)
260
       NodalBetaDer(1:n) = ListGetReal( Material, 'SSA Friction Parameter Derivative',n, NodeIndexes,Found=HaveBetaDer)
261
       DO i=1,n
262
         Vn(i)=0._dp
263
         Do j=1,DOFs
264
           Vnode(i,j)=VValues(VVar%DOFs*(VPerm(NodeIndexes(i))-1)+j)
265
           Vn(i)=Vn(i)+Vnode(i,j)*Vnode(i,j)
266
         END DO
267
         IF (Vn(i).GT.AEPS) THEN
268
           Taub(i)=Beta(i)*Vn(i)**(fm/2)
269
         ELSE
270
           Taub(i)=0._dp
271
         END IF
272
       END DO
273

274
!------------------------------------------------------------------------------
275
!    Numerical integration
276
!------------------------------------------------------------------------------
277
        IntegStuff = GaussPoints( Element )
278

279
        Taubb=0._dp
280
        DO p=1,IntegStuff % n
281
          U = IntegStuff % u(p)
282
          V = IntegStuff % v(p)
283
          W = IntegStuff % w(p)
284
!------------------------------------------------------------------------------
285
!        Basis function values & derivatives at the integration point
286
!------------------------------------------------------------------------------
287
          stat = ElementInfo( Element,ElementNodes,U,V,W,SqrtElementMetric, &
288
              Basis,dBasisdx )
289

290
          s =  SqrtElementMetric * IntegStuff % s(p)
291

292
          DO j=1,DOFs
293
           coeff=0._dp
294
           DO i=1,n
295
            coeff=coeff+Taub(i)*dbasisdx(i,j)
296
           END DO
297
           Cost=Cost+0.5*s*coeff*coeff
298

299
           coeffb=s*Lambda*coeff
300
           DO i=1,n
301
            Taubb(i)=Taubb(i)+coeffb*dbasisdx(i,j)
302
           END DO
303

304
          END DO
305

306
        End do !IP
307

308
        DO i=1,n
309
          IF (Vn(i).GT.AEPS) THEN
310
            Vnb(i)=Taubb(i)*Beta(i)*0.5*fm*Vn(i)**(fm/2-1._dp)
311
            Betab(i)=Taubb(i)*Vn(i)**(fm/2)
312
            IF (HaveBetaDer) Betab(i)=Betab(i)*NodalBetaDer(i)
313
          ELSE
314
            Vnb(i)=0._dp
315
            Betab(i)=0._dp
316
          END IF
317

318
          DO j=1,DOFs
319
            Vnodeb(i,j)=Vnb(i)*2*Vnode(i,j)
320
            VbValues(VbVar%DOFs*(VbPerm(NodeIndexes(i))-1)+j)=&
321
              VbValues(VbVar%DOFs*(VbPerm(NodeIndexes(i))-1)+j)+Vnodeb(i,j)
322
          END DO
323

324
          DJDValues(DJDPerm(NodeIndexes(i)))=&
325
              DJDValues(DJDPerm(NodeIndexes(i)))+Betab(i)
326
        END DO
327

328
    End do !Elements
329

330
    TimeVar => VariableGet( Solver % Mesh % Variables, 'Time' )
331
    IF (Parallel) THEN
332
         CALL MPI_ALLREDUCE(Cost,Cost_S,1,&
333
                  MPI_DOUBLE_PRECISION,MPI_SUM,ELMER_COMM_WORLD,ierr)
334
         CostVar => VariableGet( Solver % Mesh % Variables, TRIM(CostSolName) ,UnFoundFatal=.TRUE.)
335
         IF (ResetCost) then
336
             CostVar % Values(1)=Lambda*Cost_S
337
         Else
338
             CostVar % Values(1)=CostVar % Values(1)+Lambda*Cost_S
339
         Endif
340
         IF (Solver % Matrix % ParMatrix % ParEnv % MyPE == 0) then
341
            OPEN (12, FILE=TRIM(CostFile),POSITION='APPEND')
342
              WRITE(12,'(e13.5,2x,e15.8)') TimeVar % Values(1),Cost_S
343
            CLOSE(12)
344
         End if
345
   ELSE
346
         CostVar => VariableGet( Solver % Mesh % Variables, TRIM(CostSolName),UnFoundFatal=.TRUE. )
347
         IF (ResetCost) then
348
              CostVar % Values(1)=Lambda*Cost
349
         Else
350
              CostVar % Values(1)=CostVar % Values(1)+Lambda*Cost
351
         Endif
352
         OPEN (12, FILE=TRIM(CostFile),POSITION='APPEND')
353
           WRITE(12,'(e13.5,2x,e15.8)') TimeVar % Values(1),Cost
354
         CLOSE(12)
355
   END IF
356
   
357
   RETURN
358

359
 1000  format('#date,time,',a2,'/',a2,'/',a4,',',a2,':',a2,':',a2)
360
 1001  format('#lambda,',e15.8)
361
!------------------------------------------------------------------------------
362
END SUBROUTINE AdjointSSA_CostTaubSolver
363
!------------------------------------------------------------------------------
364
! *****************************************************************************
365

366
Product

Resources

Company
