# # # ! 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