glads_3dmesh.sif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Test GlaDS solvers              !!
!!                                 !!
!! Olivier Gagliardini June 2017   !!
!!                                 !!
!! SHMIP - test B5                 !!
!!                                 !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


check keywords warn
echo on

$Bnum = "B5"
$namerun = "glads_3dmesh_"

! All Units are in m, year, MPa
! For ice flow 
$yearindays = 365.00 
$yearinsec = yearindays*24*60*60 
$MPainPa = 1.0e6 

$Iter = 2
$dtIni = (1.0/yearindays)
$dtMin = 0.01/yearindays
$OutPut = 1

$ev = 0.0 
$Source = 0.0                 ! all input in moulins

! Common input parameters for all tests SHMIP 


$ub = 1.0e-6*yearinsec        !(1.0e-6 m/s)

!Prefactor from Cuffey and Paterson (2010) in MPa^{-3} a^{-1}
! temperate Ice (power law for SSA)
$rhoi = 910.0/(MPainPa*yearinsec^2) 
$rhow = 1000.0/(MPainPa*yearinsec^2) 
$A1 = 2.89165e-13*yearinsec*MPainPa^3
$A2 = 2.42736e-02*yearinsec*MPainPa^3

$ng = 3.0
$Aglen = 2.5e-25*yearinsec*MPainPa^3
$eta = (2.0*Aglen)^(-1.0/ng)

$gravity = 9.80*yearinsec^2  

! For the sheet
$Ar = Aglen
$alphas = 1.25 
$betas = 1.5 
$lr = 2.0 
$hr = 0.1 
$Ks = 0.005*yearinsec*(1.0/MPainPa)^(1.0-betas) 
$Hs = 0.05 ! IC for h

! 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 = 4.0

! Definition of the sqrt geometry 
! for SHMIP A, B, C, D

$Hmin = 1.0
$ function H(x) \
  import Hmin {\
  _H  = 6.0*(sqrt(x+5e3)-sqrt(5e3))+Hmin ;\
}


Header
  Mesh DB "." "mesh_B5_3d"
End

!---------------------------------------------------
!---------------- CONSTANTS ------------------------
!---------------------------------------------------

Constants
  Latent Heat = Real $Lw
  Gravity Norm = Real $gravity
  Fresh 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 "Zb"
End

!---------------------------------------------------
!---------------- SIMULATION -----------------------
!---------------------------------------------------

Simulation
  Coordinate System  = Cartesian 3D 
  Simulation Type = transient

  Timestepping Method = "bdf"
  BDF Order = 1

  Timestep Intervals(1) = $Iter
  Output Intervals(1) = $OutPut
  Timestep Sizes(1) = $dtIni 

  Steady State Max Iterations = 1 
  Steady State Min Iterations = 1
  
  Post File = $namerun$".vtu"

  max output level = 3
End

!---------------------------------------------------
!---------------- BODIES ---------------------------
!---------------------------------------------------

! This body is located at the ice/bed interface and will be used to solve 
! the sheet equation
Body 1
  Name= "Ice"
  Equation = 1
  Material = 1
  Body Force = 1
  Initial Condition = 1
End

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


!---------------------------------------------------
!---------------- INITIAL CONDITIONS ---------------
!---------------------------------------------------

!! for the sheet 
Initial Condition 1
  Sheet Thickness = Real $Hs
  Zs = Variable Coordinate 1
    Real MATC "H(tx)"
  Zb = Real 0.0 
  Effective Pressure = Variable Coordinate 1
     Real MATC "rhoi*gravity*H(tx)"
End

!---------------------------------------------------
!---------------- BODY FORCES ----------------------
!---------------------------------------------------

! source 1 cm/day
Body Force 1
  Hydraulic Potential Volume Source = Real $Source 
End

!---------------------------------------------------
!---------------- MATERIALS ------------------------
!---------------------------------------------------

!! ice material properties in MPa - m - a system 
Material 1
  Density = Real $rhoi          
  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       

  Sliding Velocity = Real $ub

  Bedrock Bump Length = Real $lr
  Bedrock Bump High = Real $hr
  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 Coordinate 1
     Real MATC "rhoi*gravity*H(tx)"
End

!---------------------------------------------------
!---------------- SOLVERS --------------------------
!---------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

  Activate Channels = Logical True
  Activate Melt from Channels = Logical 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
  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
End

Solver 2 
  ! Just here to declare the variable Sheet Thickness
  Equation = "GlaDS Thickness sheet"
  Procedure = "ElmerIceSolvers" "GlaDSsheetThickDummy"
  Variable = -dofs 1 "Sheet Thickness"
End

Solver 3 
  ! Just here for output and declare the Channel Area variable
  ! It is executed simultaneously to saving
  Exec Solver = After Saving
  Equation = "GlaDS Channel 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 = String "results"
  Channels OutPut File Name = String $namerun$"_channels"
End

Solver 4
  Exec Solver = Never
  Equation = "Navier-Stokes"
  Linear System Solver = Direct     

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

  Nonlinear System Max Iterations = 50
  Nonlinear System Convergence Tolerance  = 1.0e-5
  Nonlinear System Newton After Iterations = 50  
  Nonlinear System Newton After Tolerance = 1.0e-05
  Nonlinear System Relaxation Factor = 1.00
  Nonlinear System Reset Newton = Logical True
    
  Steady State Convergence Tolerance = Real 1.0e-4

  Stabilization Method = String Bubbles

  Exported Variable 1 = Flow Solution Loads[Stress Vector:2 CEQ Residual:1] 
  Exported Variable 2 = -dofs 1 "Zs"
  Exported Variable 3 = -dofs 1 "Zb"
  Calculate Loads = Logical True

  Flow Model = String "Stokes"
End

!---------------------------------------------------
!---------------- EQUATIONS ------------------------
!---------------------------------------------------
! Equation for the ice
Equation 1
  Active Solvers (2) = 3 4
End

Equation 2
  Active Solvers (2) = 1 2
End


!---------------------------------------------------
!---------------- BOUNDARY CONDITIONS --------------
!---------------------------------------------------

Boundary Condition 1
  Target Boundaries(2) = 1 3
  No Channel BC = Logical True
End

Boundary Condition 2
  Name = "Sym"
  Target Boundaries(1) = 2
  No Channel BC = Logical True
End

Boundary Condition 3
  Name = "front"
  Target Boundaries(1) = 4
  Hydraulic Potential = Real 0.0
  No Channel BC = Logical True
End

Boundary Condition 4
  Name = "bed"
  Target Boundaries(1) = 5
  Body Id = 2
  Velocity 1 = Real 0.0
  Velocity 2 = Real 0.0
  Velocity 3 = Real 0.0
End

Boundary Condition 5
  Name = "surface"
  Target Boundaries(1) = 6
End

! This is the boundary condition for the moulins, case B5
Boundary Condition 6
  Name = "moulins"
  Target Boundaries(100) = 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
  Moulin Storage = Logical True
  Moulin Area = Real $Am
  Moulin Flux = Real $0.9*yearinsec
End
Solver 1 :: Reference Norm = Real 8.28146474E+00 
Solver 1 :: Reference Norm Tolerance = Real 1E-04