# 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