#
#
#
! LAS: 2D ice shelf geometry 50 km by 1 km.
! define initial values
$yearinsec = 365.25*24*60*60
$rhoi = 910.0/(1.0e6*yearinsec^2)
$rhow = 1000.0/(1.0e6*yearinsec^2)
! g= 9.81*yearinsec^2 m a-2
$g = 9.7692e15 ! this should be 9.7692e15, I increased it ten times
$Ts = 248.15 !-10C
$QG=1.5779
$PI = 3.14
$H = 1000 ! total height of domain
$L = 5e4 ! total length of domain
$DivPosition = 0
$IceCorePosition = 0.0
$PeakM = 5.6 ! to balance a surface mass balance of +0.7 m yr^-1
$BackM = 0.0 ! ice shelf melt on flanks
$Thinning_a = 0 ! mean thinning due to imbalance in surface and basal mass balances
$Thinning_f = 0 ! mean thinning due to far-field ice dynamics - imbalance between flow in/out at margins.
! variable for the normal distribution of melt
$sigma = 2500 ! seems to yield a channel order few kilometers
! accumulation (defined here for choosing the timesteps)
$a0 = 0.7
$SeaLevel = H*10/11 ! in [m] relative to the modle domain
$n=3
$Factor=1
! conductivity and heat capacity function of temperature
$ function conductivity(T) { _conductivity=9.828*exp(-5.7E-03*T)}
$ function capacity(T) { _capacity=146.3+(7.253*T)}
!! Glen's flow law----- with the Arrhenius factor a function of temperature
! units are Pa^-3 yr^-1------------agrees with Paterson
$ function A(T) { _A = 1.0e-15*(0.2071*exp(0.5978*(T-273.15))+0.09833*exp(0.14747*(T-273.15)))}
!! pressuremeltingpointinwater (in SI units)
$ function pressuremeltingpointinwater(D) {\
P = 1000*D*9.81;\
if (P<0.0) P=0.0;\
beta=9.8E-08;\
_pressuremeltingpointinwater=273.15-(beta*P);\
}
!! pressuremeltingpointinice (in SI units)
$ function pressuremeltingpointinice(D) {\
P = 910*D*9.81;\
if (P<0.0) P=0.0;\
beta=9.8E-08;\
_pressuremeltingpointinice=273.15-(beta*P);\
}
Header
Mesh DB "." "mesh"
Include Path ""
Results Directory ""
End
Constants
Stefan Boltzmann = Real 0.0
Gas Constant = Real 8.314 !Joule/mol x K
! For the Buoyancy User function
Buoyancy Use Basal Melt = Logical True
Bottom Surface Name = String "hb"
Water Density = Real $rhow
End
Simulation
Coordinate System = "Cartesian 2D"
Simulation Type = Transient
Timestepping Method = "bdf"
BDF Order = 1
Timestep Intervals = 1000 ! run for 1 characteristic time
Output Intervals = 10
Timestep Sizes = Real $ (H/a0)*(1.0/1000.0) !
Steady State Min Iterations = 1
Steady State Max Iterations = 1
Set Dirichlet BCs By BC Numbering = Logical True
Output File = "Uniform1.dat"
!Restart File = "SteadyDivideInCentre_1.dat"
max output level = 3
End
Body 1
name = Ice
Equation = 1
Material = 1
Body force = 1
Initial Condition = 1
End
! The Free surface
Body 2
name = FreeSurf
Equation = 2
Body Force = 2
Material = 1
Initial Condition = 2
End
! The Free Base
Body 3
name = FreeBase
Equation = 3
Body Force = 3
Material = 1
Initial Condition = 3
End
Body Force 1
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real $ -g
Friction Heat = Logical True
End
Body Force 2
hs Accumulation Flux 1 = Real 0.0e0
hs Accumulation Flux 2 = Real $ a0
hs Accumulation Flux 3 = Real 0.0e0
End
Body Force 3
!! melting/accretion under ice/shelf
!! positive for melting
!! negative for accretion
hb Accumulation = Variable Coordinate 1
Real MATC "BackM + PeakM*exp(-1*(pow(tx,2))/(2*pow(sigma,2)))"
End
Initial Condition 1
Depth = Variable Coordinate 2
Real MATC "H-tx"
Height = Variable Coordinate 2
Real MATC "tx"
Age = Variable Depth, Height
Real MATC "tx(0)/a0"
! linear change between Ts at the surface and the pressure melting point at the base given by pressuremeltingpointinice(tx(0)+tx(1)).
Temperature = Variable Depth, Height
Real MATC "Ts*(1+tx(0)/(tx(0)+tx(1))*(pressuremeltingpointinice(tx(0)+tx(1))/Ts - 1))
End
Initial Condition 2
hs = Equals Coordinate 2
hsREF = Equals Coordinate 2
End
Initial Condition 3
hb = Equals Coordinate 2
hbREF = Equals Coordinate 2
End
Material 1
Density = Real $ rhoi
Viscosity Model = String "power law"
Viscosity Exponent = Real $ 1.0/n
Critical Shear Rate = Real $ 1.0e-3/31556926.0
Viscosity = Variable Temperature ! 10e-6/(2A)^(1/n) MPa a^1/3
Real MATC "1e-6*pow(2.0*A(tx)*Factor,-1.0/n)"
Heat Capacity = Variable Temperature
Real MATC "capacity(tx)*(31556926.0)^(2.0)"
Heat Conductivity = Variable Temperature
Real MATC "conductivity(tx)*(31556926.0)*1.0E-06"
! used by the Buoyancy User Function
Sea level = Real $SeaLevel
End
!!!!! SOLVERS
Solver 1
Equation = "Mesh Update"
Exec Solver = "Before Timestep"
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0e-6
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
End
! Flow-depth on the unstructured FEM-mesh
! not really needed, but nice for post-processing
! Just un-comment Exec Solver = Never, in order
! to skip it
!-------------------------------------------------
Solver 2
Equation = "Flowdepth"
Exec Solver = "Before TimeStep"
Procedure = File "src/Flowdepth" "FlowDepthSolver"
Variable = String "Depth"
Variable DOFs = 1
Linear System Solver = "Iterative"
! Linear System Solver = "Direct"
Linear System Direct Method = "MUMPS" ! MUMPS only works in parallel, o if in serial and using direct method, comment this and it still works.
Linear System Iterative Method = "GCR"
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
! this sets the direction
! -1 is negative z-direction (upside down)
! +1 is positive (bottom up)
Gradient = Real -1.0E00
! switch that to True, if you
! want to have free surface gradients
! to be computed
!------------------------------------
Calc Free Surface = Logical True
! NEW FlowDepth adds the variable automatically
Freesurf Name = String "FreeSurf"
! THE REST OF THIS SECTION IS NOT NEEDED WITH THE NEW FlowDepth Solver
! will contain the variable of the corresponding free surface position
!---------------------------------------------------------------------
Exported Variable 1 = String "FreeSurf"
Exported Variable 1 DOFs = 1
! next two will contain values of the free surface gradient
! nice to be used in post-processing for evaluation of the
! surface mass balance in diagnostic runs
!----------------------------------------------------------
Exported Variable 2 = String "FreeSurfGrad1"
Exported Variable 2 DOFs = 1
Exported Variable 3 = String "FreeSurfGrad2"
Exported Variable 3 DOFs = 1
End
! This solves the height above the bedrock
!------------------------------------------
Solver 3
Equation = "Flowheight" ! mind different name
Exec Solver = "Before TimeStep"
Procedure = File "src/Flowdepth" "FlowDepthSolver"
Variable = String "Height" ! mind different name for variable
Variable DOFs = 1
! Linear System Solver = "Direct"
! Linear System Direct Method = "Umfpack"
Linear System Solver = "Iterative"
Linear System Max Iterations = 500
Linear System Iterative Method = "GCR"
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
Gradient = Real 1.0E00 ! this time positive
Calc Free Surface = Logical False
End
Solver 4
Equation = "Heat Equation"
Stabilization Method = Stabilized
! Linear System Solver = "Direct"
Linear System Direct Method = "umfpack"
Linear System Solver = "Iterative"
Linear System Iterative Method = "GCR"
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-6
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU2"
Linear System Residual Output = 0
Nonlinear System Max Iterations = 1
Nonlinear System Convergence Tolerance = 1.0e-4
Nonlinear System Newton After Iterations = 1000
Nonlinear System Newton After Tolerance = 1.0e-5
Nonlinear System Relaxation Factor = 1
Steady State Convergence Tolerance = 1.0e-3
End
! the Navier-Stokes Solver, solves ice flow dynamics
!---------------------------------------------------
Solver 5
Equation = "Navier-Stokes"
Flow model = String "Stokes"
Linear System Solver = "Iterative"
Linear System Iterative Method = "GCR"
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU2"
Linear System Residual Output = 0
Nonlinear System Max Iterations = 1
Nonlinear System Convergence Tolerance = 1.0E-04
Nonlinear System Newton After Iterations = 1000
Nonlinear System Newton After Tolerance = 1.0E-05
Nonlinear System Relaxation Factor = 1
Steady State Convergence Tolerance = 1.0E-03
! Stabilization Method = [Stabilized,P2/P1,Bubbles]
Stabilization Method = Stabilized
End
Solver 6
Equation = String "Free Surface Evolution"
Exec Solver = "After TimeStep"
Variable = "hs"
Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
Linear System Solver = Iterative
Linear System Iterative Method = "GCR"
Linear System Max Iterations = 500
Linear System Preconditioning = "ILU1"
Linear System Convergence Tolerance = 1.0e-06
Optimize Bandwidth = Logical False
Nonlinear System Max Iterations = 1
Nonlinear System Min Iterations = 1
Nonlinear System Convergence Tolerance = 1.0e-05
Steady State Convergence Tolerance = 1.0e-3
Apply Dirichlet = Logical False
Flow Solution Name = String "Flow Solution"
Stabilization Method = Stabilize
Optimize Bandwidth = Logical False
Exported Variable 1 = hs Residual
Exported Variable 1 DOFS = 1
Exported Variable 2 = hsREF
Exported Variable 2 DOFS = 1
End
Solver 7
Equation = String "Free Base Evolution"
Exec Solver = "After TimeStep"
Variable = "hb"
Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
Linear System Solver = Iterative
Linear System Iterative Method = "GCR"
Linear System Max Iterations = 500
Linear System Preconditioning = "ILU1"
Linear System Convergence Tolerance = 1.0e-06
Optimize Bandwidth = Logical False
Nonlinear System Max Iterations = 1
Nonlinear System Min Iterations = 1
Nonlinear System Convergence Tolerance = 1.0e-05
Steady State Convergence Tolerance = 1.0e-3
Apply Dirichlet = Logical False
Flow Solution Name = String "Flow Solution"
Stabilization Method = Stabilize
Optimize Bandwidth = Logical False
Exported Variable 1 = hb Residual
Exported Variable 1 DOFS = 1
Exported Variable 2 = hbREF
Exported Variable 2 DOFS = 1
End
! Age Solver
!
Solver 8
Equation = "Age Equation"
Exec Solver = "After TimeStep"
Variable = String "Age"
Variable DOFs = 1
Optimize Bandwidth = Logical False
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = Diagonal
Linear System Preconditioning = NoNe
Linear System Convergence Tolerance = Real 1.0e-6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
Procedure = "./src/AgeSolver" "AgeSolver"
End
Solver 9
Equation = "Velocity gradient tensor"
Exec Solver = "After TimeStep"
Procedure = "./src/DeformationSolver" "DSolver"
Variable = String "VG[D:6 W:3]"!dxx dyy dzz dxy dxz dyz wxy wxz wyz in 3D
Variable DOFs = 9
Optimize Bandwidth = Logical False
Linear System Solver = "Iterative"
Linear System Iterative Method = "Diagonal"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0e-6
Linear System Abort Not Converged = False
Linear System Preconditioning = "None"
Linear System Residual Output = 0
End
Solver 10
Equation = "Surface Data"
Exec Solver = "After Saving"
Procedure = "SaveData" "SaveLine"
Filename = "Surface.dat"
Polyline Coordinates(2,2) = Real $ IceCorePosition -250 IceCorePosition 1250
End
Solver 11
Equation = "VTU Output"
Exec Solver = "After Saving"
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = "Output_PeakM1"
Output Format = String "vtu"
Vector Field 1 = String "Velocity"
! Vector Field 2 = String "VG"
Scalar Field 1 = String "Depth"
Scalar Field 2 = String "Height"
Scalar Field 3 = String "Pressure"
Scalar Field 4 = String "Temperature"
Scalar Field 5 = String "Age"
Scalar Field 6 = String "D 2"
Scalar Field 7 = String "D 4"
Scalar Field 8 = String "D 1"
! Scalar Field 9 = String "D 7"
End
!!!!! EQUATION
Equation 1
Active Solvers (9) = 1 2 3 4 5 8 9 10 11
Flow Solution Name = String "Flow Solution"
Convection = String "Computed"
End
Equation 2
Active Solvers(1) = 6
Flow Solution Name = String "Flow Solution"
Convection = String "Computed"
End
Equation 3
Active Solvers(1) = 7
Flow Solution Name = String "Flow Solution"
Convection = String "Computed"
End
!!!!! BOUNDARY CONDITIONS
Boundary Condition 1
name = Base
Target Boundaries = 1
Body Id = 3
Height = Real 0.0
Normal-Tangential Velocity = Logical True
Flow Force BC = Logical True
External Pressure = Variable Coordinate 2
Real Procedure "./src/Buoyancy" "SeaPressure"
Slip Coefficient 1 = Variable Coordinate 2
Real Procedure "./src/Buoyancy" "SeaSpring"
Compute Sea Pressure = Logical True
Compute Sea Spring = Logical True
Temperature = Variable Coordinate 2
Real MATC "pressuremeltingpointinwater(SeaLevel -tx)"
Mesh Update 1 = Real 0.0
Mesh Update 2 = Variable hb, hbREF
Real MATC "tx(0)-tx(1)"
Save Line = Logical True
End
Boundary Condition 2
name = Surface
Target Boundaries = 3
Body Id = 2
Depth = Real 0.0
Temperature = Real $ Ts
Age = Real 0.0
Mesh Update 1 = Real 0.0
Mesh Update 2 = Variable hs, hsREF
Real MATC "tx(0)-tx(1)"
Save Line = Logical True
End
Boundary Condition 3
name = Left
Target Boundaries = 2
Velocity 1 = 0
External Pressure = Variable Depth
Real MATC "-rhoi*g*tx"
Mesh Update 1 = Real 0.0
Mesh Update 2 = Variable hb, hbREF
Real MATC "tx(0)-tx(1)"
Save Line = Logical False
End
Boundary Condition 4
name = Right
Target Boundaries = 4
Velocity 1 = 0
External Pressure = Variable Depth
Real MATC "-rhoi*g*tx"
Mesh Update 1 = Real 0.0
Mesh Update 2 = Variable hb, hbREF
Real MATC "tx(0)-tx(1)"
Save Line = Logical False
End