mlb_tmc_trans_linsys_SI.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 ru
!!!!!!!!!!!!!!!!!!!!!!!!!!
$name="S_1995_DEMs_tmc_trans_SI"

#MINH=20.0
#directmethod="mumps" 
!---LUA BEGIN
! assert(loadfile('./icematerial_SI.lua'))()
! assert(loadfile('./accumulation_SI.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_tmc_si"






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


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! HEADER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Partitioned mesh
Header
  Mesh DB "." "outline62_lc75"
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! CONSTANTS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Constants
End


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
  Coordinate System  = "Cartesian 3D"
  Simulation Type = "Transient"
  
  Timestepping Method = "bdf"
  BDF Order = 1
  Output Intervals(1) = 1         !#OutPutIntervals
  Timestep Intervals(1) = 30      !#Intervals
  Timestep Sizes(1) = #3600.0*24.0 !#dt

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

  ! Coupled iterations between different solvers
  !---------------------------------------------
  Steady State Max Iterations = 10
  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
  Restart File = "$restartfile$.result"  
  Restart Position = 0
  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 = 3
End

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




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

  Active Coordinate = Integer 3

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

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

  ! These are needed to host the variables that are read in
  Variable = "Mesh Update"
  Mesh Velocity Variable = String "Mesh Velocity 3"
  Mesh Update Variable = String "Mesh Update"
  Mesh Velocity First Zero  = True

  Exported Variable 1 = -dofs 3 "Mesh Velocity"
  Exported Variable 2 = -dofs 1 "Mesh Update"

  Exported Variable 3 = "bedrockDEM"
  !Exported Variable 1 Mask = String "BedRock"
  Exported Variable 4 = "surfaceDEM1995"
  !Exported Variable 2 Mask = String "Surface"

  !Exported Variable 3 = -dofs 3 "Mesh Velocity"
 
End


! Computes height and depth assuming an extruded mesh.
Solver 2 
  Exec Solver = "before Timestep" 
  Equation = "HeightDepth"
  Procedure = "StructuredProjectToPlane" "StructuredProjectToPlane"
  Active Coordinate = Integer 3
  Operator 1 = depth
  Operator 2 = height
End 


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

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

  !Non-linear iteration settings:
  !------------------------------ 
  Nonlinear System Max Iterations = 1
  Nonlinear System Convergence Tolerance  = 1.0e-3
  Nonlinear System Newton After Iterations = 10
  Nonlinear System Newton After Tolerance = 1.0e-1
  ! Nonlinear System Relaxation Factor = 1.00
  ! Nonlinear System Reset Newton = Logical True
 Linear System Residual Output = 0
  include linsys.sif

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

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

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

! Some timing info
  !Boundary Assembly Timing = Logical True
  !Bulk Assembly Timing = Logical True
  !Solver Timing = Logical True
  !Linear System Timing = Logical True
End

!-----------------------------------------------------
! heat transfer limited by the pressure melting point
! as upper limit
!-----------------------------------------------------
Solver 4
  !Exec Solver = "Never"
  Equation = String "Homologous Temperature Equation"
  Procedure =  File "ElmerIceSolvers" "TemperateIceSolver"
  ! Comment next line in parallel, as EliminateDirichlet does
  ! not work in parallel
  !------------------------------------------------------------
!  Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
  Variable = String "Temp"
  Variable DOFs = 1
  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 Max Iterations = 500
  Linear System Convergence Tolerance = 1.0E-07
  Linear System Abort Not Converged = False
  Linear System Preconditioning = "ILU1"
  Linear System Residual Output = 0
  Nonlinear System Convergence Tolerance = 1.0E-06
  Nonlinear System Max Iterations = 50
  !Nonlinear System Relaxation Factor = Real 9.999E-01
  Steady State Convergence Tolerance = 1.0E-03
  ! the contact algorithm (aka Dirichlet algorithm)
  !-----------------------------------------------------
  !Apply Dirichlet = Logical True
  ! those two variables are needed in order to store
  ! the relative or homologous temperature as well
  ! as the residual
  !-------------------------------------------------
  Exported Variable 1 = String "Temp Homologous"
  Exported Variable 1 DOFs = 1
  Exported Variable 2 = String "Temp Residual"
  Exported Variable 2 DOFs = 1
  Apply Dirichlet = True
End

Solver 5
   Exec Solver = "after timestep"
   !Exec Solver = "Never"
   Equation = "Free Surface"
   Variable = String "Zs"
   Variable DOFs =  1
   ! needed for evaluating the contact pressure
   Exported Variable 1 = -dofs 1 "Zs Residual"
   ! needed for storing the initial shape (needed for updates)
   !Exported Variable 2 = -dofs 1 "RefZs"
   Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
  ! This would take the contrained points out of solution
  ! Use in serial run, only
 !  Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
   Linear System Solver = Iterative
   Linear System Max Iterations = 1500
   Linear System Iterative Method = BiCGStab
   Linear System Preconditioning = ILU1
   Linear System Convergence Tolerance = Real 1.0e-8
   !Linear System Abort Not Converged = False
   Linear System Residual Output = 1
   Nonlinear System Min Iterations = 2
   Nonlinear System Max Iterations = 100
   Nonlinear System Convergence Tolerance  = 1.0e-7
   !Nonlinear System Relaxation Factor = 0.60
   Steady State Convergence Tolerance = 1.0e-03
   Stabilization Method = Stabilized
   !Stabilization Method = Bubbles
   ! Apply contact problem    
   Apply Dirichlet = Logical True

! How much the free surface is relaxed
!  Relaxation Factor = Real 0.90      
 
End





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODIES (i.e., domains to compute on)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Body 1
  Name = "glacier"
  Equation = 1
  Material = 1
  Body Force = 1
  Initial Condition = 1
End
Body 2
  Name = "surface"
  Equation = 2
  Material = 2
  Body Force = 2
  Initial Condition = 2
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! EQUATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
  Active Solvers(4) = 1 2 3 4
  !Convection = Constant
  Convection = Computed
  Flow Solution Name = String "Flow Solution"
End

Equation 2
  Active Solvers(1) = 5
  Convection = Computed
  Flow Solution Name = String "Flow Solution"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! INITIAL CONDITIONS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
  ! Initial velocity is not really needed when using newtonian start
  !Velocity 1 = 0.0
  !Velocity 2 = 0.0
  !Velocity 3 = 0.0
  Temp = Variable Coordinate 3 ! lapserate of -10/1000 C/m
      Real
            0.0   272.14
	  100.0   272.15
         1000.0   263.15
      End 
End

Initial Condition 2
  Zs = Equals  surfaceDEM1995
!Zs = Variable  surfaceDEM1995,bedrockDEM
  !  Real lua "initzs(tx[0],tx[1])"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODY FORCE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
  Flow BodyForce 1 = 0.0
  Flow BodyForce 2 = 0.0 
  Flow BodyForce 3 = #gravity
  Velocity 1 = 0.0
  Velocity 1 Condition = Opposes glaciated
  Velocity 2 = 0.0
  Velocity 2 Condition = Opposes glaciated
  Velocity 3 = 0.0
  Velocity 3 Condition = Opposes glaciated
  !Flow Solution Passive =  Opposes glaciated
  Temp Passive  =  Opposes glaciated
  !Temp = Variable Coordinate 3 ! lapserate of -10/1000 C/m
  !    Real
  !          0.0   272.15
!	  100.0   272.14
!         1000.0   263.15
!      End 
  Temp Condition = Opposes glaciated
  
End
Body Force 2
  Zs Accumulation = Variable Coordinate 3, glaciated
     Real lua "accum(tx[0],tx[1])"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MATERIAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
  Name = "Ice"
  Density = Real #rhoi

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

  ! 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 = Equals Temp Homologous
  !Variable Temp
  !  Real lua "reltemp(tx[0])"

  ! heat transfer
    ! Heat transfer stuff
  !Temp Heat Capacity = Real #2100.0
  !Temp Heat Conductivity =  Real #2.2
  Temp Heat Capacity = Variable Temp
    Real lua "capacity(tx[0])"
  !  Real lua "capacity(tx[0])*yearinsec^(2.0)"

  Temp Heat Conductivity = Variable Temp
    Real lua "conductivity(tx[0])"
  Temp Upper Limit = Real 273.15 ! we ignore pressure melting point
  !Temp Lower Limit = Real 263.15
End

Material 2
   Name = "SMB"
  Min Zs = Variable "bedrockDEM"
    Real lua "tx[0] + MINH"
End

Boundary Condition 1
  ComputeNormal = Logical False
  Name = "sides1"
  Target Boundaries(1) = 1
  Normal-Tangential Velocity = True
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Zs = Variable bedrockDEM
    Real lua "tx[0] + MINH"
End

Boundary Condition 2
  ComputeNormal = Logical False
  Name = "sides2"
  Target Boundaries(1)  = 2
  Normal-Tangential Velocity = True
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Zs = Variable bedrockDEM
    Real lua "tx[0] + MINH"
End

!! DON'T CHANGE ORDER OF NEXT 2 BC's!
!! They are automaticaly created in internal extrusion
!! bedrock:
Boundary Condition 3
  ComputeNormal = Logical True
  Name = "bedrock"

  ! No-slip velocity conditions
  Velocity 1 = Real 0.0
  Velocity 2 = Real 0.0
  Velocity 3 = Real 0.0
  Bottom Surface = Equals "bedrockDEM"

! Mask for creating the restart fields only where needed
! Refererred by: Exported Variable 1 Mask = ...
  !Bedrock = Logical True
  Temp Flux BC = Logical True
  Temp Heat Flux = Real #0.050
  !Velocity 1 = 0.0
  !Velocity 1 Condition = Variable depth, height
  !  Real lua "(tx[0] + tx[1]) < 1.01*MINH"
  !Velocity 2 = 0.0
  !Velocity 2 Condition = Variable depth, height
  !  Real lua "(tx[0] + tx[1]) < 1.01*MINH"
  !Velocity 3 = 0.0
  !Velocity 3 Condition = Variable depth, height
  !  Real lua "(tx[0] + tx[1]) < 1.01*MINH"
End

Boundary Condition 4
  ComputeNormal = Logical False
  Name = "surface"
  Body ID = 2
  Top Surface = Equals "zs"


! Referred by: Exported Variable 2 Mask = ...
  Surface = Logical True

  Temp = Variable Coordinate 3 ! lapserate of -10/1000 C/m
      Real
            0.0   272.15
	  100.0   272.14
         1000.0   263.15
      End 
End