mlb_tmc_trans_linsys_smooth.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_vec_smooth"

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







! 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) = 50      !#Intervals
  Timestep Sizes(1) = $1.0/365.25

  ! 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 = "s_1995_dems_tmc_vec.result"
  Restart Position = 0
  Restart Before Initial Conditions = Logical True
  Interpolation Passive Coordinate = Integer 3 
  Restart Variable 1 = String "bedrockDEM"
  Restart Variable 2 = String "surfaceDEM1995"
  Restart Variable 3 = String "Velocity 1"
  Restart Variable 4 = String "Velocity 2"
  Restart Variable 5 = String "Velocity 3"
  Restart Variable 6 = String "Pressure"
  Restart Variable 7 = String "Temp"
  Restart Variable 8 = String "Temp Homologous"
  ! how verbose the solver should be
  !  3 = Only warnings
  ! 32 = Maximum verbosity
  !-------------------------------------------------------
  Max Output Level = 5
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 #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"
  !Mesh Velocity Variable = String "mv"
  Mesh Velocity Variable = String "Mesh Velocity 3"
  Mesh Velocity First Zero  = True
  Mesh Update Variable = String "dz"
  !Mesh Update Variable = String "Mesh Update 3"


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


  !Exported Variable 3 = -dofs 1 "dz"
 
End

Solver 2
  Equation = "Mesh smoother"
  Exec Solver = "Before Timestep"
  !Exec Solver = "Never"
  Procedure = File "MeshSolve" "MeshSolver"
  Variable = -dofs 3 "Mesh Update"
!  Linear System Solver = Iterative
!  Linear System Iterative Method = BiCGStab

  Linear System Solver = Direct
  Linear System Direct Method = MUMPS
  

  Linear System Preconditioning = ILU1
  Linear System Max Iterations = 1000
  Linear System Residual Output = 20
  Linear System Convergence Tolerance = 1.0e-8
  Steady State Convergence Tolerance = 1.0e-4
End

! Computes height and depth assuming an extruded mesh.
Solver 3 
  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 4
! Exec Solver = "Never"
  Equation = "Stokes-Vec"
  Procedure = "IncompressibleNSVec" "IncompressibleNSSolver"
  Div-Curl Discretization = Logical False

  Optimize Bandwidth = Logical True
  Stokes Flow = Logical True
  Stabilize = true
  Flow Model = Stokes

  !Non-linear iteration settings:
  !------------------------------ 
  Nonlinear System Max Iterations = 10
  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 = 10
 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 False

! 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 5
  !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"

  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 = 250
  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
  Exported Variable 3 = -dofs 1 "dz"
  !Exported Variable 4 = -dofs 1 "mv"
  !Exported Variable 4 = -dofs 3 "Mesh Velocity"
  !Exported Variable 4 = -dofs 3 "Mesh Update"
  Apply Dirichlet = True
End

Solver 6
   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 "Ref Zs"
   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(5) = 1 2 3 4 5
  Convection = Computed
  Flow Solution Name = String "Flow Solution"
End

Equation 2
  Active Solvers(1) = 6
  Convection = Computed
  Flow Solution Name = String "Flow Solution"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! INITIAL CONDITIONS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
  ! Initial velocity is not really needed when using newtonian start
   Mesh Velocity 1 = Real 0.0
   Mesh Velocity 2 = Real 0.0
   Mesh Velocity 3 = Real 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])"
  Ref 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

  Temp = Real 273.15
  Temp Condition =Variable Coordinate 3
      Real
            0.0      -1.0
          399.9999   -1.0
          400.        1.0
         1000.0       1.0
      End
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
End
Body Force 2
 ! Zs Accumulation = Real 0.0
 Zs Accumulation = Variable Coordinate 3, glaciated
    Real lua "accum(tx[0],tx[1])"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MATERIAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
  Name = "Ice"
  Density = Real #rhoi*(yearinsec^(-2.0))*Pa2MPa

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

  ! Nonnewtonian viscosity
  Viscosity Model = String Glen
  Glen Exponent = Real 3.0
  Critical Shear Rate =  Real 1.0E-10
  ! 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*yearinsec^(2.0)
  Temp Heat Conductivity =  Real #2.2*yearinsec*Pa2MPa

  !Temp Heat Capacity = Variable Temp
  !  Real lua "capacity(tx[0])*yearinsec^(3.0)"
  !  Real lua "capacity(tx[0])*yearinsec^(2.0)"

  !Temp Heat Conductivity = Variable Temp
  !  Real lua "conductivity(tx[0])*yearinsec*Pa2MPa"
  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
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Zs = Equals "Ref Zs"
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
End

Boundary Condition 2
  ComputeNormal = Logical False
  Name = "sides2"
  Target Boundaries(1)  = 2
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
  Zs = Equals "Ref Zs"
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
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 * yearinsec * Pa2MPa ! 200 -> 1000 mW m^-2
  !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"
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
  Mesh Update 3 = Real 0.0
End

Boundary Condition 4
  ComputeNormal = Logical False
  Name = "surface"
  Body ID = 2
  Top Surface = Equals "zs"
  Mesh Update 3 = Variable zs, Ref Zs
    Real LUA "tx[0] - tx[1]"
  
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
! Top Surface = Equals  surfaceDEM1995 !if you want to test constant geometry


! 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