ismip_coulomb_pw.sif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                       !!
!! Test for Coulomb friction law         !!
!! N is computed from a given            !! 
!! water pressure and the stress         !!
!!                                       !!
!! From test ISMIP-HOM B020              !!
!!                                       !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

check keywords warn
echo on
$L = 20.0e3
$Slope = 0.5 * pi / 180.0

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

Header
  Mesh DB "." "rectangle"
End

Constants
! No constant Needed
End

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

! More than one steady iteration is needed as the stress are 
! not known at the first iteration 
  Steady State Min Iterations = 2
  Steady State Max Iterations = 10

  Post File = "ismip_coulomb_pw.vtu"
  max output level = 3   
End

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
  Pressure = Real 0.0 
  Velocity 1 = Real 0.0 
  Velocity 2 = Real 0.0 

  Stress 1 = Real 0.0
  ! initialize with hydrostatic stress (<0)
  Stress 2 = Variable Coordinate 1
    Real MATC "rhoi*gravity*(1000.0-500.0*sin(2.0*pi*tx/L))"
  Stress 3 = Real 0.0
  Stress 4 = Real 0.0

End

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
  Density = Real $rhoi 
  Viscosity Model = String "power law"
  Viscosity = Real $eta                       
  Viscosity Exponent = Real $1.0/n                
  Critical Shear Rate = Real 1.0e-10

  Cauchy = Logical True
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1  
  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
End

Solver 2
  Equation = "Normal vector"
  Procedure = "ElmerIceSolvers" "ComputeNormalSolver"
  Variable = -dofs 2 "Normal Vector"   
  Optimize Bandwidth = Logical False 
  ComputeAll = Logical False
End

Solver 3
  Equation = "Navier-Stokes"
  Exported Variable 1 = -dofs 1 "dSdt" 
  Exported Variable 2 = -dofs 1 "dS" 
  
  Stabilization Method = String Bubbles
  Flow Model = Stokes

  Linear System Solver = Direct         
  Linear System Direct Method = umfpack

  Nonlinear System Max Iterations = 100
  Nonlinear System Convergence Tolerance  = 1.0e-5
  Nonlinear System Newton After Iterations = 5
  Nonlinear System Newton After Tolerance = 1.0e-02
  Nonlinear System Relaxation Factor = 1.00
  Nonlinear System Reset Newton = Logical True

  Steady State Convergence Tolerance = Real 1.0e-3
End

Solver 4
  Equation = String "StressSolver"
  Procedure =  File "ElmerIceSolvers" "ComputeDevStress"
  Variable = -dofs 1 -nooutput "Sij"
  Flow Solver Name = String "Flow Solution"
  Stress Variable Name = String "Stress"
  Exported Variable 1 = -dofs 4 "Stress" ! [Sxx, Syy, Szz, Sxy] in 2D
  Linear System Solver = Direct         
  Linear System Direct Method = umfpack
End
									   


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Bedrock
Boundary Condition 1
  Target Boundaries = 1
  Normal-Tangential Velocity = Logical True
 !! Water pressure given through the Stokes 'External Pressure' parameter 
 !! What should be given here is -pw (negative)
 External Pressure = Variable Coordinate 1
    Real MATC "0.3+rhoi*gravity*(1000.0-500.0*sin(2.0*pi*tx/L))"

  Velocity 1 = Real 0.0
	     
  Slip Coefficient 2 =  Variable Coordinate 1
    Real Procedure "ElmerIceUSF" "Friction_Coulomb"
  !! Parameters needed for the Coulomb Friction Law
  Friction Law Sliding Coefficient = Real $As
  Friction Law Post-Peak Exponent  = Real 1.0      
  Friction Law Maximum Value = Real 0.5           
  Friction Law PowerLaw Exponent = Real $n         
  Friction Law Linear Velocity = Real 1.0e-4
    
  ComputeNormal = Logical True

  Bottom Surface = Variable Coordinate 1
    Real MATC "-tx*tan(Slope)-1000.0+500.0*sin(2.0*pi*tx/L)"
End

! Periodic Right
Boundary Condition 2
  Target Boundaries = 2
  Periodic BC = 4 
  Periodic BC Translate(2) = Real $L 0.0  
  Periodic BC Velocity 1  = Logical True
  Periodic BC Velocity 2  = Logical True
  Periodic BC Pressure  = Logical True
  
  ComputeNormal = Logical False
End

! Upper Surface
Boundary Condition 3
  Target Boundaries = 3
  ComputeNormal = Logical False

  Top Surface = Variable Coordinate 1
    Real MATC "-tx*tan(Slope)"
End

! Periodic Left
Boundary Condition 4
  Target Boundaries = 4
  ComputeNormal = Logical False
End

Solver 3 :: Reference Norm = Real 63.28803