CoCalc -- AdjointThickness_ThicknessSolver.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/Solvers/AdjointThickness/AdjointThickness_ThicknessSolver.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
! ******************************************************************************
24
! *
25
! *  Authors: F. Gillet-Chaulet
26
! *  Web:     http://elmerice.elmerfem.org
27
! *
28
! *  Original Date: Dec. 2020
29
! * 
30
! *****************************************************************************
31
!-----------------------------------------------------------------------------
32
SUBROUTINE AdjointThickness_ThicknessSolver( Model,Solver,dt,TransientSimulation )
33
  USE DefUtils
34
  IMPLICIT NONE
35
  !------------------------------------------------------------------------------
36
  !    external variables
37
  !------------------------------------------------------------------------------
38
  TYPE(Model_t) :: Model
39
  TYPE(Solver_t):: Solver
40
  REAL(KIND=dp) :: dt
41
  LOGICAL :: TransientSimulation
42
  !------------------------------------------------------------------------------
43
  !    Local variables
44
  !------------------------------------------------------------------------------
45
  TYPE(ValueList_t), POINTER :: BodyForce, SolverParams
46

47
  INTEGER :: NMAX,NSDOFs
48
  INTEGER :: istat
49
  INTEGER :: t,i,j,n
50
  LOGICAL :: ConvectionVar, Found
51

52
  REAL(KIND=dp) :: Norm
53

54
  CHARACTER(LEN=MAX_NAME_LEN)  :: FlowSolName
55
  TYPE(Variable_t), POINTER :: FlowSol
56
  REAL(KIND=dp), POINTER :: FlowSolution(:)
57
  INTEGER, POINTER :: FlowPerm(:)
58

59
  REAL(KIND=dp), ALLOCATABLE,SAVE :: STIFF(:,:),FORCE(:),LOAD(:),Velo(:,:)
60

61
  CHARACTER(LEN=MAX_NAME_LEN)  :: SolverName= 'AdjointThickness_ThicknessSolver'
62

63
  TYPE(Nodes_t),SAVE   :: ElementNodes
64
  TYPE(Element_t),POINTER :: CurrentElement
65
  INTEGER, POINTER :: NodeIndexes(:)
66

67
  LOGICAL, SAVE :: AllocationsDone = .FALSE.
68

69
  !------------------------------------------------------------------------------
70
  !    Go
71
  !------------------------------------------------------------------------------
72
  SolverParams => GetSolverParams()
73

74
  !------------------------------------------------------------------------------
75
  !    check incompatibilities
76
  !------------------------------------------------------------------------------
77
  IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
78
          CALL FATAL(SolverName,'sorry only for cartesian coordinate system')
79
  END IF
80
  IF (TransientSimulation) THEN
81
          CALL FATAL(SolverName,'sorry works only in steady state')
82
  ENDIF
83

84
  !------------------------------------------------------------------------------
85
  !    Allocate some permanent storage, this is done first time only
86
  !------------------------------------------------------------------------------
87
  IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed) THEN
88
     NMAX = Model % MaxElementNodes
89

90
     IF ( AllocationsDone ) THEN
91
        DEALLOCATE( ElementNodes % x,    &
92
             ElementNodes % y,    &
93
             ElementNodes % z,    &
94
             FORCE,    &
95
             STIFF, &
96
             LOAD,&
97
             Velo)
98
     END IF
99

100
     ALLOCATE( ElementNodes % x( NMAX ),    &
101
          ElementNodes % y( NMAX ),    &
102
          ElementNodes % z( NMAX ),    &
103
          FORCE( NMAX ),    &
104
          STIFF( NMAX, NMAX ), &
105
          LOAD(NMAX) , &
106
          Velo( 2, NMAX ), &
107
          STAT=istat )
108
     IF ( istat /= 0 ) THEN
109
        CALL Fatal(SolverName,'Memory allocation error 1, Aborting.')
110
     END IF
111

112
     CALL Info(SolverName,'Memory allocations done' )
113
     AllocationsDone = .TRUE.
114
  END IF
115

116
  !------------------------------------------------------------------------------
117
  !    Get Flow solution
118
  !------------------------------------------------------------------------------
119
   ConvectionVar=.True.
120
   FlowSolName =  GetString(SolverParams ,'Flow Solution Name', Found)
121
   IF(.NOT.Found) THEN        
122
        WRITE(Message,'(A)') &
123
           '<Flow Solution Name> Not Found; will look for <convection velocity> in body forces'
124
        CALL Info(SolverName,Message,level=10)
125
        ConvectionVar=.False.
126
        NSDOFS=GetInteger(SolverParams ,'Convection Dimension',Found)
127
        IF(.NOT.Found) &
128
            CALL Fatal(SolverName,'if <Flow Solution Name> not given prescribe <Convection Dimension>')
129
   ELSE
130
        FlowSol => VariableGet( Solver % Mesh % Variables, FlowSolName,UnFoundFatal=.TRUE.)
131
        FlowPerm     => FlowSol % Perm
132
        NSDOFs     =  FlowSol % DOFs
133
        FlowSolution => FlowSol % Values
134
   END IF
135

136

137
   CALL DefaultInitialize()
138
   !------------------------------------------------------------------------------
139
   !    Do the assembly
140
   !------------------------------------------------------------------------------
141
   DO t=1,Solver % NumberOfActiveElements
142

143
        CurrentElement => GetActiveElement(t)
144
        n = GetElementNOFNodes(CurrentElement)
145
        NodeIndexes => CurrentElement % NodeIndexes
146

147
        ! set coords of highest occurring dimension to zero (to get correct path element)
148
        !-------------------------------------------------------------------------------
149
        ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
150
        IF (NSDOFs == 1) THEN
151
           ElementNodes % y(1:n) = 0.0
152
           ElementNodes % z(1:n) = 0.0
153
        ELSE IF (NSDOFs == 2) THEN
154
           ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
155
           ElementNodes % z(1:n) = 0.0_dp
156
        ELSE
157
           WRITE(Message,'(a,i0,a)')&
158
                'It is not possible to compute Thickness evolution if Flow Sol DOFs=',&
159
                NSDOFs, ' . Aborting'
160
           CALL Fatal( SolverName, Message) 
161
           STOP   
162
        END IF
163

164
        ! get pointers on BF
165
        !---------------------------------------------------------------
166
        BodyForce => GetBodyForce(CurrentElement)
167

168
        ! get flow soulution and velocity field from it
169
        !----------------------------------------------
170
        Velo = 0.0_dp
171
        !----------------------------------------------------
172
        ! get velocity profile
173
        IF (ConvectionVar) Then
174
          DO i=1,n
175
             j = NSDOFs*FlowPerm(NodeIndexes(i))
176
              !2D problem - 1D Thickness evolution
177
              IF(NSDOFs == 1) THEN 
178
                 Velo(1,i) = FlowSolution( j ) 
179
                 Velo(2,i) = 0.0_dp
180
              !2D problem 
181
              ELSE IF (NSDOFs == 2) THEN
182
                 Velo(1,i) = FlowSolution( j-1 ) 
183
                 Velo(2,i) = FlowSolution( j ) 
184
              ELSE
185
                 WRITE(Message,'(a)') 'Velocity should be size 1 or 2'
186
                 CALL Fatal( SolverName, Message)
187
              END IF
188
           END DO
189
       ELSE
190
        IF (ASSOCIATED( BodyForce ) ) THEN
191
          Velo(1,1:n) = ListGetReal( BodyForce, 'Convection Velocity 1',n, NodeIndexes,UnFoundFatal=.TRUE. )
192
          if (NSDOFs.eq.2) &
193
           Velo(2,1:n) = ListGetReal( BodyForce, 'Convection Velocity 2',n, NodeIndexes,UnFoundFatal=.TRUE. )
194
        END IF
195
       END IF
196
        !------------------------------------------------------------------------------
197
        !      get the accumulation/ablation rate (i.e. normal surface flux)
198
        !      from the body force section
199
        !------------------------------------------------------------------------------
200
        LOAD=0.0_dp
201
        IF (ASSOCIATED( BodyForce ) ) THEN
202
              LOAD(1:n) = LOAD(1:n) +   &
203
                      GetReal( BodyForce, 'Top Surface Accumulation', Found )
204
              LOAD(1:n) = LOAD(1:n) +   &
205
                      GetReal( BodyForce, 'Bottom Surface Accumulation', Found )
206
        END IF
207

208
        !------------------------------------------------------------------------------
209
        !      Get element local matrix, and rhs vector
210
        !------------------------------------------------------------------------------
211
        CALL LocalMatrix( STIFF, FORCE,LOAD,&
212
                Velo, NSDOFs, &
213
             CurrentElement, n, ElementNodes, NodeIndexes)
214

215
        !------------------------------------------------------------------------------
216
        CALL DefaultUpdateEquations( STIFF, FORCE )
217
        !------------------------------------------------------------------------------
218
        
219
     END DO ! End loop bulk elements
220
     CALL DefaultFinishBulkAssembly()
221

222
     !------------------------------------------------------------------------------
223
     !     Neumann & Newton boundary conditions
224
     !------------------------------------------------------------------------------
225
     !
226
     ! MIND: In weak formulation it is not possible to prescribe a contact angle on
227
     !       a boundary in this solver. This has to be taken care of in the boundary
228
     !       condition for the stress tensor in the Navier-Stokes Solver. Thus, in
229
     !       generally it does not make sense to prescribe a Neumann type of
230
     !       condition here.
231

232
     !------------------------------------------------------------------------------
233
     !    FinishAssemebly must be called after all other assembly steps, but before
234
     !    Dirichlet boundary settings. Actually no need to call it except for
235
     !    transient simulations.
236
     !------------------------------------------------------------------------------
237
     CALL DefaultFinishAssembly()
238
     CALL DefaultDirichletBCs()
239

240
     Norm = DefaultSolve()
241

242
     !------------------------------------------------------------------------------
243
   CONTAINS
244

245
     !------------------------------------------------------------------------------
246
     !==============================================================================
247
     SUBROUTINE LocalMatrix( STIFF, FORCE,LOAD,&
248
                     Velo, NSDOFs, &
249
                     Element, n, Nodes, NodeIndexes)
250
       !------------------------------------------------------------------------------
251
       !    INPUT:  LOAD(:)   nodal values of the accumulation/ablation function
252
       !            
253
       !            Element         current element
254
       !            n               number of nodes
255
       !            Nodes           current node points
256
       !
257
       !    OUTPUT: STIFF(:,:)
258
       !            FORCE(:)
259
       !------------------------------------------------------------------------------
260
       !      external variables:
261
       !      ------------------------------------------------------------------------
262
       REAL(KIND=dp) :: STIFF(:,:),FORCE(:), LOAD(:), Velo(:,:)
263

264
       INTEGER :: n, NodeIndexes(:), NSDOFs
265
       TYPE(Nodes_t) :: Nodes
266
       TYPE(Element_t), POINTER :: Element
267

268
       !------------------------------------------------------------------------------
269
       !      internal variables:
270
       !      ------------------------------------------------------------------------
271
       TYPE(GaussIntegrationPoints_t) :: IntegStuff
272
       REAL(KIND=dp) :: Basis(n),dBasisdx(n,3)
273
       REAL(KIND=dp) :: SU(n),SW(n)
274
       REAL(KIND=dp) :: Vgauss(2), UNorm,divu,Source
275
       REAL(KIND=dp) :: U,V,W,S,SqrtElementMetric
276
       REAL(KIND=dp) :: Tau,hk
277
       INTEGER :: i,j,t,p,q
278
       LOGICAL :: stat
279
       !------------------------------------------------------------------------------
280

281
       FORCE = 0.0_dp
282
       STIFF = 0.0_dp
283

284
       hK = ElementDiameter( Element, Nodes )
285

286
       !
287
       !      Numerical integration:
288
       !      ----------------------
289
       IntegStuff = GaussPoints( Element )
290

291
       SU = 0.0_dp
292
       SW = 0.0_dp
293

294
       DO t = 1,IntegStuff % n
295
          U = IntegStuff % u(t)
296
          V = IntegStuff % v(t)
297
          W = IntegStuff % w(t)
298
          S = IntegStuff % s(t)
299
          !
300
          !        Basis function values & derivatives at the integration point:
301
          !        -------------------------------------------------------------
302
          stat = ElementInfo( Element,Nodes,U,V,W,SqrtElementMetric, &
303
               Basis,dBasisdx, Bubbles=.False. )
304

305
          !        Correction from metric
306
          !        ----------------------
307
          S = S * SqrtElementMetric
308

309
          !
310
          !        Velocities and (norm of) gradient of free surface and source function 
311
          !        at Gauss point
312
          !        ---------------------------------------------------------------------
313

314
          Vgauss=0.0_dp
315
          DO i=1,NSDOFs
316
             Vgauss(i) = SUM( Basis(1:n)*(Velo(i,1:n)))
317
          END DO
318

319
          divu = 0.0_dp
320
          DO i=1,NSDOFs
321
             divu = divu +  SUM( dBasisdx(1:n,i)*(Velo(i,1:n)))
322
          END DO
323

324
          UNorm = SQRT( SUM( Vgauss(1:NSDOFs)**2 ) )
325
          Tau=0._dp
326
          IF (UNorm .GT. 0.0_dp) Tau = hK / ( 2*Unorm )
327

328
          DO p=1,n
329
             SU(p) = 0.0_dp
330
             SW(p) = 0.0_dp
331
             DO i=1,NSDOFs
332
                   SU(p) = SU(p) + Vgauss(i) * dBasisdx(p,i)
333
                   SW(p) = SW(p) + Vgauss(i) * dBasisdx(p,i)
334
             END DO
335
          END DO
336

337
          !        Stiffness matrix:
338
          !        -----------------
339
          DO p=1,n
340
             DO q=1,n
341
                DO i=1,NSDOFs
342
                   STIFF(p,q) = STIFF(p,q) + &
343
                        s * Vgauss(i) * dBasisdx(q,i) * Basis(p)
344
                END DO
345
                STIFF(p,q) =  STIFF(p,q) + s * Tau * SU(q) * SW(p)
346
                STIFF(p,q) =  STIFF(p,q) + s * divu * Basis(q) * (Basis(p) + Tau*SW(p))
347
             END DO
348
          END DO
349

350
          !        Get accumulation/ablation function 
351
          !        --------------------------------------------------------- 
352
          Source = 0.0_dp
353
          Source=SUM(Basis(1:n)*LOAD(1:n))
354

355
          !        Assemble force vector:
356
          !        ---------------------
357
          FORCE(1:n) = FORCE(1:n) &
358
               + Source * (Basis(1:n) + Tau*SW(1:n)) * s
359

360
       END DO
361

362
       !------------------------------------------------------------------------------
363
     END SUBROUTINE LocalMatrix
364

365
     !------------------------------------------------------------------------------
366
   END SUBROUTINE AdjointThickness_ThicknessSolver
367
!------------------------------------------------------------------------------
368

369
Product

Resources

Company
