!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