mlb_gw_tmc_pf_init.sif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Elmer Solver input file
!!
!! Vectorized Stokes solution for Midre Lovenbreen.
!! Including Semi-Lagrangian dating solver.
!!
!! All units are in m-MPa-year
!! Temperatures are in Kelvin
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! switch that on for debugging only
! ---------------------------------
! check keywords warn
! echo on
!!!!!!!!!!!!!!!!!!!!!!!!!
! DEFINITIONS used in run
!!!!!!!!!!!!!!!!!!!!!!!!!!
$name="S_1995_DEMs_gw_tmc_pf_init_full2"
$meshfile="outline62_lc75"
$element="p:1 b:1"
#MINH=20.0
#Tc = -1.0

#directmethod="mumps" 
!---LUA BEGIN
! assert(loadfile('./icematerial_SI.lua'))()
! assert(loadfile('./glads.lua'))()
!---LUA END


!echo on

! Note that the Mesh2MeshSolver requires solver specific
! mesh and hence this must be different than primary mesh,
! even symlink is ok.
$restartdir="outline62_lc75"
$restartfile="S_1995_DEMs_gw_glads_init.result"






! Temperature of the simulation in Celsius
! with the formula for A works only if T > -10


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! HEADER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Partitioned mesh
Header
  Mesh DB "." "$meshfile$"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! CONSTANTS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Constants
  Gravity(4) = 0 -1 0 #-gravity
  Gas Constant = Real 8.314 ! SI Units
  Stefan Boltzmann = Real 6.78E-08
  Permafrost Offset Density = Logical True
  Reference Temperature = Real 273.15
  Reference Pressure = Real #p0
  Gas Constant = Real #Rg
  Permafrost DeltaT = Real 4.0
  Permafrost eps = Real 0.99
  Avogadro Number = Real 6.022140857E23
  Undulation Amplitude = Real 0.0
  ! GlaDS
  Latent Heat = Real #Lw
  Gravity Norm = Real #-gravity
  Fresh Water Density = Real #rhow
  Water Density = Real #rhow
  Ice Density = Real #rhoi
  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 "bedrockDEM"
  ! Permafrost
  Permafrost Offset Density = Logical True
  Reference Temperature = Real #T0
  Reference Pressure = Real #p0
  Permafrost DeltaT = Real 1.0
  Permafrost eps = Real 0.99
End


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
  Coordinate System  = "Cartesian 3D"
  Simulation Type = "Steady State"             

  ! Internal extrusion parameters, may be altered.
!  Extruded Mesh Levels = Integer 9 
  Extruded Max Coordinate = Real 1.0		

  Extruded Mesh Levels = Integer 25
  Extruded BC Layers(2) = Integer 12 15
  Preserve Baseline = Logical True
  Extruded Baseline Layer = Integer 3
  Extruded Mesh Density = Variable Coordinate 1
    Real
      0.0 2.0
      0.39999 0.5
      0.49999 0.3
      0.5 0.3
      1.0 2.0
    End
    
  Simulation Type = "Steady State"
  
 ! Timestepping Method = "bdf"
 ! BDF Order = 1
 ! Timestep Intervals(1) = #Iter
  !Timestep Sizes(1) = #dtIni
  Output Intervals(1) = 1 ! #OutPut
 ! Timestep Size = Variable Timestep
 !    Real LUA "settimestepsize(tx[0])"

  ! Coupled iterations between different solvers
  !---------------------------------------------
  Steady State Max Iterations = 30
  Steady State Min Iterations = 1

  ! usually, Dirichlet BC's are initialized before everything else. Sometimes those 
  ! conditions are dependent on solutions of earlier solvers; next line ensures that
  ! this is not an issue.
  !-----------------------------------------
  Initialize Dirichlet Conditions = Logical False
  
  ! Output files
  ! ------------
  Post File = "$name$.vtu"
  !vtu: Save Bulk Only = Logical True

  !Scalars File = $name$.dat
  !scalars: Parallel Reduce = Logical True

  Output File = "$name$.result"

  ! Restart happens in Solver 1
  Restart Before Initial Conditions = Logical True
  !Interpolation Passive Coordinate = Integer 3 

  ! how verbose the solver should be
  !  3 = Only warnings
  ! 32 = Maximum verbosity
  !-------------------------------------------------------
  Max Output Level = 5
  
  Output Variable 1 = bedrockdem
  Output Variable 2 = surfacedem1995
  Output Variable 3 = velocity
  Output Variable 4 = pressure
  Output Variable 5 = flow solution loads
  Output Variable 6 = temperature
  Output Variable 7 = temperature boundary weights
  Output Variable 8 = temperature loads
  Output Variable 9 = relative temperature
  Output Variable 10 = temperature contact active
  Output Variable 11 = glaciated
  Output Variable 12 = depth
  Output Variable 13 = height
  Output Variable 14 = sheet thickness
  Output Variable 15 = water pressure
  Output Variable 16 = effective pressure
  Output Variable 17 = groundwater flux
  Output Variable 18 = PFTemp
  Output Variable 19 = GWPressure
  Output Variable 20 = PFTemp Loads
  Output Variable 21 = Xi

  Restart Variable 1 = bedrockdem
  Restart Variable 2 = surfacedem1995
  Restart Variable 3 = velocity
  Restart Variable 4 = pressure
  Restart Variable 5 = flow solution loads
  Restart Variable 6 = temperature
  Restart Variable 7 = temperature boundary weights
  Restart Variable 8 = temperature loads
  Restart Variable 9 = relative temperature
  Restart Variable 10 = temperature contact active
  Restart Variable 11 = glaciated
  Restart Variable 12 = depth
  Restart Variable 13 = height
  Restart Variable 14 = sheet thickness
  Restart Variable 15 = water pressure
  Restart Variable 16 = effective pressure

  Restart File = File "$restartfile$"
  Restart Position = 10
  Restart Time = 0.0
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SOLVER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




! Maps the constant-thickness mesh between given bedrock and surface topology
Solver 1
  Exec Solver = "before simulation"
  Equation = "MapCoordinate"
  Procedure = "StructuredMeshMapper" "StructuredMeshMapper"

  Active Coordinate = Integer 3

  Displacement Mode = Logical False
  Correct Surface = Logical True
  Minimum Height = Real #MINH

  Correct Surface Mask = String "Glaciated"
  Dot Product Tolerance = 1.0e-3

  ! These are needed to host the variables that are read in
  Variable = -nooutput "DummyMesh"

  Exported Variable 1 = "bedrockDEM"
  !Exported Variable 1 Mask = String "BedRock"

  Exported Variable 2 = "surfaceDEM1995"
  !Exported Variable 2 Mask = String "Surface"
  !Bubbles in global system = False
End

Solver 2
 Equation = "partial update"
 Procedure = "MeshSolve" "MeshSolver"
 Exec Solver = "Before Simulation"
 Variable = -dofs 3 "Partial Vertical MeshUpdate"
 Linear System Solver = Iterative

 Linear System Iterative Method = BiCGStabL
 Linear System Preconditioning = none
 !ILU1
 Linear System Max Iterations = 500
 Linear System Convergence Tolerance = 1.0e-12
 Linear System Residual Output = 50

 Steady State Convergence Tolerance = 1.0e-5
 !Bubbles in global system = False
End

! Computes height and depth assuming an extruded mesh.
Solver 3 
  Exec Solver = "before simulation" 
  Equation = "HeightDepth"
  Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"
  Active Coordinate = Integer 3
  Operator 1 = depth
  Operator 2 = height
  Dot Product Tolerance = 1.0e-3
  !Bubbles in global system = False
End 


! vectorized & threaded version of the Navier-Stokes solver.
Solver 4
  !Exec Solver = "Never"
  Equation = "Stokes-Vec"
  Procedure = "IncompressibleNSVec" "IncompressibleNSSolver"
  Div-Curl Discretization = Logical False

  !Variable = -dofs 4 "Flow Solution" [Velocity:3 Pressure:1]

  Optimize Bandwidth = Logical True
  Flow Model = Stokes
  Stabilization Method = String Stabilized

  !Non-linear iteration settings:
  !------------------------------ 
  Nonlinear System Max Iterations = 50
  Nonlinear System Convergence Tolerance  = 1.0e-5
  Nonlinear System Newton After Iterations = 5
  Nonlinear System Newton After Tolerance = 1.0e-1
  Nonlinear System Relaxation Factor = #2/3
  Nonlinear System Reset Newton = Logical True

  include linsys/gcr.sif

  ! Convergence on timelevel (not required here)
  !---------------------------------------------
  Steady State Convergence Tolerance = Real 1.0e-3

  !Relative Integration Order = -1
  Number of Integration Points = Integer 21 ! 21, 28, 44, 64, ...

  ! 1st iteration viscosity is constant
  Constant-Viscosity Start = Logical True

  Calculate Loads = True

  !Bubbles in global system = False
  Exported Variable 1 = -dofs 4 "Flow Solution Loads" 
End

!-----------------------------------------------------
! heat transfer limited by the pressure melting point
! as upper limit
!-----------------------------------------------------
Solver 5
  !Exec Solver = "Never"
  Equation = String "Homologous Temperature Equation"
  Procedure =  File "HeatSolveVec" "HeatSolver"
  ! Comment next line in parallel, as EliminateDirichlet does
  ! not work in parallel
  !------------------------------------------------------------
!  Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
  Variable = String "Temperature"
  Variable DOFs = 1
  Stabilization method = Stabilized
  Stabilize = True
  Optimize Bandwidth = Logical True
 ! Linear System Solver = "Direct"
 ! Linear System Direct Method = #directmethod
  Linear System Solver = "Iterative"
  Linear System Iterative Method = "GCR"
  Linear System GCR Restart = Integer 200
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0E-08
  Linear System Abort Not Converged = True
  Linear System Preconditioning = "ILU1"
  Linear System Residual Output = 1
  Nonlinear System Convergence Tolerance = 1.0E-06
  Nonlinear System Max Iterations = 50
  Nonlinear System Min Iterations = 2
  !Nonlinear System Relaxation Factor = Real 9.999E-01
  Steady State Convergence Tolerance = 1.0E-03
  !Steady State Relaxation Factor = Real #1.0/3.0
  ! the contact algorithm (aka Dirichlet algorithm)
  !-----------------------------------------------------
  Apply Limiter = Logical True
  Calculate Loads = Logical True
  Save Limiter = Logical True
  Calculate Boundary Weights = Logical True
  ! those two variables are needed in order to store
  ! the relative or homologous temperature as well
  ! as the residual
  !-------------------------------------------------
  Exported Variable 1 = -dofs 1 "Relative Temperature"
  Exported Variable 2 = -dofs 1 "Friction Loads"
  Nonlinear Update Exported Variables = Logical True
  Update Exported Variables = Logical True
  !Bubbles in global system = False
End

Solver 6
  Equation = "GlaDS Coupled sheet"
  Exec Solver = "never"
  Procedure = "ElmerIceSolvers" "GlaDSCoupledSolver"
  Variable = -dofs 1 "Hydraulic Potential"

  Use GroundedMask = Logical True
  Zero Sheet At GL = Logical False
  
  Activate Channels = Logical True ! True
  Activate Melt from Channels = Logical False !True
  Neglect Sheet Thickness in Potential = Logical True

  ! choices are EXPLICIT, CRANK-NICOLSON, IMPLICIT
  Channels Integration method = String "Crank-Nicolson"
  Sheet Integration method = String "Crank-Nicolson"

  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 2 "Sheet Discharge"
  Exported Variable 7 = -dofs 1 "Sheet Storage"
  !Exported Variable 8 = -dofs 1 "Flux from Moulins"

  Linear System Solver = Direct     

  !Replace UMFPACK with MUMPS IF you are using
  !multiple partitions 
  Linear System Direct Method = UMFPACK

  Nonlinear System Max Iterations = 10
  Nonlinear System Convergence Tolerance  = 1.0e-6 !1.0e-6
  Nonlinear System Relaxation Factor = 1.00

  Coupled Max Iterations = Integer 10
  Coupled Convergence Tolerance = Real 1.0e-3

  Steady State Convergence Tolerance = 1.0e-03
  !Bubbles in global system = False
End

Solver 7
  ! Just here to declare the variable Sheet Thickness
  Equation = "GlaDS Thickness sheet"
  Exec Solver = "Before timestep"
  Exec Solver = "Never"
  Procedure = "ElmerIceSolvers" "GlaDSsheetThickDummy"
  Variable = -dofs 1 "Sheet Thickness"
  Exported Variable 1 = -dofs 1 "Meltrate"

  Nonlinear Update Exported Variables = Logical True
  Update Exported Variables = Logical True

  !Bubbles in global system = False
End

Solver 8
  ! Just here for output and declare the Channel Area variable
  ! It is executed simultaneously to saving
  !Exec Solver = After Timestep
  !Exec Solver = "never"
  Equation = "GlaDS Channel OutPut"
  Exec Solver = "never"
  Exec Interval = #OutPut
  Procedure = "ElmerIceSolvers" "GlaDSchannelOut"
  Variable = -dofs 1 "Channel Area"
! Define that the variable is define on the edges only
  Element = "n:0 e:1"
  
  Exported Variable 1 = -dofs 1 "Channel Flux"

  VTU OutPutFile = Logical True 
  VTU BinaryFile = Logical False

  Channels OutPut Directory Name = File "$meshfile$/channels"
  Channels OutPut File Name = String "$name$_channels"
  !Bubbles in global system = False
End
!==============================================================================
! /// Heat transfer in bedrock ///
!==============================================================================
!==============================================================================
! /// Initialize porosity accordin to material file  ///
!==============================================================================
Solver 9
  Equation = "Porosity Initialization"
  Exec Solver = "Before Simulation"  
  Procedure = "ElmerIceSolvers" "PorosityInit"
  Porosity Variable = String "Porosity"
  Nodal Porosity = Logical True
  Bubbles in global system = False
  Element = $element  
End

Solver 10
  Equation = "Permafrost"
  !Exec Solver = "Never"
  Element = $element ! has to be same as all other permafrost solvers (IP variables!)
  !Exec Solver = "before simulation"
  !Exec Solver = "never"
  !Variable = "DummyTemp"
  Variable = -dofs 1 "PFTemp"
  Procedure = "ElmerIceSolvers" "PermaFrostHeatTransfer"
  !Linear System Solver = Direct
  !Linear System Direct Method = #directmethod
  Output Xi = Logical True
  
  Linear System Solver = Iterative
  Linear System Iterative Method  = IDRS
  Linear System Max Iterations = 5000
  Linear System Preconditioning = "ILU2"
  Nonlinear System Convergence Tolerance = 1.0e-06
  Nonlinear System Max Iterations = 1
  Nonlinear System Relaxation Factor = $2.0/3.0
  Linear System Convergence Tolerance = 1e-8
  Steady State Convergence Tolerance = 1e-03
  !Exported Variable 1 = -dofs 1 -IP "Xi"
  Exported Variable 1 = -dofs 1 "Salinity"
  Exported Variable 2 = -dofs 1 "Porosity"
  !Exported variable 3 = -dofs 1 "PFTemp Loads"
  Calculate Loads = True
  Initialize Steady State = Logical True
  !Groundwater Flux Variable = String "Groundwater Flux"
  Bubbles in global system = False
  Computed Groundwater Flux = Logical True
End

!==============================================================================
! /// Groundwater flow ///
!==============================================================================
Solver 11
  Equation = "Darcy"
  !Exec Solver = "Before Timestep"
  !Exec Solver = "Never"
  Element = $element ! has to be same as all other permafrost solvers (IP variables!)
  Variable = -dofs 1 "GWPressure"
  Procedure = "ElmerIceSolvers" "PermafrostGroundwaterFlow"  
  !Linear System Solver = Direct
  !Linear System Direct Method = #directmethod
  Linear System Solver = Iterative
  Linear System Iterative Method  = IDRS
  Linear System Preconditioning = "ILU2"
  Linear System Max Iterations = 5000
  Linear System Convergence Tolerance = 1e-09
  Nonlinear System Max Iterations = 30
  Nonlinear System Convergence Tolerance = 1e-05
  Steady State Convergence Tolerance = 1e-02
  Linear System Nullify Guess = Logical True
  Initialize Steady State = Logical True
  !Linear System Save = True
  !Linear System Save Continuous Numbering = Logical True
  !Exported Variable 1 = -dofs 1 Reference Pressure 
  Compute Freshwater Head = Logical True
  !Ground Stress Invariant Variable Name = String "Stress Invariant"
  !Ground Deformation Variable Name = String "bedrock deformation"
  !Steady State Relaxation Factor = $1/10
  Exported Variable 1 = -dofs 3 -IP "Groundwater Flux"
  !Project Type 1 = String "nodal"
  !Project Variable 1 = String "Nodal Groundwater Flux"
  Exported Variable 2 = -dofs 1 "Reference GWPressure"
  Exported Variable 3 = -dofs 1 -IP "Reference Offset Density"
  Exported Variable 4 = -dofs 1 -IP "Freshwater Head"
  Exported Variable 5 = -dofs 3 "Nodal Groundwater Flux"
  !Exported Variable 7 = -dofs 1 -IP "Xi"
  Bubbles in global system = False
  Update Exported Variables = Logical True
End


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODIES (i.e., domains to compute on)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
  Name= "bedrock"
  Equation = 1
  Material = 1
  Initial Condition = 1
  Body Force = 1
End

Body 2
  Name= "sediment"
  Equation = 1
  Material = 2
  Initial Condition = 1
  Body Force = 1
End

Body 3
  Name = "glacier"
  Equation = 2
  Material = 3
  Body Force = 2
  Initial Condition = 2
End

Body 4
  Name = "bedrock interface"
  Equation = 3
  Material = 3
  Body Force = 2
  Initial Condition = 2
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! EQUATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
  Name = "bedrock equations"
  Active Solvers(5) =  1 2 9 10 11
  Convection = "none"
End

Equation 2
  Name = "ice equations"
  Active Solvers(5) = 1 2 3 4 5
  Convection = "Computed"
  Flow Solution Name = String "Flow Solution"
End

Equation 3
  Name = "interface equations"
  Active Solvers(3) =  6 7 8
  Convection = "Computed"
  Flow Solution Name = String "Flow Solution"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! INITIAL CONDITIONS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
 PFTemp = Variable height
   Real Lua "tx[0]*(274.15 - 279.15)/110.0 + 279.15"
 Salinity = Real 0.0
! GWPressure = Real 0.0
 GWPressure = Variable height
     Real LUA "-gravity*rhow*(110.0 - tx[0]) + p0"
! Sheet Thickness = Real #Hs
End

Initial Condition 2
  ! Initial velocity is not really needed when using newtonian start
  !Velocity 1 = 0.0
  !Velocity 2 = 0.0
  !Velocity 3 = 0.0
  !Pressure = 0.0
  !Temperature = Variable Coordinate 3 ! lapserate of -10/1000 C/m
  !    Real
  !          0.0   272.14
  !	  100.0   272.15
  !       1000.0   263.15
  !    End
  !Sheet Thickness = Real #Hs
End


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODY FORCE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force  1
  Stress Bodyforce at IP = Logical True
  Stress Bodyforce 2 = Variable PFTemp, GWPressure, Porosity, Salinity, Xi
    Real Procedure "ElmerIceSolvers" "GetElasticityForce"
  Stress Bodyforce 1 = 0.0
  Stress Pressure = Equals GWPressure
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
End

Body Force 2
  Flow BodyForce 1 = 0.0
  Flow BodyForce 2 = 0.0 
  Flow BodyForce 3 = #gravity
  Temperature Upper Limit = Real 273.15
  Temperature = Real 273.15
  Temperature Condition = Variable Coordinate 3 
      Real
            0.0      -1.0
	  399.9999   -1.0
	  400.        1.0
         1000.0       1.0
      End
  Relative Temperature = Variable Temperature
     Real Lua "min(tx[0] - 273.15, 0.0)"

  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Channel Area = Real 0.0
  Channel Area Condition = Opposes glaciated
  Channel Flux = Real 0.0
  Channel Flux Condition = Opposes glaciated
  Channel Area Passive = Opposes glaciated
  Channel Flux Passive = Opposes glaciated
  Hydraulic Potential = Real 0.0
  !Hydraulic Potential = Variable Coordinate 3
  !   Real Lua "-rhow*gravity*tx[0]"
  Hydraulic Potential Condition = Opposes glaciated
  Hydraulic Potential Passive = Opposes glaciated
  Sheet Thickness = Real 0.001
  Sheet Thickness Condition = Opposes glaciated

  Hydraulic Potential Volume Source = Variable "temperature contact load", "temperature boundary weights", "temperature contact active"
  Real LUA "getmeltrate(tx[0],tx[1],tx[2]) + Source*tx[2]"
   Meltrate =  Variable "temperature contact load", "temperature boundary weights", "temperature contact active"
       Real LUA "getmeltrate(tx[0],tx[1],tx[2]) + Source*tx[2]"
 End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MATERIAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
  Name = "rock material"

  Rock Material File = File "./toymodel.db"
  Rock Material ID = Integer 2
  Hydraulic Conductivity Limit = Real 1.0E-14
  Swap Tensor = Logical True
  Conductivity Arithmetic Mean Weight = Real 0.5
  Radiogenic Reference Depth = Real 10000.0
  Youngs Modulus = Variable Porosity, Xi
     Real Procedure "ElmerIceSolvers" "GetEG"
  Poisson Ratio =   Variable Porosity, Xi
     Real  Procedure "ElmerIceSolvers" "GetNuG"
  Youngs Modulus at IP = Logical True
  Poisson Ratio at IP = Logical True
  Porosity Variable = Equals "Porosity"
  Temperature Variable = Equals "PFTemp"
  Salinity Variable = Equals "Salinity" 
  Pressure Variable = Equals "GWPressure"
  
  Mesh Elastic Modulus = Real 100.0
  Mesh Poisson Ratio = Real 0.0

  Viscosity = 10.0
  Density = 3000.0
End
Material 2
  Name = "sediment material"

  Rock Material File = File "./toymodel.db"
  Rock Material ID = Integer 1
  Hydraulic Conductivity Limit = Real 1.0E-14
  Swap Tensor = Logical True
  Conductivity Arithmetic Mean Weight = Real 0.5
  Radiogenic Reference Depth = Real 10000.0
  Youngs Modulus = Variable Porosity, Xi
     Real Procedure "ElmerIceSolvers" "GetEG"
  Poisson Ratio =   Variable Porosity, Xi
     Real  Procedure "ElmerIceSolvers" "GetNuG"
  Youngs Modulus at IP = Logical True
  Poisson Ratio at IP = Logical True
  Porosity Variable = Equals "Porosity"
  Temperature Variable = Equals "PFTemp"
  Salinity Variable = Equals "Salinity" 
  Pressure Variable = Equals "GWPressure"
  Mesh Elastic Modulus = Real 100.0
  Mesh Poisson Ratio = Real 0.0


  Mesh Elastic Modulus = Real 100.0
  Mesh Poisson Ratio = Real 0.0

  Viscosity = 10.0
  Density = 3000.0
 ! Nonnewtonian viscosity
  Viscosity Model = String Glen
  Glen Exponent = Real 3.0
  Critical Shear Rate =  Real 1.0E-16
  ! Paterson value in MPa^-3a^-1
  Limit Temperature = Real -10.0
  ! those are already transferred to Mpa-m-a system
  Rate Factor 1 = Real #A1 
  Rate Factor 2 = Real #A2 
  Activation Energy 1 = Real #Q1    
  Activation Energy 2 = Real #Q2
  !------------------------------------
  Glen Enhancement Factor = Real 1.0
  ! GlaDS
  ! For the sheet 
  Sheet Conductivity = Real #Ks 
  Sheet flow exponent alpha = Real #alphas
  Sheet flow exponent beta = Real #betas
  Englacial Void Ratio = Real #ev       

  Sliding Velocity = Variable Velocity 1, Velocity 2, Velocity 3
     Real LUA "math.sqrt(tx[0]*tx[0] + tx[1]*tx[1] + tx[2]*tx[2])"

  Bedrock Bump Length = Real #lr
  Bedrock Bump High = Real #hr
  !Sheet Closure Coefficient = Real #Ar
  Sheet Closure Coefficient = Variable "Relative Temperature"
     Real Lua "arrhenius(tx[0])"

! For the Channels
  Channel Conductivity = Real #Kc 
  Channel flow exponent alpha = Real #alphac
  Channel flow exponent beta = Real #betac
  !Channel Closure Coefficient = Real #Ac
   Channel Closure Coefficient = Variable "Relative Temperature"
     Real Lua "arrhenius(tx[0])"
  Sheet Width Over Channel = Real #lc
  Pressure Melting Coefficient = Real #Ct
  Water Heat Capacity = Real #Cw

! For both
  Ice Normal Stress = Variable depth
     Real Lua "-rhoi*gravity*tx[0]"
End

Material 3
  Name = "Ice"
  Density = Real #rhoi

  ! First viscosity with newtonian fluid
  ! happens to give velocities of proper size
  Viscosity = Real 1.0E+15 

  ! Nonnewtonian viscosity
  Viscosity Model = String Glen
  Glen Exponent = Real 3.0
  Critical Shear Rate =  Real 1.0E-16
  ! Paterson value in MPa^-3a^-1
  Limit Temperature = Real -10.0
  ! those are already transferred to Mpa-m-a system
  Rate Factor 1 = Real #A1 
  Rate Factor 2 = Real #A2 
  Activation Energy 1 = Real #Q1    
  Activation Energy 2 = Real #Q2
  !------------------------------------
  Glen Enhancement Factor = Real 1.0
  Relative Temperature = Variable Temperature
    Real lua "tx[0] - 273.15"
  ! Heat transfer stuff
  Heat Capacity = Variable Temperature
    Real lua "capacity(tx[0])"
  

   Heat Conductivity = Variable Temperature
    Real lua "conductivity(tx[0])"
    
  Pressure Scaling Factor = Real 1.0
  Temperature Upper Limit = Real 273.15 ! we ignore pressure melting point

  ! GlaDS
  ! For the sheet 
  Sheet Conductivity = Real #Ks 
  Sheet flow exponent alpha = Real #alphas
  Sheet flow exponent beta = Real #betas
  Englacial Void Ratio = Real #ev       

  Sliding Velocity = Variable Velocity 1, Velocity 2, Velocity 3
     Real LUA "math.sqrt(tx[0]*tx[0] + tx[1]*tx[1] + tx[2]*tx[2])"

  Bedrock Bump Length = Real #lr
  Bedrock Bump High = Real #hr
  !Sheet Closure Coefficient = Real #Ar
  Sheet Closure Coefficient = Variable "Relative Temperature"
     Real Lua "arrhenius(tx[0])"
! For the Channels
  Channel Conductivity = Real #Kc 
  Channel flow exponent alpha = Real #alphac
  Channel flow exponent beta = Real #betac
  !Sheet Closure Coefficient = Real #Ar
  !Channel Closure Coefficient = Real #Ac  
  Channel Closure Coefficient = Variable "Relative Temperature"
     Real Lua "arrhenius(tx[0])"
  Sheet Width Over Channel = Real #lc
  Pressure Melting Coefficient = Real #Ct
  Water Heat Capacity = Real #Cw

! For both
  Ice Normal Stress = Variable depth
     Real Lua "-rhoi*gravity*tx[0]"
End
!---------------------------------------------------
!---------------- BOUNDARY CONDITIONS --------------
!---------------------------------------------------
Boundary Condition 1
  ComputeNormal = Logical False
  Name = "sideline"
  Target Boundaries(1) = 1
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  No Channel BC = Logical True
  Hydraulic Potential =  Real 0.0
  Hydraulic Potential Conditions = Opposes  glaciated
End

Boundary Condition 2
  ComputeNormal = Logical False
  Name = "sidesline2"
  Target Boundaries(1)  = 2
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  No Channel BC = Logical True
  Hydraulic Potential =  Real 0.0
  Hydraulic Potential Conditions = Opposes glaciated
End

Boundary Condition 3
  Name="side rock1"
  Target Boundaries(1) = 3
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0 
End

Boundary Condition 4
  Name="side rock2"
  Target Boundaries(1) = 4
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0 
End
Boundary Condition 5
  Name="side sediment1"
  Target Boundaries(1) = 5
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  No Channel BC = Logical True
End

Boundary Condition 6
  Name="side sediment2"
  Target Boundaries(1) = 6
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  No Channel BC = Logical True
End

Boundary Condition 7
  ComputeNormal = Logical False
  Name = "side ice1"
  Target Boundaries(1) = 7
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Partial Vertical MeshUpdate 3 = Real 0.0
  No Channel BC = Logical True
End

Boundary Condition 8
  ComputeNormal = Logical False
  Name = "side ice2"
  Target Boundaries(1)  = 8
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Partial Vertical MeshUpdate 3 = Real 0.0
  No Channel BC = Logical True
End

!! DON'T CHANGE ORDER OF NEXT 2 BC's!
!! They are automaticaly created in internal extrusion
!! bedrock:
Boundary Condition 9
  Name = "bedrock lower"
  !Heat Flux = Real 0.05
  Permafrost Target Body = Integer 1
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Partial Vertical MeshUpdate 3 = Real 0.0
  Bottom Surface = Variable bedrockDEM
    Real Lua "tx[0] - 100.0"

  !PFTemp = Real 277.15

  Heat Flux = Real 0.05
End
Boundary Condition  10
  Name = "rock - sediment"
  Partial Vertical MeshUpdate 3 = Variable Coordinate 3, bedrockDEM
     Real Lua "tx[1] - tx[0] - sedimentheight"
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
End
Boundary Condition 11
  ComputeNormal = Logical True
  Name = "bedrock"

  Body ID = 4

  ! No-slip velocity conditions
  !Velocity 1 = Real 0.0
  !Velocity 2 = Real 0.0
  !Velocity 3 = Real 0.0

  ! slip conditions
  Normal-Tangential Velocity = True
  Velocity 1 = Real 0.0
  Weertman Friction Coefficient = Real #0.1*(Pa2MPa^(-1))*yearinsec^(1.0/3.0)
  Weertman Exponent = Real #1.0/3.0
  Weertman Linear Velocity = Real 1.0E-10
  Friction Newton Epsilon = Real 1.0E-03
  Friction Newton Linearization = Logical True
  
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Partial Vertical MeshUpdate 3 = Real 0.0

  Mid Surface = Equals "bedrockDEM"

! Mask for creating the restart fields only where needed
! Refererred by: Exported Variable 1 Mask = ...
  Bedrock = Logical True
  Temperature Flux BC = Logical True
  !Heat Flux = Real #0.050  !!!! SWITCH OFF FOR COUPLES PF an TMC
  ! 50 mW m^-2
  Friction Loads = Variable "Flow Solution Loads 1",  "Flow Solution Loads 2",  "Flow Solution Loads 3", "Velocity 1", "Velocity 2", "Velocity 3"
    Real Lua "frictionloads(tx[0],tx[1],tx[2],tx[3],tx[4],tx[5])"
  Temperature Load = Variable "Flow Solution Loads 1",  "Flow Solution Loads 2",  "Flow Solution Loads 3", "Velocity 1", "Velocity 2", "Velocity 3", "PFTemp Loads"
    Real Lua "frictionloads(tx[0],tx[1],tx[2],tx[3],tx[4],tx[5]) - tx[7]"
  Calculate Temperature Boundary Weights = Logical True

  Permafrost Target Body = Integer 2
  GWPressure = Variable Water Pressure
     Real Lua "max(tx[0],p0)"

  External Temperature = Equals "Temperature"
  Heat Transfer Coefficient = Variable "Sheet Thickness", "water pressure", "depth"
      Real Lua "min(htcoeff(kappaw0,tx[0]), kappaw0/0.001)"
  !    Real Lua "max(htcoeff(kappaw0,tx[0])*tx[1]/(-tx[2]*gravity*rhoi),0.00001)"
  !PFTemp = Real #TSurf0
  !PFTemp Condition =  Variable depth
  !  Real LUA "1.01*MINH - tx[0]"

End

Boundary Condition 12
  ComputeNormal = Logical False
  Name = "surface"
  Top Surface = Equals "surfaceDEM1995"

! Referred by: Exported Variable 2 Mask = ...
  Surface = Logical True
  !Target Boundaries(1) = 4
  Temperature = Variable Coordinate 3 ! lapserate of -10/1000 C/m
      Real
            0.0   273.15
	  100.0   272.14
         1000.0   263.15
      End
  Partial Vertical MeshUpdate 1 = Real 0.0
  Partial Vertical MeshUpdate 2 = Real 0.0
  Partial Vertical MeshUpdate 3 = Real 0.0
End