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