!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Test GlaDS solvers !!
!! !!
!! Olivier Gagliardini June 2017 !!
!! !!
!! SHMIP - test B5 !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
check keywords warn
echo on
$Bnum = "B5"
$namerun = "B5_shmip_"
! All Units are in m, year, MPa
! For ice flow
$yearindays = 365.00
$yearinsec = yearindays*24*60*60
$MPainPa = 1.0e6
$Iter = 2
$dtIni = (1.0/yearindays)
$dtMin = 0.01/yearindays
$OutPut = 1
$ev = 0.0
$Source = 0.0 ! all input in moulins
! Common input parameters for all tests SHMIP
$ub = 1.0e-6*yearinsec !(1.0e-6 m/s)
!Prefactor from Cuffey and Paterson (2010) in MPa^{-3} a^{-1}
! temperate Ice (power law for SSA)
$rhoi = 910.0/(MPainPa*yearinsec^2)
$rhow = 1000.0/(MPainPa*yearinsec^2)
$A1 = 2.89165e-13*yearinsec*MPainPa^3
$A2 = 2.42736e-02*yearinsec*MPainPa^3
$ng = 3.0
$Aglen = 2.5e-25*yearinsec*MPainPa^3
$eta = (2.0*Aglen)^(-1.0/ng)
$gravity = 9.80*yearinsec^2
! For the sheet
$Ar = Aglen
$alphas = 1.25
$betas = 1.5
$lr = 2.0
$hr = 0.1
$Ks = 0.005*yearinsec*(1.0/MPainPa)^(1.0-betas)
$Hs = 0.05 ! IC for h
! For the Channels
$alphac = 1.25
$betac = 1.5
$Kc = 0.1*yearinsec*(1.0/MPainPa)^(1.0-betac)
$Ac = Aglen
$lc = 2.0
$Ct = -7.5e-8*MPainPa
$Cw = 4220.0*yearinsec^2
$Lw = 334000.0*yearinsec^2
! For the Moulins
$Am = 4.0
! Definition of the sqrt geometry
! for SHMIP A, B, C, D
$Hmin = 1.0
$ function H(x) \
import Hmin {\
_H = 6.0*(sqrt(x+5e3)-sqrt(5e3))+Hmin ;\
}
Header
Mesh DB "." "mesh_B5"
End
!---------------------------------------------------
!---------------- CONSTANTS ------------------------
!---------------------------------------------------
Constants
Latent Heat = Real $Lw
Gravity Norm = Real $gravity
Fresh Water Density = Real $rhow
Ice Density = Real $rhoi
Sheet Thickness Variable Name = String "Sheet Thickness"
Hydraulic Potential Variable Name = String "Hydraulic Potential"
Channel Area Variable Name = String "Channel Area"
Bedrock Variable Name = String "Zb"
End
!---------------------------------------------------
!---------------- SIMULATION -----------------------
!---------------------------------------------------
Simulation
Coordinate System = Cartesian 2D
Simulation Type = transient
Timestepping Method = "bdf"
BDF Order = 1
Timestep Intervals(1) = $Iter
Output Intervals(1) = $OutPut
Timestep Sizes(1) = $dtIni
Steady State Max Iterations = 1
Steady State Min Iterations = 1
! Output File = "$namerun".result"
Post File = "$namerun".vtu"
max output level = 3
End
!---------------------------------------------------
!---------------- BODIES ---------------------------
!---------------------------------------------------
! This body is located at the ice/bed interface and will be used to solve
! the sheet equation
Body 1
Name= "sheet"
Equation = 1
Material = 1
Body Force = 1
Initial Condition = 1
End
!---------------------------------------------------
!---------------- INITIAL CONDITIONS ---------------
!---------------------------------------------------
!! for the sheet
Initial Condition 1
Sheet Thickness = Real $Hs
Zs = Variable Coordinate 1
Real MATC "H(tx)"
Zb = Real 0.0
Effective Pressure = Variable Coordinate 1
Real MATC "rhoi*gravity*H(tx)"
End
!---------------------------------------------------
!---------------- BODY FORCES ----------------------
!---------------------------------------------------
! source 1 cm/day
Body Force 1
Hydraulic Potential Volume Source = REal $Source
End
!---------------------------------------------------
!---------------- MATERIALS ------------------------
!---------------------------------------------------
!! ice material properties in MPa - m - a system
Material 1
Density = Real $rhoi
Glen Exponent = Real $ng
! For the sheet
Sheet Conductivity = Real $Ks
Sheet flow exponent alpha = Real $alphas
Sheet flow exponent beta = Real $betas
Englacial Void Ratio = Real $ev
Sliding Velocity = Real $ub
Bedrock Bump Length = Real $lr
Bedrock Bump High = Real $hr
Sheet Closure Coefficient = Real $Ar
! For the Channels
Channel Conductivity = Real $Kc
Channel flow exponent alpha = Real $alphac
Channel flow exponent beta = Real $betac
Channel Closure Coefficient = Real $Ac
Sheet Width Over Channel = Real $lc
Pressure Melting Coefficient = Real $Ct
Water Heat Capacity = Real $Cw
! For both
Ice Normal Stress = Variable Coordinate 1
Real MATC "rhoi*gravity*H(tx)"
End
!---------------------------------------------------
!---------------- SOLVERS --------------------------
!---------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
Equation = "GlaDS Coupled sheet"
Procedure = "ElmerIceSolvers" "GlaDSCoupledSolver"
Variable = -dofs 1 "Hydraulic Potential"
Activate Channels = Logical True
Activate Melt from Channels = Logical True
Neglect Sheet Thickness in Potential = Logical True
! choices are EXPLICIT, CRANK-NICOLSON, IMPLICIT
Channels Integration method = String "Crank-Nicolson"
Sheet Integration method = String "Crank-Nicolson"
Exported Variable 1 = -dofs 1 "Vclose"
Exported Variable 2 = -dofs 1 "Wopen"
Exported Variable 3 = -dofs 1 "Normal Stress"
Exported Variable 4 = -dofs 1 "Water Pressure"
Exported Variable 5 = -dofs 1 "Effective Pressure"
Exported Variable 6 = -dofs 2 "Sheet Discharge"
Exported Variable 7 = -dofs 1 "Sheet Storage"
Exported Variable 8 = -dofs 1 "Zs"
Exported Variable 9 = -dofs 1 "Zb"
Exported Variable 10 = -dofs 1 "Flux from Moulins"
Linear System Solver = Direct
Linear System Direct Method = umfpack
Nonlinear System Max Iterations = 10
Nonlinear System Convergence Tolerance = 1.0e-6
Nonlinear System Relaxation Factor = 1.00
Coupled Max Iterations = Integer 10
Coupled Convergence Tolerance = Real 1.0e-3
Steady State Convergence Tolerance = 1.0e-03
End
Solver 2
! Just here to declare the variable Sheet Thickness
Equation = "GlaDS Thickness sheet"
Procedure = "ElmerIceSolvers" "GlaDSsheetThickDummy"
Variable = -dofs 1 "Sheet Thickness"
End
Solver 3
! Just here for output and declare the Channel Area variable
! It is executed simultaneously to saving
Exec Solver = After Saving
Equation = "GlaDS Channel OutPut"
Procedure = "ElmerIceSolvers" "GlaDSchannelOut"
Variable = -dofs 1 "Channel Area"
! Define that the variable is define on the edges only
Element = "n:0 e:1"
Exported Variable 1 = -dofs 1 "Channel Flux"
VTU OutPutFile = Logical False
VTU BinaryFile = Logical False
Channels OutPut Directory Name = String "results"
Channels OutPut File Name = String "$namerun"_channels"
End
!---------------------------------------------------
!---------------- EQUATIONS ------------------------
!---------------------------------------------------
! Equation for the ice
Equation 1
Active Solvers (3) = 1 2 3
End
!---------------------------------------------------
!---------------- BOUNDARY CONDITIONS --------------
!---------------------------------------------------
Boundary Condition 1
Target Boundaries(2) = 1 3
No Channel BC = Logical True
End
Boundary Condition 2
Name = "Sym"
Target Boundaries(1) = 2
No Channel BC = Logical True
End
Boundary Condition 3
Name = "front"
Target Boundaries(1) = 4
Hydraulic Potential = Real 0.0
No Channel BC = Logical True
End
! This is the boundary condition for the moulins, case B5
Boundary Condition 4
Name = "moulins"
Target Boundaries(100) = 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
Moulin Storage = Logical True
Moulin Area = Real $Am
Moulin Flux = Real $0.9*yearinsec
End
Solver 1 :: Reference Norm = Real 8.28199040E+00
Solver 1 :: Reference Norm Tolerance = Real 1E-04
