ismip_SSA_3D.sif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                       !!
!! ISMIP-HOM D020 solved in prognostic   !!
!!   using SSA and thickness solvers     !!
!!           done with a 3D Mesh         !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


$Step = "SSA_3D"
$L = 20.0e3
$Slope = 0.1 * pi / 180.0

$yearinsec = 365.25*24*60*60
$rhoi = 900.0/(1.0e6*yearinsec^2)   
$gravity = -9.81*yearinsec^2
$n = 3.0
$eta = (2.0*100.0)^(-1.0/n)

Header
  Mesh DB "." "rectangle"
End

Constants
! No constant Needed
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
  Coordinate System  = Cartesian 2D 
  Simulation Type = Transient     

  Timestepping Method = "bdf"
  BDF Order = 1
  Output Intervals = 1
  Timestep Intervals =  100 
  Timestep Sizes = 1.0 

  Steady State Min Iterations = 1
  Steady State Max Iterations = 1

  Post File = "test_$Step$.vtu"

  max output level = 3

!!! If using a 3D mesh it as to be extruded along vertical direction
!!! and BC of the 2D mesh (baseline) as to be preserved to impose Neumann BC cond. to the SSA.
  Extruded Mesh Levels = 3
  Preserve Baseline = Logical True


End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The bulk
Body 1
  Equation = 1
  Body Force = 1
  Material = 1
  Initial Condition = 1
End

! the bottom BC where SSA and Thickness are solved
Body 2
  Equation = 2
  Body Force = 1
  Material = 1
  Initial Condition = 1
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
   SSAVelocity 1 = Real 0.0
   SSAVelocity 2 = Real 0.0

! Thickness
  H = Real 1000.0

! free surface elevation
  Zs = Variable Coordinate 1
    Real MATC "-tx*tan(Slope)"
! as Zs is computed on the bottom surface it has to be exported on the bulk if we want to use it to update the 3D mesh
  Bottom Zs = Variable Coordinate 1
      Real MATC "-tx*tan(Slope)"

! bedrock elevation
  Zb = Variable Coordinate 1
    Real MATC "-tx*tan(Slope)-1000.0"

  Viscosity = REAL $eta
  Density = REAL $rhoi
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
  Flow BodyForce 1 = Real 0.0                          
  Flow BodyForce 2 = Real 0.0
  Flow BodyForce 3 = Real $gravity 

  Top Surface Accumulation = Real 0.0
  Bottom Surface Accumulation = Real 0.0


  Zs = Variable Zb , H
     REAL MATC "tx(0)+tx(1)"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1

! Material properties
  Viscosity Exponent = Real $1.0/n                
  Critical Shear Rate = Real 1.0e-10

  SSA Mean Viscosity = Variable "int Viscosity", thickness
       REAL MATC "tx(0)/tx(1)"
  SSA Mean Density = Variable "int Density", thickness
       REAL MATC "tx(0)/tx(1)"

! slip coeffs
  SSA Friction Law = String "linear"
  SSA Friction Parameter = Variable coordinate 1 , Coordinate 2
     Real  MATC "1.0e-3*(1.0 + sin(2.0*pi* tx(0) / L)*sin(2.0*pi* tx(1) / L))  
End
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
  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

End

Solver 2
!  StructuredProjectToPlane: used to compute integrated viscosity and density
!   thickness will also be used to compute mean viscosity and density
  Equation = "HeightDepth"
  Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"
  Active Coordinate = Integer 3

  Operator 1 = depth
  Operator 2 = height
  Operator 3 = thickness

  !! compute the integrated horizontal Viscosity and Density
  Variable 4 = Viscosity
  Operator 4 = int

  Variable 5 = Density
  Operator 5 = int

End

Solver 3
  Equation = "SSA"
  Variable = -dofs 2 "SSAVelocity"

  Procedure = "ElmerIceSolvers" "SSABasalSolver"

  Linear System Solver = Direct         
  Linear System Direct Method = umfpack

  Nonlinear System Max Iterations = 100
  Nonlinear System Convergence Tolerance  = 1.0e-12
  Nonlinear System Newton After Iterations = 5
  Nonlinear System Newton After Tolerance = 1.0e-05
  Nonlinear System Relaxation Factor = 1.00

End

!! Thickness solver
Solver 4
   Equation = "Thickness"
   Variable = -dofs 1 "H"

   Exported Variable 1 = -dofs 1 "H Residual"

   Exported Variable 2 = -dofs 1 "dHdt"

   Compute dHdT = Logical True

  Procedure = "ElmerIceSolvers" "ThicknessSolver"
   Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"

   Linear System Solver = Iterative
   Linear System Max Iterations = 1500
   Linear System Iterative Method = BiCGStab
   Linear System Preconditioning = ILU0
   Linear System Convergence Tolerance = Real 1.0e-12
   Linear System Abort Not Converged = False
   Linear System Residual Output = 1500

! equation is linear if no min/max
   Nonlinear System Max Iterations = 50
   Nonlinear System Convergence Tolerance  = 1.0e-6
   Nonlinear System Relaxation Factor = 1.00


  Stabilization Method = stabilized
  Apply Dirichlet = Logical False

!!  the convection velocity (mean horizontal velocity)
     Flow Solution Name = String "SSAVelocity"
End

!! This solver is just used to update Zs (look for the value of Zs in Body forces
Solver 5
   Equation = "UpdateExport"
   Procedure = "ElmerIceSolvers" "UpdateExport"
   Variable = -nooutput "dumy"

   Exported Variable 1 = -dofs 1 "Zs"
End


Solver 6
!  StructuredProjectToPlane: used to export Zs that has been computed on the bottom surf
  Equation = "HeightDepth"
  Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"
  Active Coordinate = Integer 3


 !! compute the integrated horizontal Viscosity and Density
  Variable 1 = Zs
  Operator 1 = bottom

End

Solver 7
  Exec Solver = never
  Equation = "Navier-Stokes"

   Stabilization Method = String Bubbles
   Flow Model = Stokes

   Exported Variable 1 = -dofs 1 "Viscosity"
   Exported Variable 2 = -dofs 1 "int Viscosity"
   Exported Variable 3 = -dofs 1 "Density"
   Exported Variable 4 = -dofs 1 "int Density"
   Exported Variable 5 = -dofs 1 "thickness"
   Exported Variable 6 = -dofs 1 "bottom Zs"

   Exported Variable 7 = -dofs 1 "dSdt"
   Exported Variable 8 = -dofs 1 "dS"

   Exported Variable 9 = -dofs 1 "Zb"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
   Active Solvers(4) = 1 2 6 7
End

Equation 2
  Active Solvers(3) = 3 4 5 
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Target 1 to 4 are the preserved 2D baseline
Boundary Condition 1
  Target Boundaries = 1

End

! Periodic Right
Boundary Condition 2
  Target Boundaries = 2
  Periodic BC = 4 
  Periodic BC Translate(2) = Real $L 0.0  
  Periodic BC SSAVelocity 1  = Logical True
  Periodic BC SSAVelocity 2  = Logical True
  Periodic BC H  = Logical True
  
End

Boundary Condition 3
  Target Boundaries = 3

  Periodic BC = 1 
  Periodic BC Translate(2) = Real 0.0 $L 
  Periodic BC SSAVelocity 1  = Logical True
  Periodic BC SSAVelocity 2  = Logical True
  Periodic BC H  = Logical True

End

Boundary Condition 4
  Target Boundaries = 4
End

!! Target 5 to 8 are the extruded BCs
Boundary Condition 5
  Target Boundaries = 5
End
Boundary Condition 6
  Target Boundaries = 6
End
Boundary Condition 7
  Target Boundaries = 7
End
Boundary Condition 8
  Target Boundaries = 8
End

!! bottom and surf come next
Boundary Condition 9
  Body Id = 2
  Bottom Surface = Equals Zb
End

Boundary Condition 10
  Top surface = Equals "bottom Zs"
End