!echo on
Header
CHECK KEYWORDS Warn
Mesh DB "." "EISMINT"
Include Path ""
Results Directory ""
End
!!!!!!!!!!!!!!!!!!!!!!
!!! User Functions !!!
!!!!!!!!!!!!!!!!!!!!!!
!! conductivity
$ function conductivity(T) { _conductivity=9.828*exp(-5.7E-03*T)}
!! capacity
$ function capacity(T) { _capacity=146.3+(7.253*T)}
!!accumulation rate
$ function accum(X) {\
lapserate = (.49/525.0);\
ela = 1500.0;\
asl = -ela*lapserate;\
if (X(0) < 400000)\
{_accum = .5;}\
else\
{ _accum = 1.0e-5*(450000 - X(0));}\
}
!!! The temperature initialization condition
!!! This keeps the solver from applying the lapse rate BC to the first timestep
$ function surftemp(X) {\
if (X(2) < 2)\
{_surftemp = 260.0;}\
else\
{_surftemp = 270.0 - 10.0*X(1)/1000.0;}\
}
!! Glen's flow law (using power law)
!-----------------
$ function glen(Th) {\
EF = 1.0;\
AF = getArrheniusFactor(Th);\
_glen = (2.0*EF*AF)^(-1.0/3.0);\
}
!! Arrhenius factor needed by glen Th is the HOMOLOGOUS Temperature
!! (in SI units)
!---------------------------------
$ function getArrheniusFactor(Th){ \
if (Th<-10) {_getArrheniusFactor=3.985E-13 * exp( -60.0E03/(8.314 * (273.15 + Th)));}\
else {\
if (Th>0) _getArrheniusFactor=1.916E03 * exp( -139.0E03/(8.314 * (273.15)));\
else _getArrheniusFactor=1.916E03 * exp( -139.0E03/(8.314 * (273.15 + Th)));}\
}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Simulation Params !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Max Output Level = 4
Coordinate System = "Cartesian 2D"
Coordinate Mapping(3) = 1 2 3
Simulation Type = "Transient"
Steady State Max Iterations = 1
Timestepping Method = "BDF"
BDF Order = 2
Timestep Sizes = 1.0
Timestep Intervals = 15000
Output Intervals = 100
!Output File = "Eismint_freesurfacediroff.result"
!Post File = "Stokes_prognostic_ELA400_SMB_noflow.vtu"
Post File = "Eismint_acculim.vtu"
Initialize Dirichlet Conditions = Logical False
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Constants !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Constants
Stefan Boltzmann = 5.67e-08
Gas Constant = Real 8.314 !Joule/mol x K
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Bodies !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
Name = "Glacier"
Body Force = 1
Equation = 1
Material = 1
Initial Condition = 1
End
Body 2
Name = "Surface"
Body Force = 2
Equation = 2
Material = 2
Initial Condition = 2
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Equations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
Name = "Equation1"
Convection = "computed"
Flow Solution Name = String "Flow Solution"
Active Solvers(6) = 1 2 4 5 6 7
End
Equation 2
Name = "Equation2"
Convection = "computed" !!! CHANGE TO THIS ONE TO GET REASONABLE RESULTS
!Convection = "none"
Active Solvers(1) = 3
Flow Solution Name = String "Flow Solution"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Initial Conditions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
Velocity 1 = 0.0
Velocity 2 = 0.0
Pressure = 0.0
Temp = Real 260.0
Mesh Velocity 1 = Real 0.0
Mesh Velocity 2 = Real 0.0
Mesh Velocity 3 = Real 0.0
! Free surface
Zs = Equals Coordinate 2
RefZs = Equals Coordinate 2
End
Initial Condition 2
!Free Surface
Zs = Equals Coordinate 2
RefZs = Equals Coordinate 2
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! Solvers !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Calculate the Height and Depth
Solver 1
Equation = "HeightDepth"
Exec Solver = "Before Timestep"
Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"
Active Coordinate = Integer 2
Operator 1 = depth
Operator 2 = height
End
!!! Navier-Stokes Equations
Solver 2
Equation = "Navier-Stokes"
Flow Model = String "Stokes"
Optimize Bandwidth = Logical True
Linear System Solver = Direct
Linear System Direct Method = "UMFPACK"
! Linear System Solver = "Iterative"
! Linear System Iterative Method = "GCR" !"BICGStab"
Linear System Max Iterations = 5000
Linear System Convergence Tolerance = 1.0E-06
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU1"
Linear System Residual Output = 1
Steady State Convergence Tolerance = 1.0E-03
! Stabilization Method = [Stabilized,P2/P1,Bubbles]
Stabilization Method = Stabilized
Nonlinear System Convergence Tolerance = 1.0E-05
Nonlinear System Convergence Measure = Solution
Nonlinear System Max Iterations = 100 ! 1 try without non-lin iters
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0E-03
Nonlinear System Reset Newton = Logical True
Exported Variable 1 = -dofs 3 "Mesh Velocity"
Exported Variable 2 = -dofs 3 "Mesh Update"
Nonlinear System Relaxation Factor = 0.8
End
!!! Solves for the free surface
Solver 3
Exec Solver = always
Equation = "Free Surface"
Variable = String "Zs"
Variable DOFs = 1
! needed for evaluating the contact pressure
Exported Variable 1 = -dofs 1 "Zs Residual"
! needed for storing the initial shape (needed for updates)
Exported Variable 2 = -dofs 1 "RefZs"
Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
! This would take the constrained points out of solution
! Use in serial run, only
! Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
Linear System Solver = Iterative
Linear System Max Iterations = 1500
Linear System Iterative Method = BiCGStab
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0e-7
Linear System Abort Not Converged = False
Linear System Residual Output = 1
Nonlinear System Max Iterations = 100
Nonlinear System Convergence Tolerance = 1.0e-6
Nonlinear System Relaxation Factor = 0.60
Steady State Convergence Tolerance = 1.0e-03
Stabilization Method = Bubbles
! Apply contact problem
!Apply Dirichlet = Logical True
! How much the free surface is relaxed
! Relaxation Factor = Real 0.9
End
!!! Maps the Structured Mesh
Solver 4
Exec Solver = "after timestep"
Equation = "MapCoordinate"
Procedure = "StructuredMeshMapper" "StructuredMeshMapper"
Active Coordinate = Integer 2 ! the mesh-update is y-direction
! For time being this is currently externally allocated
Mesh Velocity Variable = String "Mesh Velocity 2"
! The 1st value is special as the mesh velocity could be unrelistically high
Mesh Velocity First Zero = Logical True
Top Surface Variable = String "Zs"
Dot Product Tolerance = Real 0.01
End
!!! Deformational Heat
Solver 5
Equation = DeformationalHeat
Variable = String "StrHeat"
Variable DOFs = 1
procedure = "ElmerIceSolvers" "DeformationalHeatSolver"
Linear System Solver = direct
Linear System direct Method = umfpack
Recompute Stabilization = Logical True
End
Solver 6
! Exec Solver = "Never"
Equation = String "Homologous Temperature Equation"
Procedure = File "ElmerIceSolvers" "TemperateIceSolver"
! Comment next line in parallel, as EliminateDirichlet does
! not work in parallel
!------------------------------------------------------------
! Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
Variable = String "Temp"
Variable DOFs = 1
Stabilize = True
Optimize Bandwidth = Logical True
Linear System Solver = "Iterative"
Linear System Direct Method = UMFPACK
! 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 = "ILU1"
Linear System Residual Output = 0
Nonlinear System Convergence Tolerance = 1.0E-07
Nonlinear System Max Iterations = 100
Nonlinear System Relaxation Factor = Real 9.999E-01
Steady State Convergence Tolerance = 1.0E-04
! the contact algorithm (aka Dirichlet algorithm)
!-----------------------------------------------------
Apply Dirichlet = Logical True
! those two variables are needed in order to store
! the relative or homologous temperature as well
! as the residual
!-------------------------------------------------
Exported Variable 1 = String "Temp Homologous"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "Temp Residual"
Exported Variable 2 DOFs = 1
End
Solver 7
Equation = String "Divergence Solver"
Procedure = "DivergenceSolver" DivergenceSolver"
Target Variable = Velocity
Linear System Solver = "Iterative"
Linear System Iterative Method = "GCR"
Linear System Preconditioning = None
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! Materials !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
Name = "ice"
Density = Real $910.0*1.0E-06*(31556926.0)^(-2.0)
! Glen's flow law (using Glen)
!-----------------
! viscosity stuff
!----------------
!! call in scaled units (m-MPa-years)
Viscosity = Variable Temp Homologous
Real MATC "glen(tx)*31556926.0^(-1.0/3.0)*1.0E-06"
Critical Shear Rate = $.1E-09 * 31556926.0
!! this holds for both unit systems
Viscosity Model = String "power law"
Viscosity Exponent = $1.0/3.0
! Heat transfer stuff
Temp Heat Capacity = Variable Temp
Real MATC "capacity(tx)*(31556926.0)^(2.0)"
Temp Heat Conductivity = Variable Temp
Real MATC "conductivity(tx)*31556926.0*1.0E-06"
Temp Upper Limit = Variable Depth
Real MATC "273.15 - 9.8E-08 * tx * 910.0 * 9.81" !-> this is the correction of the presure melting point with respect to the hydrostatic overburden at the point
End
Material 2
Density = Real $910.0*1.0E-06*(31556926.0)^(-2.0)
Min Zs = Variable RefZs
Real MATC "tx - 0.1"
Max Zs = Variable RefZs
Real MATC "tx + 5000.0"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! Body Forces !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Name = "BodyForce1"
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real -9.7696e15 !gravity in MPa - a - m
!! passive below threshold for more stable ice front
Flow Solution Passive = Variable depth, height
Real MATC "((tx(0) + tx(1)) < 20.0)"
!! Strain Heating
Temp Volume Source = Equals StrHeat
End
Body Force 2
Name = "Climate"
Zs Accumulation Flux 1 = Real 0.0e0
Zs Accumulation Flux 2 = Variable Coordinate 1, Coordinate 2, time
Real MATC "accum(tx)"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! Boundary Conditions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Boundary Condition 1
Name = "bedrock"
Target Boundaries = 1
Bottom Surface = Equals Coordinate 2
Normal-Tangential Velocity = True
Flow Force BC = Logical True
Velocity 1 = Real 0.0e0
Velocity 2 = Real 0.0e0
Temp Flux BC = Logical True
Temp Heat Flux = Real $ 0.042 * (31556926.0)*1.0E-06 ! 20 mW m-2
End
Boundary Condition 2
Name = "sides"
Target Boundaries(2) = 3 4
Velocity 1 = Real 0.0e0
End
Boundary Condition 3
Name = "surface"
Top Surface = Equals "Zs"
Target Boundaries = 2
Body ID = 2 !!! THIS IS ESSENTIAL: the body the free surface solver is being rnu on
Temp = Variable Coordinate 1, Coordinate 2, time
Real MATC "surftemp(tx)"
End