#
echo on
!Conductivity
$ function conductivity(T) {_conductivity=9.828*exp(-5.7e-03*T)}
!Capacity
$ function capacity(T) {_capacity=146.3+(7.253*T)}
!Pressure melting point
$ function pressuremeltingpoint(PIN) {\
P = PIN;\
if (P<0.0) P=0.0;\
beta=9.8E-08*1.0E06;\
_pressuremeltingpoint=273.15-(beta*P);\
}
check keywords warn
$namerun = "teterousse"
$yearindays = 365.25
$yearinsec = 365.25*24*60*60
$Rho = 9.03760022565934E-019
$RhoFW = 1.00417780285104E-018
$RhoWS = 1.03430313693657E-018
$g = 9.76960357522560E015 !equivalent to 9.81*yearinsec
$MPainPa = 1.0E6
$ub = 4000.0/yearinsec !Need to link this up to the NS solver
$A1 = 2.89165E-13*yearinsec*MPainPa^3
$A2 = 2.42736E-2*yearinsec*MPainPa^3
$ng = 3.0
$Aglen = 2.5E-25*yearinsec*MPainPa^3 !Need to link this up to actual viscosity
$eta = (2.0*Aglen)^(-1.0/ng)
!For the Sheet
$Ar = Aglen
$alphas = 1.25
$betas = 1.5
$lr = 5.0
$hr = 0.1
$Ks = 0.005*yearinsec*(1.0/MPainPa)^(1.0-betas)
$Hs = 0.01
$ev = 0.01
!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 = 0
!Source Term
$Source = 2.0e-10*yearinsec
!Regularisation Parameter
!$Lambda = 1.0E10
!For Glen
$GlenN = 3.0
$SL = 0.0
$Iter = 100
$dt = 0.1/yearindays
$OutPut = 1
$OutPutdt = 1.0/yearindays
$FinalTime = Iter*dt
Header
Mesh DB "." "teterousse"
!Include path "include"
!Results Directory "./Results"
End
Constants
Water Density = Real $RhoWS
Latent Heat = Real $Lw
Gravity Norm = Real $g
Ice Density = Real $Rho
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 "ZbDEM"
End
Simulation
Coordinate System = Cartesian 3D
Simulation Type = transient
Timestepping Method = "bdf"
BDF Order = 1
Timestep Intervals = $Iter
Output Intervals = $OutPut
Timestep Sizes = $dt
!Restart File = "$restart".result" !"
!Restart Position = 0
!Restart Time = 0.0
!Restart Before Initial Conditions = Logical True
!Restart Variable 1 = BedDEM
!Restart Variable 2 = SurfDem
!Restart Variable 3 = vx
!Restart Variable 4 = vy
!Restart Variable 5 = vmag
!Restart Variable 6 = Coordinate 1
!Restart Variable 7 = Coordinate 2
!Restart Variable 8 = Coordinate 3
Steady State Max Iterations = 10
Steady State Min Iterations = 1
Set Dirichlet BCs By BC Numbering = Logical True
max output level = 100
Output File = "$namerun".result"
Post File = "$namerun".vtu"
Output Coordinates = Logical True
Extruded Mesh Levels = Integer 15
End
!!!----------------------- BODIES -----------------------!!!
Body 1
name = ice
Equation = 1
Material = 1
Body Force = 1
Initial Condition = 1
End
Body 2
name = sheet
Equation = 2
Material = 2
Body Force = 2
Initial Condition = 2
End
!!!----------------------- MATERIALS -----------------------!!!
Material 1
Density = Real $Rho
Viscosity Model = String "Glen"
!Have to set viscosity to dummy value to remove warning output
Viscosity = Real $1.0E13*yearinsec*1.0E-6
!n in the Flow Law
Glen Exponent = Real $ng
!Rate Factors (from Paterson, in MPa^-3a^-1)
Rate Factor 1 = Real 1.258E13
Rate Factor 2 = Real 6.046E28
!These ones in SI units instead - fine, as long as gas constant also is.
Activation Energy 1 = Real 60E3
Activation Energy 2 = Real 139E3
Glen Enhancement Factor = Real 1.0
!Need a variable to evaluate the Arrhenius Law - ideally, should be temperature !relative to PMP. This will plug in right value from the TemperateIceSolver.
! Temperature Field Variable = String "Temp Homologous"
!Temperature to switch between regimes in flow law
Limit Temperature = Real -10.0
!Can just specify a constant temperature instead to simplify things.
Constant Temperature = Real -4.0
Latent Heat = Real $3.35E05*(31556926.0)^(2.0)
Critical Shear Rate = Real 1.0E-10
!Need heat capacity - here as MATC function of temperature itself
Temp Heat Capacity = Variable Temp
Real MATC "capacity(tx)*(3.1556926.0)^2"
!Same for heat conductivity
Temp Heat Conductivity = Variable Temp
Real MATC "conductivity(tx)*(3.1556926.0)*1.0E-06"
!Upper temperature limit as PMP
Temp Upper Limit = Variable Pressure
Real MATC "pressuremeltingpoint(tx)"
!Lower limit - may as well say 0K.
Temp Lower Limit = Real 0.0
Sea Level = Real $SL
Cauchy = Logical True
Youngs Modulus = Real 1.0
Poisson Ratio = Real 0.3
!Minimum value of Free Surface Variable (i.e. how thin the ice can get)
Min Zs = Real 100.0
End
Material 2
Density = Real $Rho
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 !Realistic Store value?
Sliding Velocity = Variable Velocity 1, Velocity 2
Real MATC "sqrt((tx(0)^2)+(tx(1)^2))"!Equals vmag!Equals Flow Solution
Bedrock Bump Length = Real $lr !What's a realistic Store value?
Bedrock Bump High = Real $hr !Realistic Store value
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 Stress 3
Real MATC "abs(tx)"
!Ice Normal Stress = Variable ZsDEM, ZbDEM
! Real MATC "abs(Rho*gravity*(tx(0)-tx(1)))"
Cauchy = Logical False
End
!!!----------------------- BODY FORCES -----------------------!!!
Body Force 1
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real 0.0
Flow BodyForce 3 = Real $-g
End
!For the sheet
Body Force 2
Hydraulic Potential Volume Source = Real $Source
End
!!!-------------------------SOLVERS------------------------!!!
Solver 1
Exec Solver = "Before Simulation"
Equation = "Read DEMs"
Procedure = "ElmerIceSolvers" "Grid2DInterpolator"
! Bedrock DEM
Variable 1 = String "ZbDEM"
Variable 1 data file = File "BedDEM.xyz"
Variable 1 x0 = Real 947700.0d0
Variable 1 y0 = Real 2104850.0d0
Variable 1 lx = Real 600.0
Variable 1 ly = Real 350.0
Variable 1 Nx = Integer 301
Variable 1 Ny = Integer 176
Variable 1 Invert = Logical False
Variable 1 Fill = Logical False
Variable 1 Position Tol = Real 1.0e-1
Variable 1 No Data = Real -9999
Variable 1 No Data Tol = Real 1.0
! Surface DEM
Variable 2 = String "ZsDEM"
Variable 2 data file = File "SurfDEM.xyz"
Variable 2 x0 = Real 947700.0d0
Variable 2 y0 = Real 2104850.0d0
Variable 2 lx = Real 800.0
Variable 2 ly = Real 350.0
Variable 2 Nx = Integer 268
Variable 2 Ny = Integer 118
Variable 2 Invert = Logical False
Variable 2 Fill = Logical False
Variable 2 Position Tol = Real 1.0e-1
Variable 2 No Data = Real -9999.0
Variable 2 No Data Tol = Real 1.0
End
Solver 2
! Exec Solver = "Never"
Equation = "MapCoordinate"
Procedure = "StructuredMeshMapper" "StructuredMeshMapper"
Active Coordinate = Integer 3
!Mesh Velocity Variable = String "dSdt"
!Mesh Update Variable = String "dS"
Mesh Velocity First Zero = Logical True
Top Surface Variable Name = String "ZsDEM"
Bottom Surface Variable Name = String "ZbDEM"
Displacement Mode = Logical False
Correct Surface = Logical True
Minimum Height = Real 1.0
End
Solver 3
Equation = "Depth"
Exec Solver = "Before Timestep"
Procedure = File "ElmerIceSolvers" "FlowDepthSolver"
Variable = String "Depth"
Variable DOFs = 1
!These variable need to be exported from somewhere to be used in first solver
Exported Variable 1 = "ZbDEM"
Exported Variable 2 = "ZsDEM"
Gradient = Real -1.0E0
Calc Free Surface = Logical False
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End
Solver 4
Equation = "Elevation"
Exec Solver = "Before Timestep"
Procedure = File "ElmerIceSolvers" "FlowDepthSolver"
Variable = String "Elevation"
Variable DOFs = 1
Gradient = Real 1.0E0
Calc Free Surface = Logical False
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End
Solver 5
Equation = "Navier-Stokes"
Exec Solver = "Before Timestep"
Stabilize = Logical True
Flow Model = Stokes
Variable = Flow Solution[Velocity:3 Pressure:1]
Solver Timing = Logical True
!Mandatory to save bulk stiffness matrix
Calculate Loads = Logical True
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStabl"
BiCGStabl Polynomial Degree = 4
Linear System Max Iterations = 1000
Linear System Convergence Tolerance = Real 1.0E-6
Linear System Abort Not Converged = False
Linear System Preconditioning = ILU0
Linear System Residual Output = 0
Nonlinear System Max Iterations = 100
Nonlinear System Convergence Tolerance = 1.0E-5
Nonlinear System Newton After Iterations = 10
Nonlinear System Newton After Tolerance = 5.0E-3
Nonlinear System Newton Max Tolerance = Real 1.0E-2
Nonlinear System Newton Max Iterations = Integer 20
Nonlinear System Reset Newton = Logical True
Steady State Convergence Tolerance = Real 1.0E-12
End
Solver 6
Equation = "Distance"
Exec Solver = "Before Timestep"
Procedure = "DistanceSolve" "DistanceSolver1"
Variable = "Distance"
Nonlinear System Max Iterations = 200
Nonlinear System Convergence Tolerance = 1.0E-5
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Direct Method = UMFPack
Steady State Convergence Tolerance = 1.0E-4
! Exported Variable 1 = "vx"
! Exported Variable 2 = "vy"
! Exported Variable 3 = "SurfDEM"
! Exported Variable 4 = "BedDEM"
! Exported Variable 1 = "Zs Top"
Exported Variable 1 = "Reference Zs Top"
Exported Variable 2 = "Zs"
Exported Variable 3 = "Zb"
End
Solver 7
Equation = "GlaDS Coupled sheet"
!Exec Solver = "Never"
Procedure = "../../SOFTWARE/GlaDSCoupledSolver" "GlaDSCoupledSolver"
Variable = -dofs 1 "Hydraulic Potential"
Activate Channels = Logical True
Activate Melt from Channels = Logical True
! arr added
Neglect Sheet Thickness in Potential = Logical True
Channels Integration method = String "Crank-Nicolson"
Sheet Integration method = String "Crank-Nicolson"
!arr added
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 1 "Sheet Discharge"
Exported Variable 7 = -dofs 1 "Sheet Storage"
Exported Variable 8 = -dofs 1 "Zs"
Exported Variable 9 = -dofs 1 "Zb"
Linear System Solver = Direct
Linear System Direct Method = umfpack
Nonlinear System Max Iterations = 30
Nonlinear System Convergence Tolerance = 1.0E-6
Nonlinear System Relaxation Factor = 1.00
Coupled Max Iterations = Integer 40
Coupled Convergence Tolerance = Real 1.0E-3
Steady State Convergence Tolerance = 1.0E-3
End
Solver 8
!Purely for declaring sheet thickness variable
Equation = "GlaDS Thickness sheet"
!Exec Solver = "Never"
Procedure = "../../SOFTWARE/GlaDSCoupledSolver" "GlaDSsheetThickDummy"
Variable = -dofs 1 "Sheet Thickness"
End
Solver 9
!Output and declare Channel Area variable
Equation = "GlaDS Channel OutPut"
!Exec Solver = "Never"
Procedure = "../../SOFTWARE/GlaDSchannelSolver" "GlaDSchannelOut"
Variable = -dofs 1 "Channel Area"
!Define variable on edges only
Element = "n:0 e:1"
Exported Variable 1 = -dofs 1 "Channel Flux"
VTK OutPutFile = Logical True
ASCII OutPutFile = Logical True
Channels OutPut Directory Name = String "./channel_output"
Channels OutPut File Name = String "$namerun"_channels"
OutPut Time = Real $OutPutdt
Final Time = Real $FinalTime
Steady State Convergence Tolerance = 1.0E-3
End
Solver 10
Equation = String "StressSolver"
Exec Solver = "Before Timestep"
Procedure = File "ElmerIceSolvers" "ComputeDevStress"
Variable = -nooutput "Sij" !This is just a dummy, so no output needed
Variable DOFs = 1
Stress Variable Name = String "Stress"
!The name of the variable containing the flow solution (U, V, W, Pressure)
Flow Solver Name = String "Flow Solution"
Exported Variable 1 = "Stress" ![Sxx, Syy, Szz, Sxy] in 2D
![Sxx, Syy, Szz, Sxy, Syz, Szx] in 3D
Exported Variable 1 DOFs = 6 !4 in 2D
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0E-1
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
End
Solver 11
Equation = "ResultOutput"
Exec Solver = "After Saving"
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = "$namerun"_" !"
Vtu Format = Logical True
Binary Output = True
Single Precision = True
Save Geometry IDs = True
Vector Field 1 = "Velocity"
Vector Field 2 = "Stress"
Scalar Field 1 = "Depth"
Scalar Field 2 = "Elevation"
Scalar Field 3 = "Distance"
! Scalar Field 4 = "Temp"
Scalar Field 4 = "Zs"
Scalar Field 5 = "Reference Zs Top"
! Scalar Field 7 = "smb"
! Scalar Field 7 = "vx"
! Scalar Field 8 = "vy"
! Scalar Field 9 = "vmag"
Scalar Field 6 = "Zb"
Scalar Field 7 = "Pressure"
Scalar Field 8 = "Sheet Thickness"
Scalar Field 9 = "Vclose"
Scalar Field 10 = "Wopen"
Scalar Field 11 = "Effective Pressure"
Scalar Field 12 = "Sheet Discharge"
Scalar Field 13 = "Sheet Storage"
Scalar Field 14 = "Hydraulic Potential"
Scalar Field 15 = "Normal Stress"
Scalar Field 16 = "ZbDEM"
Scalar Field 17 = "ZsDEM"
End
!!!----------------------- EQUATION -----------------------!!!
Equation 1
Active Solvers (8) = 1 2 3 4 5 6 10 11
End
Equation 2
Active Solvers (3) = 7 8 9
End
!!!----------------------- Initial Conditions -----------------------!!!
! IC ICE
Initial Condition 1
Zs = Equals ZsDEM
Reference Zs Top = Equals ZsDEM
Zb = Equals ZbDEM
Mesh Update 1 = Real 0.0
Mesh Update 2 = Real 0.0
Mesh Update 3 = Real 0.0
! Pressure = Real 0.0
! Velocity 1 = Real 0.0
! Velocity 2 = Real 0.0
! Velocity 3 = Real 0.0
End
!IC water sheet
Initial Condition 2
Sheet Thickness = Real $Hs
End
!!!!!!!!!!!!!!!!!!!!! BOUNDARY CONDTIONS !!!!!!!!!!!!!!!!!!!!!
!Walls = 2, Terminus = 1, Upper limit = 8
Boundary Condition 1
name = "Inflow"
Target Boundaries(1) = 8
! Body ID = 2
No Channel BC = Logical True
Inflow Mask = Logical True
Flow Force BC = Logical True
Velocity 1 = Real 1.0
Velocity 2 = Real 1.0
Zs = Equals Reference Zs Top
End
Boundary Condition 2
name = "Terminus"
Target Boundaries(1) = 1
! Body ID = 2
Distance = Real 0.0
Hydraulic Potential = Real 0.0
No Channel BC = Logical True
End
Boundary Condition 3
name = "Walls"
Target Boundaries(1) = 2
! Body ID = 2
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
Velocity 2 Condition = Variable Elevation
Real MATC "-(tx-1.0)"
Velocity 3 = Real 0.0
Velocity 3 Condition = Variable Elevation
Real MATC "-(tx-1.0)"
Slip Coefficient 2 = Real 1.0E-2
Slip Coefficient 3 = Real 1.0E-2
No Channel BC = Logical True
Normal-Tangential Velocity = Logical True
Mass Consistent Normals = Logical True
ComputeNormal = Logical True
End
Boundary Condition 4
name = "bed"
Body ID = 2
Bottom Surface Mask = Logical True
Elevation = Real 0.0
Save Line = Logical True
Flow Force BC = Logical True
Slip Coefficient 2 = Real 1.0E-2
Slip Coefficient 3 = Real 1.0E-2
Normal-Tangential Velocity = Logical True
Mass Consistent Normals = Logical True
ComputeNormal = Logical True
Velocity 1 = Real 0.0
End
Boundary Condition 5
name = "Surface"
Top Surface Mask = Logical True
Depth = Real 0.0
Save Line = Logical True
Flow Force BC = Logical True
Temp = Real -4.0 !Variable Coordinate 1
Mesh Update 1 = Real 0.0
Mesh Update 2 = Real 0.0
Mesh Update 3 = Variable Zs
Real Procedure "ElmerIceUSF" "ZsTopMzsIni"
End