mlb_tmc_linsys.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 run
!!!!!!!!!!!!!!!!!!!!!!!!!!
$name="S_1995_DEMs_tmc_loads"
#MINH=20.0
#Tc = -1.0
#directmethod="mumps" 
!---LUA BEGIN
! assert(loadfile('./icematerial.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_init_2d.result"






! 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 = "Steady State"             

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

  ! Coupled iterations between different solvers
  !---------------------------------------------
  Steady State Max Iterations = 20
  Steady State Min Iterations = 1

  Output Intervals = 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 happens in Solver 1
  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 = 5
End

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

Solver 1
  Exec Solver = before all
  Equation = "InterpSolver"

! Perform the interpolation solver before all other solvers, only once.
  Exec Solver = Before all
  Procedure = "Mesh2MeshSolver" "Mesh2MeshSolver"

  ! Restart is always in a serial mesh
  Mesh = -single $restartdir
  Restart File = $restartfile

  ! We use the primary mesh an single process and unextruded
  ! We want to have unique copy of it.   
  Mesh Enforce Local Copy = Logical True

  ! These are the variables for restart 
  Restart Position = Integer 0
  Restart Time = Real 2000.0
  Restart Variable 1 = String "bedrockDEM"
  Restart Variable 2 = String "surfaceDEM1995"

! Continue even if some node is not found, deal with problems!
  Restart Error Continue = Logical True  

! These are the variables to be mapped from 2D ThisMesh to 3D TargetMesh'
! All existing variables are automatically mapped. 
!  Variable 1 = "bedrockDEM"
!  Variable 2 = "surfaceDEM1995"

! Ensures that we perform interpolation on plane 
  Interpolation Passive Coordinate = Integer 3 
End 


! Maps the constant-thickness mesh between given bedrock and surface topology
Solver 2 
  Exec Solver = "before simulation"
  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"

  Exported Variable 1 = "bedrockDEM"
  !Exported Variable 1 Mask = String "BedRock"

  Exported Variable 2 = "surfaceDEM1995"
  !Exported Variable 2 Mask = String "Surface"
End


! Computes height and depth assuming an extruded mesh.
Solver 3 
  Exec Solver = "before simulation" 
  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
  Flow Model = Stokes
  Stabilization Method = String Stabilized

  !Non-linear iteration settings:
  !------------------------------ 
  Nonlinear System Max Iterations = 1
  Nonlinear System Convergence Tolerance  = 1.0e-5
  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

  include linsys/gcr.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
  Calculate Loads = True
  !Calculate Weights = True
  Calculate Boundary Weights = True

  Exported Variable 1 = -dofs 4 Flow Solution Loads
!  Exported Variable 2 = Flow Solution Boundary Weights
End

!-----------------------------------------------------
! heat transfer limited by the pressure melting point
! as upper limit
!-----------------------------------------------------
Solver 5
 !Exec Solver = "Never"
  Equation = String "Homologous Temperature Equation"
  Procedure =  File "HeatSolveVec" "HeatSolver"
  ! Comment next line in parallel, as EliminateDirichlet does
  ! not work in parallel
  !------------------------------------------------------------
!  Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
  Variable = String "Temperature"
  Variable DOFs = 1
  Stabilization method = Stabilized
  Stabilize = True
  Optimize Bandwidth = Logical True
 ! Linear System Solver = "Direct"
 ! Linear System Direct Method = #directsolver
  Linear System Solver = "Iterative"
  Linear System Iterative Method = "GCR"
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0E-08
  Linear System Abort Not Converged = True
  Linear System Preconditioning = "ILU1"
  Linear System Residual Output = 1
  Nonlinear System Convergence Tolerance = 1.0E-04
  Nonlinear System Max Iterations = 50
  !Nonlinear System Relaxation Factor = Real 9.999E-01
  Steady State Convergence Tolerance = 1.0E-03
  !Steady State Relaxation Factor = Real #1.0/3.0
  ! the contact algorithm (aka Dirichlet algorithm)
  !-----------------------------------------------------
  Apply Limiter = Logical True
  Calculate Loads = Logical True
    Save Limiter = Logical True
  Calculate Boundary Weights = Logical True
  ! those two variables are needed in order to store
  ! the relative or homologous temperature as well
  ! as the residual
  !-------------------------------------------------
  Exported Variable 1 = -dofs 1 "Relative Temperature"
!  Exported Variable 2 = -dofs 1 "Temperature Loads"
!  Exported Variable 3 = -dofs 1 "Temperature Boundary Weights"
  Exported Variable 2 = -dofs 1 "Friction Loads"
  Nonlinear Update Exported Variables = Logical True
  Update Exported Variables = Logical True
End





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

Body 1
  Name = "glacier"
  Equation = 1
  Material = 1
  Body Force = 1
  Initial Condition = 1
End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! EQUATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
  Active Solvers(5) = 1 2 3 4 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
  Temperature = Variable Coordinate 3 ! lapserate of -10/1000 C/m
      Real
            0.0   272.14
	  100.0   272.15
         1000.0   263.15
      End 
End


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODY FORCE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
  Flow BodyForce 1 = 0.0
  Flow BodyForce 2 = 0.0 
  Flow BodyForce 3 = #gravity
  Temperature Upper Limit = Real 273.15
  Temperature = Real 273.15
  Temperature Condition =Variable Coordinate 3 
      Real
            0.0      -1.0
	  399.9999   -1.0
	  400.        1.0
         1000.0       1.0
      End
  Relative Temperature = Variable Temperature, Pressure
     Real Procedure "ElmerIceUSF" "RelativeTemperature" 
  !Variable Coordinate 3
  !  Real lua "settemperate(tx[0])"
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 = Variable Temperature
    Real lua "tx[0] - 273.15"
  ! Heat transfer stuff
  Heat Capacity = Variable Temperature
    Real lua "capacity(tx[0])*yearinsec^(2.0)"
  !  Real lua "capacity(tx[0])*yearinsec^(2.0)"

  Heat Conductivity = Variable Temperature
    Real lua "conductivity(tx[0])*yearinsec*Pa2MPa"
    
  Pressure Scaling Factor = Real #1.0/Pa2MPa
  Temperature Upper Limit = Real 273.15 ! we ignore pressure melting point
     
End

Boundary Condition 1
  ComputeNormal = Logical False
  Name = "sides1"
  Target Boundaries(1) = 1
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 0
End

Boundary Condition 2
  ComputeNormal = Logical False
  Name = "sides2"
  Target Boundaries(1)  = 2
  ! no slip
  !-------------------------
  Velocity 1 = 0
  Velocity 2 = 0
  Velocity 3 = 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

  ! slip conditions
  Normal-Tangential Velocity = True
  Velocity 1 = Real 0.0
  Weertman Friction Coefficient = Real 0.1
  Weertman Exponent = Real $1.0/3.0
  Weertman Linear Velocity = Real 0.00001
  Friction Newton Epsilon = Real 1.0E-03
  Friction Newton Linearization = Logical True


  Bottom Surface = Equals "bedrockDEM"

! Mask for creating the restart fields only where needed
! Refererred by: Exported Variable 1 Mask = ...
  Bedrock = Logical True
  Temperature Flux BC = Logical True
  Heat Flux = Real #0.050 * yearinsec * Pa2MPa
  ! 50 mW m^-2
  Friction Loads = Variable Flow Solution Loads 1, Flow Solution Loads 2, Flow Solution Loads 3, Velocity 1, Velocity 2, Velocity 3
    Real Lua "frictionloads(tx[0],tx[1],tx[2],tx[3],tx[4],tx[5])"
  Target Boundaries(1) = 3
  Temperature Load = Variable Flow Solution Loads 1, Flow Solution Loads 2, Flow Solution Loads 3, Velocity 1, Velocity 2, Velocity 3
    Real Lua "frictionloads(tx[0],tx[1],tx[2],tx[3],tx[4],tx[5])"
  Calculate Temperature Boundary Weights = Logical True
End

Boundary Condition 4
  ComputeNormal = Logical False
  Name = "surface"
  Top Surface = Equals "surfaceDEM1995"

! Referred by: Exported Variable 2 Mask = ...
  Surface = Logical True
Target Boundaries(1) = 4
  Temperature = Variable Coordinate 3 ! lapserate of -10/1000 C/m
      Real
            0.0   272.15
	  100.0   272.14
         1000.0   263.15
      End 
End