So I tried the different suggestions and that did help a bit. I am now able to converge but the result is clearly wrong. The computational time also increased quite a lot (which I was also expecting with the boundary condition). So I am now running Elmer in parallel with Hypre, but that should hopefully not make any difference on the results.
The new setup uses a Weertman sliding law with a constant coeffienent of 10.0 (i.e. no sliding). It is having a hard time resolving the velocity near the margin. I have attached a screenshot (top view) of Flow Solution 1. A similar pattern is seen in the other Flow Solutions.
I have tried to increase the area without ice (i.e. increase the distance from the ice filled area to the sides) but that didn't help. I have also tried to increase the amount of ice in the "ice-free" areas from 1 to 10 meters. Finally I have increased the number of vertical layers form 4 to 10 but also no luck.
I could not get the "normal" Navier-Stokes solver working. I get the following error (perhaps a compiler problem with the binary?),
Code: Select all
Loading user function library: [ElmerIceSolvers]...[IntegrateVertically_Init]
Loading user function library: [ElmerIceSolvers]...[IntegrateVertically]
OptimizeBandwidth: ---------------------------------------------------------
OptimizeBandwidth: Computing matrix structure for: integratevertically...done.
OptimizeBandwidth: Half bandwidth without optimization: 10208
OptimizeBandwidth:
OptimizeBandwidth: Bandwidth Optimization ...done.
OptimizeBandwidth: Half bandwidth after optimization: 645
OptimizeBandwidth: ---------------------------------------------------------
Loading user function library: [ElmerIceSolvers]...[ComputeNormalSolver_Init]
Loading user function library: [ElmerIceSolvers]...[ComputeNormalSolver]
OptimizeBandwidth: ---------------------------------------------------------
OptimizeBandwidth: Computing matrix structure for: normal vector...done.
OptimizeBandwidth: Half bandwidth without optimization: 10208
OptimizeBandwidth: ---------------------------------------------------------
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver_Init]
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver]
Floating point exception
I have made a lot of changed to the boundary conditions after the suggestions from Thomas and Matina. The domain is a square cartesian grid where the surface and bed topography is loaded in using the Grid2DInterpolator solver. Perhaps I have labed the boundaries in the wrong order? As I understand it normally 1, 2, 3 and 4 are the sides. Then 5 is the bed and 6 is the surface.
I have attached my .SIF file.
Code: Select all
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! HEADER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
$yearinsec = 365.25*24*60*60
$rhoi = 900.0/(1.0e6*yearinsec^2)
$gravity = -9.81*yearinsec^2
$n = 3.0
$eta = (2.0*100.0)^(-1.0/n)
Header
Mesh DB "." "square-part"
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MATC stuff
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
$ function depthwise(D) {\
if (D(0)-D(1) > 5.0) {\
_depthwise = 2.5E-2;\
} else {\
_depthwise = 10.0E00;\
}\
}
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SIMULATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Simulation
Coordinate System = "Cartesian 3D"
Coordinate Mapping(3) = Integer 1 2 3
Simulation Type ="Steady State"
Extruded Mesh Levels = 4
! Iterations between different solvers
! ------------------------------------
Steady State Max Iterations = 1 ! one round isenough, as no coupling between solvers
Steady State Min Iterations = 1
Initialize Dirichlet Conditions = Logical True
! Output files
! ------------
Binary Output = Logical True
Post File = "init.vtk"
Output File = "init.result"
! how verbose the solver should be
! 1 = very low (Finnish style = crucial feedback, only)
! 42 = the answer to all and everything
!-------------------------------------------------------
max output level = 9
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SOLVER
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Flow-depth on the unstructured FEM-mesh
! not really needed, but nice for post-processing
! Just un-comment Exec Solver = Never, in order
! to skip it
!-------------------------------------------------
Solver 1
Exec Solver = Before Simulation
Equation = "Read DEM"
Procedure = "ElmerIceSolvers" "Grid2DInterpolator"
Variable 1 = String "bedrockDEM"
! Variable 1 Position tol = 1e-3
Variable 1 data file = File "bed_data.dat"
Variable 1 x0 = REal 5516.08d0
Variable 1 y0 = REal 12642.84d0
Variable 1 lx = REal 18000.00d0
Variable 1 ly = REal 11334.78d0
Variable 1 Nx = Integer 125
Variable 1 Ny = Integer 79
Variable 1 position tol = REal 1e-1
Variable 2 = String "ZsDEM"
! Variable 2 Position tol = 1e-3
Variable 2 data file = File "surface_data.dat"
Variable 2 x0 = REal 5516.08d0
Variable 2 y0 = REal 12642.84d0
Variable 2 lx = REal 18000.00d0
Variable 2 ly = REal 11334.78d0
Variable 2 Nx = Integer 125
Variable 2 Ny = Integer 79
Variable 2 position tol = REal 1e-1
End
Solver 2
Equation = "MapCoordinate"
Procedure = "StructuredMeshMapper" "StructuredMeshMapper"
Active Coordinate = Integer 3
Mesh Velocity Variable = String "dSdt"
Mesh Update Variable = String "dS"
Mesh Velocity First Zero = Logical True
End
Solver 3
Equation = "Velocity Preconditioning"
Procedure = "DummyRoutine" "DummyRoutine"
Variable = -dofs 3 "V"
Variable Output = False
Exec Solver = "before simulation"
! Exec Solver = "Never"
Element = "p:1 b:4"
Bubbles in Global System = False
Variable Output = False
Bubbles in Global System = False
Linear System Row Equilibration = Logical True
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStabL
Linear System Max Iterations = 250
Linear System Convergence Tolerance = 1.0e-6
Linear System Preconditioning = ILU0
Skip Compute Nonlinear Change = Logical True
Back Rotate N-T Solution = Logical False
End
Solver 4
Equation = "Pressure Preconditioning"
Procedure = "DummyRoutine" "DummyRoutine"
Variable = -dofs 1 "P"
Variable Output = False
! Exec Solver = "Never"
Exec Solver = "before simulation"
Element = "p:1 b:4"
Bubbles in Global System = False
Linear System Symmetric = True
Linear System Scaling = True
Linear System Solver = iterative
Linear System Iterative Method = BiCGStab ! Changed
! Linear System Iterative Method = CG
Linear System Max Iterations = 1000
Linear System Convergence Tolerance = 1.0e-6
Linear System Preconditioning = Diagonal
Linear System Residual Output = 10
Skip Compute Nonlinear Change = Logical True
Back Rotate N-T Solution = Logical False
Linear System Timing = True
End
Solver 5
Equation = "Stokes"
! Exec Solver = "Never"
Procedure = "ParStokes" "StokesSolver"
Element = "p:1 b:4"
Bubbles in Global System = False
Variable = String "Flow Solution"
Variable Dofs = 4
Convective = Logical False
Exported Variable 1 = -dofs 1 "dSdt"
Exported Variable 2 = -dofs 1 "dS"
Exported Variable 3 = -dofs 1 "ZsDEM"
Exported Variable 4 = -dofs 1 "bedrockDEM"
Block Diagonal A = Logical True
Use Velocity Laplacian = Logical True
!-----------------------------------------------
! Keywords related to the block preconditioning
!-----------------------------------------------
Block Preconditioning = Logical True
Linear System Scaling = Logical True
Linear System Row Equilibration = Logical True
Linear System Solver = "Iterative"
Linear System GCR Restart = Integer 100
Linear System Max Iterations = 200
! Linear System Adaptive Tolerance = Logical True
! Linear System Relative Tolerance = Real 1e-2
! Linear System Base Tolerance = Real 1.0e-3
Linear System Convergence Tolerance = 1.0e-6
Nonlinear System Max Iterations = 50
Nonlinear System Convergence Tolerance = 1.0e-5
Nonlinear System Newton After Tolerance = 1.0e-3
End
Solver 6
Equation = "Strain Rate"
Exec Solver = "Always"
Procedure = "ElmerIceSolvers" "ComputeStrainRate"
! this is just a dummy, hence no output is needed
!-----------------------------------------------------------------------
Variable = -nooutput "Eij"
Variable DOFs = 1
Exported Variable 1 = "StrainRate"
Exported Variable 1 DOFs = 7
! the name of the variable containing the flow solution (U,V,W,Pressure)
!-----------------------------------------------------------------------
Flow Solver Name = String "Flow Solution"
! the name of the strain-rate solution (default is 'StrainRate')
StrainRate Variable Name = String "StrainRate"
! Linear System Solver = Direct
! Linear System Direct Method = umfpack
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0E-09
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 1
End
Solver 7
Equation = String "StressSolver"
Exec Solver = Always
Procedure = File "ElmerIceSolvers" "ComputeDevStress"
! this is just a dummy, hence no output is needed
!-----------------------------------------------------------------------
Variable = -nooutput "Sij"
Variable DOFs = 1
! the name of the variable containing the flow solution (U,V,W,Pressure)
!------------------------------------------------------------------------
Flow Solver Name = String "Flow Solution"
!-----------------------------------------------------------------------
Exported Variable 1 = "Stress" ! [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
! no default value anymore for "Stress Variable Name"
Stress Variable Name = String "Stress"
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0E-05
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 1
End
Solver 8
Exec Solver = "Before Simulation"
Exec Solver = Never
Equation = "Normal vector"
Variable = "Normal Vector"
! in 3dimensional simulations we have 3 entries
Variable DOFs = 3
!NB: does not need to actually solve a matrix
! hence no BW optimization needed
Optimize Bandwidth = Logical False
Procedure = "ElmerIceSolvers" "ComputeNormalSolver"
! if set to True, all boundary normals would be computed by default
ComputeAll = Logical False
Exported Variable 1 = -dofs 3 "Normal Vector"
End
Solver 9
Exec Solver = Always
Exec Interval = 1
Equation = "result output"
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Exported Variable 1 = -dofs 1 "dSdt"
Exported Variable 2 = -dofs 1 "dS"
Exported Variable 3 = -dofs 1 "ZsDEM"
Exported Variable 4 = -dofs 1 "bedrockDEM"
Output File Name = "Smokey.vtu"
Output Format = String vtu
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(9) = 1 2 3 4 5 6 7 8 9
Convection = Computed ! we have a computed velocity field
Flow Solution Name = String "Flow Solution" ! and that is its name
! NS Convect = False ! and the acceleration terms are skipped (Stokes)
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! INITIAL CONDITIONS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Initial Condition 1
Pressure = Real 0.0
Velocity 1 = Real 0.0
Velocity 2 = Real 0.0
Velocity 3 = Real 0.0
Flow Solution 1 = Real 0.0
Flow Solution 2 = Real 0.0
Flow Solution 3 = Real 0.0
Flow Solution 4 = Real 0.0
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BODY FORCE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Body Force 1
Flow BodyForce 1 = 0.0
Flow BodyForce 2 = 0.0
Flow BodyForce 3 = $gravity
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MATERIAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Material 1
Cauchy = Logical True
Viscosity Model = String "power law"
Viscosity Exponent = Real $1.0/3.0
Critical Shear Rate = Real 1.0e-10
Viscosity = REAL $eta
Density = REAL $rhoi
End
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BOUNDARY CONDITIONS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! the artificial side wall:
Boundary Condition 1
Target Boundaries(4) = Integer 1 2 3 4
Name = "sides"
! no slip
! ------------------------
Flow Solution 1 = real 0.0
Flow Solution 2 = real 0.0
Flow Solution 3 = real 0.0
Flow Solution 4 = Variable ZsDEM, bedrockDEM
Real MATC " rhoi * gravity * (tx(0)-tx(1))"
V 1 = Real 0.0
V 2 = Real 0.0
V 3 = Real 0.0
End
! surface:
Boundary Condition 2
Target Boundaries = Integer 6
Name = "surf"
Top Surface = Variable ZsDEM, bedrockDEM
Real MATC "if (tx(0)>tx(1)+1.0) {tx(0)} else {tx(1)+1.0}"
End
!! bedrock:
Boundary Condition 3
Target Boundaries = Integer 5
Name = "bedrock"
Bottom Surface = Equals bedrockDEM
Normal-Tangential V = Logical True
Normal-Tangential Flow Solution = Logical True
Flow Force BC = Logical True
V 1 = Real 0.0
Flow Solution 1 = Real 0.0
! Slip Coefficient 2 = Variable Coordinate 1
! Real Procedure "ElmerIceUSF" "Friction_Coulomb"
! Slip Coefficient 3 = Variable Coordinate 1
! Real Procedure "ElmerIceUSF" "Friction_Coulomb"
!! Parameters needed for the Coulomb Friction Law
!! Friction Law Sliding Coefficient = Real 4.1613e5
! Friction Law Post-Peak Exponent = Real 1.0 !(q=1)
!! Friction Law Maximum Value = Real 1.0 !(C=1)
! Friction Law Maximum Value = Variable ZsDEM, bedrockDEM
! Real MATC "depthwise(tx)"
! Friction Law PowerLaw Exponent = Real 3.0 !(m = n = 3 Glen's law)
! Friction Law Linear Velocity = Real 0.01
Slip Coefficient 2 = Variable Coordinate 1
Real Procedure "ElmerIceUSF" "Sliding_Weertman"
Slip Coefficient 3 = Variable Coordinate 1
Real Procedure "ElmerIceUSF" "Sliding_Weertman"
Weertman Friction Coefficient = Real 10.0
! Weertman Friction Coefficient = Real 2E-2
! Weertman Friction Coefficient = Variable ZsDEM, bedrockDEM
! Real MATC "depthwise(tx)"
Weertman Exponent = Real $1.0/3.0
Weertman Linear Velocity = Real 0.00001
End
Any clue where the problem could be? I have run out of new things to try.
Thanks,
Christian