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: Olivier Gagliardini
26! *  Email:   gagliar@lgge.obs.ujf-grenoble.fr
27! *  Web:     http://elmerice.elmerfem.org
28! *
29! *  Original Date: 23 Jully 2009
30! *
31! *****************************************************************************
32!> Solver to inquire the velocity and isotropic pressure from the SIA solution
33!> Exported Variables SIAFlow, dof=dim
34!> Needs to first compute the Depth and the FreeSurfGrad using FlowDepth Solver
35! *****************************************************************************
36SUBROUTINE SIASolver( Model,Solver,dt,TransientSimulation )
37!------------------------------------------------------------------------------
38!******************************************************************************
39!
40!  Solve the SIA Flow solution !
41!
42!  ARGUMENTS:
43!
44!  TYPE(Model_t) :: Model,
45!     INPUT: All model information (mesh, materials, BCs, etc...)
46!
47!  TYPE(Solver_t) :: Solver
48!     INPUT: Linear & nonlinear equation solver options
49!
50!  REAL(KIND=dp) :: dt,
51!     INPUT: Timestep size for time dependent simulations
52!
53!  LOGICAL :: TransientSimulation
54!     INPUT: Steady state or transient simulation
55!
56!******************************************************************************
57  USE DefUtils
58
59  IMPLICIT NONE
60!------------------------------------------------------------------------------
61  TYPE(Solver_t) :: Solver
62  TYPE(Model_t) :: Model
63
64  REAL(KIND=dp) :: dt
65  LOGICAL :: TransientSimulation
66!------------------------------------------------------------------------------
67! Local variables
68!------------------------------------------------------------------------------
69  TYPE(Element_t),POINTER :: CurrentElement, Element
70  TYPE(Matrix_t),POINTER  :: StiffMatrix
71  TYPE(ValueList_t), POINTER :: SolverParams, BodyForce, Material
72  TYPE(Variable_t), POINTER :: PointerToVariable, Grad1Sol, Grad2Sol, &
73                               DepthSol, VeloSol
74
75  LOGICAL :: AllocationsDone = .FALSE., Found, LimitVelocity,UnFoundFatal=.TRUE.
76
77  INTEGER :: i, j, n, m, t, istat, DIM, p, Indexes(128), COMP, constrainedvelocities
78  INTEGER, POINTER :: Permutation(:), VeloPerm(:), &
79       DepthPerm(:), GradSurface1Perm(:), GradSurface2Perm(:), &
80       NodeIndexes(:)
81
82  REAL(KIND=dp), POINTER :: ForceVector(:)
83  REAL(KIND=dp), POINTER :: VariableValues(:), Depth(:), GradSurface1(:), &
84                            GradSurface2(:), Velocity(:), PrevVelo(:,:)
85  REAL(KIND=dp) :: Norm, cn, dd
86
87  REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:), FORCE(:), &
88           NodalGravity(:), NodalViscosity(:), NodalDensity(:), &
89           NodalDepth(:), NodalSurfGrad1(:), NodalSurfGrad2(:), &
90           NodalU(:), NodalV(:)
91  REAL(KIND=dp) :: VelocityCutOff, VeloAbs, VelocityDirection(2)
92
93  CHARACTER(LEN=MAX_NAME_LEN) :: SolverName
94
95
96  SAVE STIFF, LOAD, FORCE, AllocationsDone, DIM, SolverName
97  SAVE NodalGravity, NodalViscosity, NodalDensity, &
98           NodalDepth, NodalSurfGrad1, NodalSurfGrad2, &
99           NodalU, NodalV
100!------------------------------------------------------------------------------
101  PointerToVariable => Solver % Variable
102  Permutation  => PointerToVariable % Perm
103  VariableValues => PointerToVariable % Values
104
105
106  !--------------------------------------------------------------
107  !Allocate some permanent storage, this is done first time only:
108  !--------------------------------------------------------------
109  IF ( (.NOT. AllocationsDone) .OR. Solver % Mesh % Changed  ) THEN
110     WRITE(SolverName, '(A)') 'SIASolver'
111     N = Solver % Mesh % MaxElementNodes ! just big enough for elemental arrays
112     M = Model % Mesh % NumberOfNodes
113     IF (AllocationsDone) DEALLOCATE(FORCE, LOAD, STIFF, NodalGravity, &
114                       NodalViscosity, NodalDensity, NodalDepth, &
115                       NodalSurfGrad1, NodalSurfGrad2, NodalU, NodalV )
116
117     ALLOCATE( FORCE(N), LOAD(N), STIFF(N,N), &
118          NodalGravity(N), NodalDensity(N), NodalViscosity(N), &
119          NodalDepth(N), NodalSurfGrad1(N), NodalSurfGrad2(N), &
120          NodalU(N), NodalV(N), STAT=istat )
121     IF ( istat /= 0 ) THEN
122        CALL Fatal( SolverName, 'Memory allocation error.' )
123     END IF
124
125
126     AllocationsDone = .TRUE.
127     CALL INFO( SolverName, 'Memory allocation done.',Level=1 )
128  END IF
129
130!------------------------------------------------------------------------------
131!    Get variables needed for solution
132!------------------------------------------------------------------------------
133  DIM = CoordinateSystemDimension()
134
135  VeloSol => VariableGet( Solver % Mesh % Variables, 'SIAFlow',UnFoundFatal=UnFoundFatal)
136  Velocity => VeloSol % Values
137  VeloPerm => VeloSol % Perm
138  PrevVelo => veloSol % PrevValues
139
140  DepthSol => VariableGet( Solver % Mesh % Variables, 'Depth',UnFoundFatal=UnFoundFatal)
141  Depth => DepthSol % Values
142  DepthPerm => DepthSol % Perm
143
144  Grad1Sol => VariableGet( Solver % Mesh % Variables, 'FreeSurfGrad1',UnFoundFatal=UnFoundFatal)
145  GradSurface1 => Grad1Sol % Values
146  GradSurface1Perm => Grad1Sol % Perm
147
148  IF (dim > 2) THEN
149     Grad2Sol => VariableGet( Solver % Mesh % Variables, 'FreeSurfGrad2',UnFoundFatal=UnFoundFatal)
150     GradSurface2 => Grad2Sol % Values
151     GradSurface2Perm => Grad2Sol % Perm
152  END IF
153
154  VelocityCutOff = ListGetConstReal( Solver % Values, 'Velocity Cutoff', Found)
155  IF (.NOT.Found) THEN
156     LimitVelocity = .FALSE.
157     CALL INFO(SolverName,'No Velocity Cutoff has been set',Level=1)
158  ELSE
159     LimitVelocity = .TRUE.
160     WRITE(Message, '(A,E10.4)') 'Velocity Cutoff set to ', VelocityCutOff
161     CALL INFO(SolverName,Message, Level=1)
162  END IF
163
164
165     StiffMatrix => Solver % Matrix
166     ForceVector => StiffMatrix % RHS
167
168
169  ! Loop over the velocity components and pressure
170  ! If DIM = 2 u, w, p
171  ! If DIM = 3 u, v, w, p
172  !-----------------------------------------------
173  DO  COMP = 1, DIM+1
174
175! No non-linear iteration, no time dependency
176  VariableValues = 0.0d0
177  Norm = Solver % Variable % Norm
178
179
180
181  !Initialize the system and do the assembly:
182  !------------------------------------------
183  CALL DefaultInitialize()
184  ! bulk assembly
185  DO t=1,Solver % NumberOfActiveElements
186     Element => GetActiveElement(t)
187     IF (ParEnv % myPe .NE. Element % partIndex) CYCLE
188     n = GetElementNOFNodes()
189
190     NodeIndexes => Element % NodeIndexes
191
192     ! Read the gravity in the Body Force Section
193     BodyForce => GetBodyForce()
194     NodalGravity = 0.0_dp
195     IF ( ASSOCIATED( BodyForce ) ) THEN
196           IF (DIM==2) THEN
197           NodalGravity(1:n) = ListGetReal( &
198                   BodyForce, 'Flow BodyForce 2', n, NodeIndexes, Found)
199           ELSE
200           NodalGravity(1:n) = ListGetReal( &
201                   BodyForce, 'Flow BodyForce 3', n, NodeIndexes, Found)
202           END IF
203     END IF
204
205     ! Read the Viscosity eta, density, and exponent m in Material Section
206     ! Same definition as NS Solver in Elmer - n=1/m , A = 1/ (2 eta^n)
207     Material => GetMaterial()
208
209     NodalDensity = 0.0D0
210     NodalDensity(1:n) = ListGetReal( &
211         Material, 'Density', n, NodeIndexes, Found )
212
213     NodalViscosity = 0.0D0
214     NodalViscosity(1:n) = ListGetReal( &
215         Material, 'Viscosity', n, NodeIndexes, Found )
216
217     cn = ListGetConstReal( Material, 'Viscosity Exponent',Found)
218
219     ! Get the Nodal value of Depth, FreeSurfGrad1 and FreeSurfGrad2
220     NodalDepth(1:n) = Depth(DepthPerm(NodeIndexes(1:n)))
221     NodalSurfGrad1(1:n) = GradSurface1(GradSurface1Perm(NodeIndexes(1:n)))
222     NodalSurfGrad2 = 0.0D0
223     IF (DIM==3) NodalSurfGrad2(1:n) = GradSurface2(GradSurface2Perm(NodeIndexes(1:n)))
224
225     IF (COMP==1) THEN     ! u
226        CALL LocalMatrixUV (  STIFF, FORCE, Element, n, NodalGravity, &
227        NodalDensity, NodalViscosity, NodalDepth, NodalSurfGrad1, &
228        NodalSurfGrad2, cn, COMP)
229
230     ELSE IF (COMP==DIM) THEN  ! w
231        NodalU(1:n) = Velocity((DIM+1)*(VeloPerm(NodeIndexes(1:n))-1)+1)
232        NodalV = 0.0D0
233        IF (DIM==3) NodalV(1:n) = Velocity((DIM+1)*(VeloPerm(NodeIndexes(1:n))-1)+2)
234        CALL LocalMatrixW (  STIFF, FORCE, Element, n, NodalU, NodalV )
235
236     ELSE IF (COMP==DIM+1) THEN ! p
237        CALL LocalMatrixP (  STIFF, FORCE, Element, n )
238
239     ELSE               ! v if dim=3
240        CALL LocalMatrixUV (  STIFF, FORCE, Element, n, NodalGravity, &
241        NodalDensity, NodalViscosity, NodalDepth, NodalSurfGrad1, &
242        NodalSurfGrad2, cn, COMP)
243
244     END IF
245
246     CALL DefaultUpdateEquations( STIFF, FORCE )
247  END DO
248
249  ! Neumann conditions only for w and p
250  IF (COMP .GE. DIM) THEN
251  DO t=1,Solver % Mesh % NUmberOfBoundaryElements
252     Element => GetBoundaryElement(t)
253     IF ( GetElementFamily() == 1 ) CYCLE
254     NodeIndexes => Element % NodeIndexes
255     IF (ParEnv % myPe .NE. Element % partIndex) CYCLE
256     n = GetElementNOFNodes()
257     STIFF = 0.0D00
258     FORCE = 0.0D00
259
260     IF (COMP==DIM) THEN
261     ! only for the surface nodes
262        dd = SUM(ABS(Depth(Depthperm(NodeIndexes(1:n)))))
263        IF (dd < 1.0e-6) THEN
264           NodalU(1:n) = Velocity((DIM+1)*(VeloPerm(NodeIndexes(1:n))-1)+1)
265           NodalV = 0.0D0
266           IF (DIM==3) NodalV(1:n) = Velocity((DIM+1)*(VeloPerm(NodeIndexes(1:n))-1)+2)
267           CALL LocalMatrixBCW (  STIFF, FORCE, Element, n, NodalU, NodalV )
268        END IF
269     ELSE IF (COMP==DIM+1) THEN
270            CALL LocalMatrixBCP(  STIFF, FORCE, Element, n, NodalDensity, &
271                    NodalGravity )
272     END IF
273     CALL DefaultUpdateEquations( STIFF, FORCE )
274  END DO
275  END IF
276
277  CALL DefaultFinishAssembly()
278
279  ! Dirichlet
280     CALL SetDirichletBoundaries( Model, StiffMatrix, ForceVector, &
281          ComponentName('SIAFlow',COMP), 1,1, Permutation )
282
283  CALL EnforceDirichletConditions( Solver, Solver % Matrix , Solver % Matrix % RHS)
284
285  !Solve the system
286  Norm = DefaultSolve()
287
288  ! Save the solution on to the right variable
289  constrainedvelocities = 0
290  DO i = 1, Model % Mesh % NumberOfNodes
291     IF (VeloPerm(i)>0) THEN
292        Velocity ((DIM+1)*(VeloPerm(i)-1) + COMP) = VariableValues(Permutation(i))
293        IF (LimitVelocity .AND. (COMP==(DIM-1))) THEN
294           VeloAbs = (Velocity ((DIM+1)*(VeloPerm(i)-1) + COMP))**2
295           IF (COMP==2) VeloAbs = VeloAbs + (Velocity ((DIM+1)*(VeloPerm(i)-1) + COMP - 1))**2
296           VeloAbs = SQRT(VeloAbs)
297           IF (VeloAbs > VelocityCutOff) THEN
298              constrainedvelocities = constrainedvelocities + 1
299              DO j= 1,DIM-1
300                 Velocity ((DIM+1)*(VeloPerm(i)-1) + j)  = VelocityCutOff * Velocity ((DIM+1)*(VeloPerm(i)-1) + j) / VeloAbs
301              END DO
302           END IF
303
304        END IF
305     END IF
306  END DO
307
308END DO ! Loop p
309
310CONTAINS
311
312!------------------------------------------------------------------------------
313  SUBROUTINE LocalMatrixUV(  STIFF, FORCE, Element, n, gravity, &
314           Density, Viscosity, LocalDepth, SurfGrad1, SurfGrad2, cm, cp)
315!------------------------------------------------------------------------------
316    REAL(KIND=dp) :: STIFF(:,:), FORCE(:), gravity(:), Density(:), &
317                     Viscosity(:), LocalDepth(:), SurfGrad1(:), SurfGrad2(:)
318    INTEGER :: n, cp
319    REAL(KIND=dp) :: cm
320    TYPE(Element_t), POINTER :: Element
321!------------------------------------------------------------------------------
322    REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3), detJ
323    REAL(KIND=dp) :: g, rho, eta, d, dhdx, dhdy, dU2dz2, nn, detadz
324    REAL(KIND=dp) :: grad1, grad2
325    LOGICAL :: Stat
326    INTEGER :: t, p, q , dim
327    TYPE(GaussIntegrationPoints_t) :: IP
328
329    TYPE(Nodes_t) :: Nodes
330    SAVE Nodes
331!------------------------------------------------------------------------------
332    CALL GetElementNodes( Nodes )
333    STIFF = 0.0d0
334    FORCE = 0.0d0
335
336    dim = CoordinateSystemDimension()
337
338    nn = 1.0/cm
339
340    IP = GaussPoints( Element )
341    DO t=1,IP % n
342       stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
343        IP % W(t),  detJ, Basis, dBasisdx, ddBasisddx, .FALSE. )
344
345       DO p=1,n
346         DO q=1,n
347           STIFF(p,q) = STIFF(p,q) + IP % S(t) * detJ * dBasisdx(q,dim)*dBasisdx(p,dim)
348         END DO
349       END DO
350       ! Point value
351       g = ABS(SUM( Gravity(1:n) * Basis(1:n) ))
352       rho = SUM( Density(1:n) * Basis(1:n) )
353       eta = SUM( Viscosity(1:n) * Basis(1:n) )
354       detadz = SUM( Viscosity(1:n) * dBasisdx(1:n,dim) )
355       grad1 = SUM( SurfGrad1(1:n) * Basis(1:n) )
356       grad2 = SUM( SurfGrad2(1:n) * Basis(1:n) )
357       d = SUM( LocalDepth(1:n) * Basis(1:n) )
358
359
360       IF (cp==1) dU2dz2 = nn * (rho * g / eta)**nn * grad1 * (1.0 + d * detadz / eta)
361       IF (cp==2) dU2dz2 = nn * (rho * g / eta)**nn * grad2 * (1.0 + d * detadz / eta)
362
363       ! Non linear case
364       IF (nn > 1.0) THEN
365               dU2dz2 = dU2dz2 * (d * SQRT(grad1**2.0 + grad2**2.0))**(nn-1.0)
366       END IF
367
368       FORCE(1:n) = FORCE(1:n) - dU2dz2 * IP % s(t) * detJ * Basis(1:n)
369
370    END DO
371!------------------------------------------------------------------------------
372  END SUBROUTINE LocalMatrixUV
373!------------------------------------------------------------------------------
374
375!------------------------------------------------------------------------------
376  SUBROUTINE LocalMatrixW(  STIFF, FORCE, Element, n, VeloU, VeloV)
377!------------------------------------------------------------------------------
378    REAL(KIND=dp) :: STIFF(:,:), FORCE(:), VeloU(:), VeloV(:)
379    INTEGER :: n
380    TYPE(Element_t), POINTER :: Element
381!------------------------------------------------------------------------------
382    REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3), detJ, &
383                     dU2dxz, dV2dyz
384    LOGICAL :: Stat
385    INTEGER :: t, p,q , DIM
386    TYPE(GaussIntegrationPoints_t) :: IP
387
388    TYPE(Nodes_t) :: Nodes
389    SAVE Nodes
390!------------------------------------------------------------------------------
391    CALL GetElementNodes( Nodes )
392    STIFF = 0.0d0
393    FORCE = 0.0d0
394
395    DIM = CoordinateSystemDimension()
396
397    IP = GaussPoints( Element )
398    DO t=1,IP % n
399       stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
400        IP % W(t),  detJ, Basis, dBasisdx, ddBasisddx, .TRUE. )
401
402       DO p=1,n
403         DO q=1,n
404           STIFF(p,q) = STIFF(p,q) + IP % S(t) * detJ * dBasisdx(q,dim)*dBasisdx(p,dim)
405         END DO
406       END DO
407
408       dU2dxz = SUM(VeloU(1:n)*ddBasisddx(1:n,1,DIM))
409       dV2dyz = 0.0d0
410       IF (DIM==3) dV2dyz = SUM(VeloV(1:n)*ddBasisddx(1:n,2,3))
411
412
413       FORCE(1:n) = FORCE(1:n) + (dU2dxz + dV2dyz) * IP % s(t) * detJ * Basis(1:n)
414
415    END DO
416
417!------------------------------------------------------------------------------
418  END SUBROUTINE LocalMatrixW
419
420!------------------------------------------------------------------------------
421  SUBROUTINE LocalMatrixP(  STIFF, FORCE, Element, n)
422!------------------------------------------------------------------------------
423    REAL(KIND=dp) :: STIFF(:,:), FORCE(:)
424    INTEGER :: n
425    TYPE(Element_t), POINTER :: Element
426!------------------------------------------------------------------------------
427    REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), ddBasisddx(n,3,3), detJ
428    LOGICAL :: Stat
429    INTEGER :: t, p,q ,dim
430    TYPE(GaussIntegrationPoints_t) :: IP
431
432    TYPE(Nodes_t) :: Nodes
433    SAVE Nodes
434!------------------------------------------------------------------------------
435    CALL GetElementNodes( Nodes )
436    STIFF = 0.0d0
437    FORCE = 0.0d0
438
439    dim = CoordinateSystemDimension()
440
441    IP = GaussPoints( Element )
442    DO t=1,IP % n
443       stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
444        IP % W(t),  detJ, Basis, dBasisdx, ddBasisddx, .FALSE. )
445
446       DO p=1,n
447         DO q=1,n
448           STIFF(p,q) = STIFF(p,q) + IP % S(t) * detJ * dBasisdx(q,dim)*dBasisdx(p,dim)
449         END DO
450       END DO
451    END DO
452
453!------------------------------------------------------------------------------
454  END SUBROUTINE LocalMatrixP
455!------------------------------------------------------------------------------
456
457!------------------------------------------------------------------------------
458  SUBROUTINE LocalMatrixBCW(  STIFF, FORCE, Element, n, VeloU, VeloV )
459!------------------------------------------------------------------------------
460    REAL(KIND=dp) :: STIFF(:,:), FORCE(:), veloU(:), veloV(:)
461    INTEGER :: n
462    TYPE(Element_t), POINTER :: Element
463!------------------------------------------------------------------------------
464    REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),ddBasisddx(n,3,3), &
465                      DetJ, Normal(3), grad, dUdx, dVdy
466    LOGICAL :: Stat
467    INTEGER :: t, DIM
468    TYPE(GaussIntegrationPoints_t) :: IP
469
470    TYPE(Nodes_t) :: Nodes
471    SAVE Nodes
472!------------------------------------------------------------------------------
473    CALL GetElementNodes( Nodes )
474    STIFF = 0.0d0
475    FORCE = 0.0d0
476
477    DIM = CoordinateSystemDimension()
478
479    IP = GaussPoints( Element )
480    DO t=1,IP % n
481      stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
482       IP % W(t),  detJ, Basis, dBasisdx, ddBasisddx, .FALSE. )
483
484       dUdx = SUM( VeloU(1:n) * dBasisdx(1:n,1) )
485       dVdy = 0.0e0
486       IF (DIM==3) dVdy = SUM( VeloV(1:n) * dBasisdx(1:n,2) )
487
488       grad = - (dUdx + dVdy)
489
490      Normal = NormalVector( Element, Nodes, IP % U(t), IP % V(t), .TRUE.)
491      FORCE(1:n) = FORCE(1:n) + grad * IP % s(t) * DetJ * Normal(dim) * Basis(1:n)
492    END DO
493!------------------------------------------------------------------------------
494  END SUBROUTINE LocalMatrixBCW
495!------------------------------------------------------------------------------
496
497!------------------------------------------------------------------------------
498  SUBROUTINE LocalMatrixBCP(  STIFF, FORCE, Element, n, Density, &
499                      Gravity)
500!------------------------------------------------------------------------------
501    REAL(KIND=dp) :: STIFF(:,:), FORCE(:), density(:), Gravity(:)
502    INTEGER :: n
503    TYPE(Element_t), POINTER :: Element
504!------------------------------------------------------------------------------
505    REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),ddBasisddx(n,3,3), &
506                      DetJ,Normal(3), rho, g, grad
507    LOGICAL :: Stat
508    INTEGER :: t, dim
509    TYPE(GaussIntegrationPoints_t) :: IP
510
511    TYPE(Nodes_t) :: Nodes
512    SAVE Nodes
513!------------------------------------------------------------------------------
514    CALL GetElementNodes( Nodes )
515    STIFF = 0.0d0
516    FORCE = 0.0d0
517
518    dim = CoordinateSystemDimension()
519
520    IP = GaussPoints( Element )
521    DO t=1,IP % n
522      stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
523       IP % W(t),  detJ, Basis, dBasisdx, ddBasisddx, .FALSE. )
524
525       g = ABS(SUM( Gravity(1:n) * Basis(1:n) ))
526       rho = SUM( Density(1:n) * Basis(1:n) )
527
528       grad = - rho * g
529
530      Normal = NormalVector( Element, Nodes, IP % U(t), IP % V(t), .TRUE.)
531      FORCE(1:n) = FORCE(1:n) + grad * IP % s(t) * DetJ * Normal(dim) * Basis(1:n)
532    END DO
533!------------------------------------------------------------------------------
534  END SUBROUTINE LocalMatrixBCP
535!------------------------------------------------------------------------------
536END SUBROUTINE SIASolver
537!------------------------------------------------------------------------------
538!> Allows to have the SIA Velocity as primary variables (and not exported one)
539!> Allow then to have access to the Previous values to construct a stable
540!> time discretization scheme.
541SUBROUTINE SIAVariable( Model,Solver,dt,TransientSimulation )
542!------------------------------------------------------------------------------
543
544  USE DefUtils
545  IMPLICIT NONE
546!------------------------------------------------------------------------------
547  TYPE(Solver_t) :: Solver
548  TYPE(Model_t) :: Model
549
550  REAL(KIND=dp) :: dt
551  LOGICAL :: TransientSimulation
552!------------------------------------------------------------------------------
553END SUBROUTINE SIAVariable
554!------------------------------------------------------------------------------
555
556
557