!!--------------------------------------------------------!! ! Matt Trevers - auto generated .sif file for frontal buoyancy ! Based on the MISMIP SETUP ! Experiment MISMIP 1a - step 4 ! starting from the Schoof solution !!--------------------------------------------------------!! check keywords warn echo on $namerun = "dw900_h980_g3_n0_coulomb_" ! ! working units are MPa, a, m ! $yearinsec = 365.25*24*60*60 $rhoi = 918.0/(1.0e6*yearinsec^2) $rhow = 1028.0/(1.0e6*yearinsec^2) $A = 4.6416e-25*yearinsec*1.0e18 $n = 3.0 $eta = 1.0/(2.0*A)^(1.0/n) $gravity = -9.8*yearinsec^2 $As = 1*1.0e18*yearinsec !Change this value to closely match the velocity profile to Coulomb ! The cryostatic pressure at the back of the glacier ! For gradient = 3 degrees, ht = 980.00m, length = 10000.00m ! Height at rear = 1504.08m (assume this isn't the actual back $function cryopressure(y) import rhoi, gravity {\ depth=604.08-y(0);\ _cryopressure=abs(depth*rhoi*gravity);\ if (y(0)<-900.00)\ {_cryopressure=0.0;}\ if (y(0)>604.08)\ {_cryopressure=0.0;}\ } $function velocityin(y) import rhoi, gravity, n {\ u0 = 0;\ B = 200000.0;\ alpha = tan(4.0*pi/180.0);\ hL = 1504.08;\ hLZ = (604.08-y(0));\ _velocityin=u0+0.5*(((rhoi*gravity*alpha/B)^3)*((hL^4)-(hLZ^4)));\ } Header Mesh DB "." "dw900_h980_g3_n0" End Constants Water Density = Real $rhow End !--------------------------------------------------- !---------------- SIMULATION ----------------------- !--------------------------------------------------- Simulation Coordinate System = Cartesian 2D Simulation Type = Steady state Timestepping Method = "bdf" BDF Order = 1 Timestep Intervals = 200 Output Intervals = 1 Timestep Sizes = 0.005 Initialize Dirichlet Conditions = Logical False Steady State Max Iterations = 1 Steady State Min Iterations = 1 Post File = "$namerun".vtu" max output level = 3 End !--------------------------------------------------- !---------------- BODIES --------------------------- !--------------------------------------------------- ! the ice Body 1 Name = "ice" Equation = 1 Body Force = 1 Material = 1 Initial Condition = 1 End ! The upper surface Body 2 Name= "top free surface" Equation = 2 Material = 1 Body Force = 2 Initial Condition = 2 End ! the lower surface Body 3 Name= "free surface sea/ice-shelf" Equation = 3 Material = 1 Body Force = 3 Initial Condition = 3 End !--------------------------------------------------- !---------------- INITIAL CONDITIONS --------------- !--------------------------------------------------- !! for ice Initial Condition 1 Pressure = Real 0.0 Velocity 1 = Real 0.0 Velocity 2 = Real 0.0 Bedrock = Real -900.0 End !! for top free surface Initial Condition 2 Zs = Equals Coordinate 2 End !! for free surface sea/ice-shelf Initial Condition 3 Zb = Equals Coordinate 2 End !--------------------------------------------------- !---------------- BODY FORCES ---------------------- !--------------------------------------------------- Body Force 1 Flow BodyForce 1 = Real 0.0 Flow BodyForce 2 = Real $gravity End !! accumulation flux in m/year - no accumulation! Body Force 2 Zs Accumulation Flux 1 = Real 0.0e0 Zs Accumulation Flux 2 = Real 0.0e0 !m/a End !! no melting/accretion under ice/shelf Body Force 3 Zb Accumulation = Real 0.0e0 ! BDrag = Variable Coordinate 2 ! Real Procedure "ElmerIceUSF" "Friction_Coulomb" End !--------------------------------------------------- !---------------- MATERIALS ------------------------ !--------------------------------------------------- !! ice material properties in MPa - m - a system Material 1 ! Viscosity Model = String "power law" ! Density = Real $rhoi ! Viscosity = Real $eta ! Viscosity Exponent = Real $1.0/n ! Critical Shear Rate = Real 1.0e-15 Density = Real $rhoi ! Viscosity Model = String "power law" ! Viscosity = Real $eta ! Viscosity Exponent = Real $1.0/n Critical Shear Rate = Real 1.0e-15 !---------------- ! viscosity stuff !---------------- Viscosity Model = String "Glen" ! Viscosity has to be set to a dummy value ! to avoid warning output from Elmer Viscosity = Real 1.0 ! Dummy variable Glen Exponent = Real 3.0 ! Rate factors (Paterson value in MPa^-3a^-1) Rate Factor 1 = Real 1.258e13 ! 1.258e13 Rate Factor 2 = Real 6.046e28 ! 8e27 ! these are in SI units - no problem, as long as ! the gas constant also is Activation Energy 1 = Real 60.0e3 Activation Energy 2 = Real 139.0e3 Glen Enhancement Factor = Real 1.0 ! the variable taken to evaluate the Arrhenius law ! in general this should be the temperature relative ! to pressure melting point. The suggestion below plugs ! in the correct value obtained with TemperateIceSolver ! Temperature Field Variable = String "Temp Homologous" ! the temperature to switch between the ! two regimes in the flow law Limit Temperature = Real -10.0 ! In case there is no temperature variable Constant Temperature = Real -9.0 Sea level = Real 0.0 Min Zs = Variable "Bottom Zb" Real MATC "tx + 10.0" Max Zs = Real 1.0e6 !! Bed condition Min Zb = Equals Bedrock Max Zb = Real 0.0 Cauchy = Logical True ! Youngs Modulus = 9e9 Poisson Ratio = 0.31 End !--------------------------------------------------- !---------------- SOLVERS -------------------------- !--------------------------------------------------- !! Initialisation of the Grounded Mask Solver 1 Exec Solver = Before All Equation = GroundedMaskIni Procedure = "ElmerIceSolvers" "GroundedSolver" Variable = GroundedMask Variable DOFs = 1 Toler = Real 1.0e-3 Bedrock Variable = String "Bedrock" End Solver 2 Exec Solver = "Never" Equation = "MapCoordinate" Procedure = "StructuredMeshMapper" "StructuredMeshMapper" Active Coordinate = Integer 2 Mesh Velocity Variable = String "dSdt" Mesh Update Variable = String "dS" Mesh Velocity First Zero = Logical True Top Surface Variable Name = String "Zs" Bottom Surface Variable Name = String "Zb" Displacement Mode = Logical False Correct Surface = Logical True Minimum Height = Real 1.0 End Solver 3 Equation = "NormalVector" Procedure = "ElmerIceSolvers" "ComputeNormalSolver" Variable = String "Normal Vector" Variable DOFs = 2 ComputeAll = Logical False Optimize Bandwidth = Logical False End Solver 4 Equation = Fw Procedure = "ElmerIceSolvers" "GetHydrostaticLoads" Variable = Fw[Fwater:2] Variable DOFs = 2 End Solver 5 Equation = "Navier-Stokes" Linear System Solver = Direct Linear System Direct Method = Mumps Nonlinear System Max Iterations = 150 Nonlinear System Convergence Tolerance = 1.0e-06 Nonlinear System Newton After Iterations = 50 Nonlinear System Newton After Tolerance = 1.0e-06 Nonlinear System Relaxation Factor = 1.0 Nonlinear System Reset Newton = Logical True Steady State Convergence Tolerance = Real 1.0e-05 Stabilization Method = String Bubbles Exported Variable 1 = Flow Solution Loads[Stress Vector:2 CEQ Residual:1] Calculate Loads = Logical True Exported Variable 2 = -dofs 1 "dSdt" Exported Variable 3 = -dofs 1 "dS" Exported Variable 4 = -dofs 1 "Bedrock" Flow Model = String "Stokes" End Solver 6 Equation = String "StressSolver" Procedure = File "ElmerIceSolvers" "ComputeDevStress" ! this is just a dummy, hence no output is needed !----------------------------------------------------------------------- Variable = -nooutput "Sij" Variable DOFs = 1 ! the name of the variable containing the flow solution (U,V,W,Pressure) !----------------------------------------------------------------------- Flow Solver Name = String "Flow Solution" ! no default value anymore for "Stress Variable Name" !----------------------------------------------------------------------- Exported Variable 1 = "Stress" ! [Sxx, Syy, Szz, Sxy] in 2D ! [Sxx, Syy, Szz, Sxy, Syz, Szx] in 3D Exported Variable 1 DOFs = 4 ! 4 in 2D, 6 in 3D Stress Variable Name = String "Stress" Linear System Solver = "Iterative" Linear System Iterative Method = "BiCGStab" Linear System Max Iterations = 300 Linear System Convergence Tolerance = 1.0E-06 Linear System Abort Not Converged = True Linear System Preconditioning = "ILU0" Linear System Residual Output = 1 End Solver 7 Equation = "HeightDepth" Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane" Active Coordinate = Integer 2 Dot Product Tolerance = Real 1.0e-4 Operator 1 = Depth Operator 2 = Height ! Export Zb on the Upper surface Variable 3 = Zb Operator 3 = Bottom End Solver 8 Exec Solver = "Never" Equation = "Free Surface Top" Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver" Variable = "Zs" Variable DOFs = 1 Exported Variable 1 = "Zs Residual" Exported Variable 1 DOFs = 1 Before Linsolve = "EliminateDirichlet" "EliminateDirichlet" Linear System Solver = Iterative Linear System Direct Method = UMFPACK Linear System Max Iterations = 1 Linear System Iterative Method = BiCGStab Linear System Preconditioning = ILU0 Linear System Convergence Tolerance = Real 1.0e-06 Linear System Abort Not Converged = False Linear System Residual Output = 1 Nonlinear System Max Iterations = 100 Nonlinear System Convergence Tolerance = 1.0e-05 Nonlinear System Relaxation Factor = 0.1 Steady State Convergence Tolerance = 1.0e-03 Stabilization Method = Stabilized Apply Dirichlet = Logical True Relaxation Factor = Real 0.1 End Solver 9 Exec Solver = "Never" Equation = "Free Surface Sea/Shelf" Procedure = "./MyFreeSurfaceSolver" "FreeSurfaceSolver" Variable = "Zb" Variable DOFS = 1 Exported Variable 1 = "Zb Residual" Exported Variable 1 DOFs = 1 Before Linsolve = "EliminateDirichlet" "EliminateDirichlet" Linear System Solver = Iterative Linear System Direct Method = UMFPACK Linear System Max Iterations = 1 Linear System Iterative Method = BiCGStab Linear System Preconditioning = ILU0 Linear System Convergence Tolerance = Real 1.0e-06 Linear System Abort Not Converged = False Linear System Residual Output = 1 Nonlinear System Max Iterations = 100 Nonlinear System Convergence Tolerance = 1.0e-04 Nonlinear System Relaxation Factor = 0.1 Steady State Convergence Tolerance = 1.0e-02 Stabilization Method = Stabilized Apply Dirichlet = Logical True Relaxation Factor = Real 0.1 End !! Compute the Mask Solver 10 Equation = GroundedMask Procedure = "ElmerIceSolvers" "GroundedSolver" Variable = GroundedMask Variable DOFs = 1 Toler = Real 1.0e-3 Bedrock Variable = String "Bedrock" End !--------------------------------------------------- !---------------- EQUATIONS ------------------------ !--------------------------------------------------- Equation 1 Active Solvers (5) = 2 3 5 6 7 End Equation 2 Active Solvers(1) = 8 Flow Solution Name = String "Flow Solution" Convection = String Computed End Equation 3 Active Solvers(4) = 1 4 9 10 Flow Solution Name = String "Flow Solution" Convection = String Computed End !--------------------------------------------------- !---------------- BOUNDARY CONDITIONS -------------- !--------------------------------------------------- !! BC Bedrock + Shelf Boundary Condition 1 Name = "bottom" Target Boundaries = 1 Body Id = 3 Normal-Tangential Velocity = Logical True Flow Force BC = Logical True ! ! Condition where the bed is stuck ! Zb = Equals Bedrock Zb Condition = Variable GroundedMask Real MATC "tx + 0.5" ! ! Bedrock conditions ! Slip Coefficient 2 = Variable Coordinate 1 Real Procedure "ElmerIceUSF" "SlidCoef_Contact" Sliding Law = String "Coulomb" !! Parameters needed for the Coulomb Friction Law Friction Law Sliding Coefficient = Real $As Friction Law Post-Peak Exponent = Real 1.0 !(q=1) Friction Law Maximum Value = Real 1.0 !(C=1) Friction Law PowerLaw Exponent = Real 3.0 !(m = n = 3 Glen's law) Friction Law Linear Velocity = Real 1.0 ! Options are 'Last Grounded' (default), 'First Floating' or 'Discontinuous' Grounding Line Definition = String "Discontinuous" Test Contact Tolerance = real 1.0e-3 Velocity 1 = Real 0.0 Velocity 1 Condition = Variable GroundedMask Real MATC "tx + 0.5" ! ! Shelf conditions ! External Pressure = Variable Coordinate 2 Real Procedure "ElmerIceUSF" "SeaPressure" Slip Coefficient 1 = Variable Coordinate 2 Real Procedure "ElmerIceUSF" "SeaSpring" ComputeNormal Condition = Variable GroundedMask Real MATC "tx + 0.5" Compute Sea Pressure = Logical True Compute Sea Spring = Logical True End !! BC Lateral Ice-Shelf (air or sea contact) Boundary Condition 2 Name = "front" Target Boundaries = 4 Flow Force BC = Logical True External Pressure = Variable Coordinate 2 Real Procedure "ElmerIceUSF" "SeaPressure" Compute Sea Pressure = Logical True ComputeNormal = Logical False End !! BC Free surface Top Boundary Condition 3 Name = "top" Target Boundaries = 3 Body Id = 2 ComputeNormal = Logical False Depth = Real 0.0 End !! Symetry axis Boundary Condition 4 Name = "back" Target Boundaries = 2 Velocity 1 = Variable Coordinate 2 Real MATC "-velocityin(tx)" ComputeNormal = Logical False External Pressure = Variable Coordinate 2 Real MATC "cryopressure(tx)" End