$restartfile ="ASF_ti_1.sif"
$resultfile ="ASF_ti_2.sif"
$outputfile ="ASF_ti_2.sif"
$Lambda = 10^6
check keywords warn
!echo on
Header
Mesh DB "../." "ASFic3_B3_v2_2"
End
$gravity = 9.7696e15
$rhoi = 9.1501e-19 ! replaces: MATC "910.0*1.0E-06*(31556926.0)^(-2.0)" or 9.1376e-19
Constants
! needed for rheology
Gas Constant = Real 8.314 !Joule/mol x K
End
!! water pressure at glacier front
$ function waterpressure(Z) {\
rhow = 1012.0;\
waterline = 0.0;\
G = 9.81;\
_waterpressure = 0.0;\
if (Z>waterline) {\
_waterpressure = 0.0;\
}else {\
_waterpressure = 1.0 * rhow * G * (waterline - Z);\
}\
}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Simulation !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Coordinate System = Cartesian 3D
Simulation Type = Steady State
Steady State Min Iterations = 1
Steady State Max Iterations = 50!25
Output Intervals = 1
Output File = $resultfile".result"
Post File = $outputfile".ep"
Restart File = $restartfile".result"
Restart Position = 0
Initialize Dirichlet Conditions = Logical False
max output level = 9
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! BODY !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body 1
Name = "glacier1"
Equation = 1
Body Force = 1
Material = 1
Initial Condition = 1
End
Body 2
Name = "Bed"
Equation = 2
Initial Condition = 1
Material = 1
End
Body 3
Name = "Surf"
Equation = 3
Initial Condition = 1
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! INIT !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
Pressure = Real 0.0
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
Velocity 3 = Real 0.0
VeloD 1 = Real 0.0
VeloD 2 = Real 0.0
VeloD 3 = Real 0.0
VeloD 4 = Real 0.0
VsurfIni 1 = Real -999
VsurfIni 2 = Real -999
Depth = Real 0.0
DJDB = Real 0.
CostValue = Real 0.
Beta = Real -2.0
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! BODY FORCE !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Flow BodyForce 1 = Real 0.0
Flow BodyForce 2 = Real 0.0
Flow BodyForce 3 = Real -$gravity
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! MATERIAL !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
Density = Real $rhoi
Critical Shear Rate = Real 1.0e-10
Sliding = equals beta
! we want to have the Cauchy stress (for ComputeDevStress solver)
!----------------------------------
Cauchy = Logical True
!mesh update solver stuff
!------------------------
Mesh Elastic Modulus = 1.0
Mesh Poisson Ratio = 0.3
! Glen's flow law (using Glen)
!----------------
! viscosity stuff
!----------------
Viscosity Model = String "Glen"
! Viscosity has to be set to a dummy value
! to avoid warning output from Elmer
Viscosity = Real 1.0
Glen Exponent = Real 3.0
! Rate factors (Paterson value in MPa^-3a^-1)
Rate Factor 1 = Real 1.258e13
Rate Factor 2 = Real 6.046e28
! these are in SI units - no problem, as long as
! the gas constant also is
Activation Energy 1 = Real 60e3
Activation Energy 2 = Real 139e3
Glen Enhancement Factor = Real 1.0
! the variable taken to evaluate the Arrhenius law
! in general this should be the temperature relative
! to pressure melting point. The suggestion below plugs
! in the correct value obtained with TemperateIceSolver
Temperature Field Variable = String "Temperature Homologous"
! the temperature to switch between the
! two regimes in the flow law
Limit Temperature = Real -10.0
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! SOLVER !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Solver 1
Equation = "Mesh Update"
Exec Solver = "Before Simulation"
Linear System Solver = Iterative
Linear System Max Iterations = 500
!Linear System Direct Method = Mumps
Linear System Iterative Method = BiCGStab
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0e-6
Steady State Convergence Tolerance = 1.0e-04
Linear System Residual Output = 1
Optimize Bandwidth = Logical False
End
Solver 2
Equation = "Navier-Stokes"
Stabilization Method = Stabilized ! Alternative: Stabilization Method = String Bubbles
flow model = Stokes
Linear System Solver = Direct
Linear System Direct Method = mumps
Nonlinear System Max Iterations = 100
Nonlinear System Convergence Tolerance = 1.0e-5
Nonlinear System Newton After Iterations = 200
Nonlinear System Newton After Tolerance = 1.0e-7
Nonlinear System Relaxation Factor = 1.00
Steady State Convergence Tolerance = Real 1.0e-6
Exported Variable 1 = String "Beta"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "DJDB"
Exported Variable 2 DOFs = 1
Exported Variable 3 = String "CostValue"
Exported Variable 3 DOFs = 1
Exported Variable 4 = String "Temperature"
Exported Variable 4 DOFs = 1
End
Solver 3
Equation = "NS-Dirichlet"
Variable = VeloD
Variable Dofs = 4
procedure = "FlowSolve" "FlowSolver"
Stabilization Method = Stabilized
flow model = Stokes
Linear System Solver = Direct
Linear System Direct Method = mumps
Nonlinear System Max Iterations = 100
Nonlinear System Convergence Tolerance = 1.0e-5
Nonlinear System Newton After Iterations = 200
Nonlinear System Newton After Tolerance = 1.0e-7
Nonlinear System Relaxation Factor = 1.00
Steady State Convergence Tolerance = Real 1.0e-6
End
!!!! Compute Stresses from Neumann Problem
Solver 4
Equation = String "StressSolverN"
! Procedure = File "ComputeDevStressNS" "ComputeDevStress"
Procedure = File "ElmerIceSolvers" "ComputeDevStress"
Variable = -nooutput "SijN"
Variable DOFs = 1
Flow Solver Name = String "Flow Solution"
Stress Variable Name = String "StressN"
Exported Variable 1 = -nooutput "StressN" ! [Sxx, Syy, Szz, Sxy] in 2D
! [Sxx, Syy, Szz, Sxy, Syz, Szx] in 3D
Exported Variable 1 DOFs = 6 ! 4 in 2D, 6 in 3D
Linear System Solver = Direct
Linear System Direct Method = mumps
End
!!!! Compute Stresses from Dirichlet Problem
Solver 5
Equation = String "StressSolverD"
! Procedure = File "ComputeDevStressNS" "ComputeDevStress"
Procedure = File "ElmerIceSolvers" "ComputeDevStress"
Variable = -nooutput "SijD"
Variable DOFs = 1
Flow Solver Name = String "VeloD"
Stress Variable Name = String "StressD"
Exported Variable 1 = -nooutput "StressD" ! [Sxx, Syy, Szz, Sxy] in 2D
! [Sxx, Syy, Szz, Sxy, Syz, Szx] in 3D
Exported Variable 1 DOFs = 6 ! 4 in 2D, 6 in 3D
Linear System Solver = Direct
Linear System Direct Method = mumps
End
!!!! Compute Cost
Solver 6
Equation = "Cost"
Variable = -nooutput "CostV"
Variable DOFs = 1
! procedure = "CostSolver_ArthernGud" "CostSolver"
procedure = "ElmerIceSolvers" "CostSolver_Robin"
Cost Variable Name = String "CostValue" ! Name of Cost Variable
Neumann Solution Name = String "Flow Solution"
Neumann Stress Name = String "StressN"
Dirichlet Solution Name = String "VeloD"
Dirichlet Stress Name = String "StressD"
Optimized Variable Name = String "Beta"
Lambda = Real $Lambda ! Regularization Coef
Cost Filename = File "CostRobin_"$outputfile".dat"
end
!!!!! Compute Derivative of Cost function / Beta
Solver 7
! Exec Solver = Never
Equation = "DJDBeta"
Variable = -nooutput "DJDBV"
Variable DOFs = 1
! procedure = "DJDBeta_ArthernGud" "DJDBeta"
procedure = "ElmerIceSolvers" "DJDBeta_Robin"
Neumann Solution Name = String "Flow Solution"
Dirichlet Solution Name = String "VeloD"
Optimized Variable Name = String "Beta" ! Name of Beta variable
Gradient Variable Name = String "DJDB" ! Name of gradient variable
Lambda = Real $Lambda ! Regularization Coef
PowerFormulation = Logical True ! Slip coef define as 10^Beta
Beta2Formulation = Logical False ! Slip coef not define as Beta^2
end
!!!!! Optimization procedure
Solver 8
Equation = "UpdateBeta_m1qn3"
Variable = -nooutput "UBV"
Variable DOFs = 1
! procedure = "UpBeta_m1qn3Parallel" "UpdateBeta"
procedure = "ElmerIceSolvers" "Optimize_m1qn3Parallel"
Cost Variable Name = String "CostValue"
Optimized Variable Name = String "Beta"
Gradient Variable Name = String "DJDB"
Min Value = Real -20.0 ! Min Value for Beta (better to leave large values)
Max Value = Real +20.0 ! Max Value for Beta
M1QN3 dxmin = Real 1.0e-10
M1QN3 epsg = Real 1.e-6
M1QN3 niter = Integer 200
M1QN3 nsim = Integer 200
M1QN3 impres = Integer 5
M1QN3 DIS Mode = Logical False
M1QN3 df1 = Real 0.2
M1QN3 normtype = String "dfn"
M1QN3 OutputFile = File "M1QN3Robin_"$outputfile".out"
M1QN3 ndz = Integer 10
end
Solver 9
Exec Solver = Before Simulation
Equation = "Pointwise Data"
Variable = -nooutput "Dummy"
Variable DOFs = 1
! Procedure = "pointwise" "InterpolatePointValue"
Procedure = "../../Solver/ElmerSolvers/pointwise" "InterpolatePointValue"
Nonlinear System Max Iterations = 1
Variable 1 = String "VsurfIni 1"
Variable Data 1 = String "v95s_p11s.x"
Variable 1 Supporting Points = Integer 3
Variable 1 Dimensions = Integer 2
Variable 1 Directions(2) = Integer 1 2
Exported Variable 1 = VsurfIni 1
Exported Variable 1 DOFS = Integer 1
Variable 2 = String "VsurfIni 2"
Variable Data 2 = String "v95s_p11s.y"
Variable 2 Supporting Points = Integer 3
Variable 2 Dimensions = Integer 2
Variable 2 Directions(2) = Integer 1 2
Exported Variable 2 = VsurfIni 2
Exported Variable 2 DOFS = Integer 1
End
Solver 10
Equation = "Flowdepth"
Exec Solver = "Before Simulation"
Procedure = File "ElmerIceSolvers" "FlowDepthSolver"
Variable = String "Depth"
Variable DOFs = 1
Linear System Solver = "Direct"
Linear System Direct Method = "MUMPS"
! Linear System Iterative Method = "BiCGStab"
! Linear System Max Iterations = 300
! Linear System Convergence Tolerance = 1.0E-09
! Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 1
! this sets the direction
! -1 is negative z-direction (upside down)
! +1 is positive (downside up)
Gradient = Real -1.0E00
! switch that to True, if you want to have
! free surface gradients to be computed
!------------------------------------
Calc Free Surface = Logical True
! the name for the exported (if not existing) added variable
! the gradients will be stored in variables with the base
! name given and "Grad1" and (in 3 dimensions) "Grad2" added,
! so in our case "FreeSurfGrad1" and "FreeSurfGrad2"
! again, if those variables did not exist, they will be
! automatically created
!-----------------------------------------------------------
Freesurf Name = String "FreeSurf"
End
Solver 11
Exec Solver = After Timestep
Equation = "SaveLine"
Procedure = File "SaveData" "SaveLine"
Filename = $outputfile"I.dat"
File Append = Logical False
End
Solver 12
Equation = "SaveMaterials"
Exec Solver = After TimeStep
Procedure = File "SaveData" "SaveMaterials"
Parameter 1 = String "Viscosity"
Parameter 2 = String "Temperature"
End
Solver 13
Exec Solver = String "after timestep"
exec interval = 1
Equation = String "ResultOutput"
Procedure = File "ResultOutputSolve" "ResultOutputSolver"
Output File Name = String "pv."$outputfile"."
Output Format = String "vtu"
Vtu Format = Logical True
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! EQUATION !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Equation 1
!ice
Active Solvers(10) = 1 2 3 4 5 6 10 11 12 13
NS Convect = Logical False
End
Equation 2
!bed
Active Solvers(2) = 7 8
NS Convect = Logical False
End
Equation 3
!surf
Active Solvers(1) = 9
NS Convect = Logical False
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! BOUNDARY !!
!! 1 bed, 2 surf, 3 border !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Bedrock
Boundary Condition 1
Name = "bed"
Target Boundaries(1) = 1
Body Id = 2
Save Line = Logical True
Flow Force BC = Logical True
Normal-Tangential Velocity = Logical True
Normal-Tangential VeloD = Logical True
Velocity 1 = Real 0.0e0
VeloD 1 = Real 0.0e0
Slip Coefficient 2 = Variable Beta
Real MATC "10^tx(0)"
Slip Coefficient 3 = Variable Beta
Real MATC "10^tx(0)"
End
! Upper Surface
Boundary Condition 2
Name = "surface"
Target Boundaries = 2
Save Line = Logical True
Body Id = 3
Depth = Real 0.0
VeloD 1 = Equals VsurfIni 1
VeloD 2 = Equals VsurfIni 2
End
!! the artificial side wall:
Boundary Condition 3
Name = "sides"
Target Boundaries(2) = 3 4
External Pressure = Variable Coordinate 3
Real MATC "-1.0*waterpressure(tx)*1.0E-06"
End