trouble going from diagnostic to prognostic

Extension of Elmer in computational glaciology
Post Reply
marsice
Posts: 4
Joined: 08 Mar 2015, 06:25
Antispam: Yes

trouble going from diagnostic to prognostic

Post by marsice »

Hello,

I am trying to do a prognostic case of a glacier viscously flowing over time with no accumulation. My diagnostic case works perfectly. I have tried to follow the wiki in changing it a prognostic case but it doesn't seem to work - the solver runs but the free surface doesn't seem to actually be free. I suspect I am missing something obvious but can't figure it out. The problem may be in my boundary condition for the free surface:

Boundary Condition 3
Name = "surface"
Target Boundaries = 2
Body ID = 2
Top Surface = Variable Coordinate 1
Real cubic
include "icescarp70degground.dat"
End
Depth = Real 0.0
End

Where "icescarp70degground.dat" is the file containing the geometry of the initial glacier that is to flow. Should this instead be in the initial condition somewhere? I have attached the prognostic case sif (which works, thanks for the helpful documentation and notes!) and the diagnostic case sif.
Attachments
icescarp70degground.dat
geometry of "glacier"
(322.27 KiB) Downloaded 1945 times
Mars_scarp_time.sif
prognostic case, free surface doesn't work
(7.59 KiB) Downloaded 1990 times
Mars_scarp.sif
diagnostic case, which works
(5.62 KiB) Downloaded 1990 times
marsice
Posts: 4
Joined: 08 Mar 2015, 06:25
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by marsice »

Also attached are the msh files and the bedrock dat file... thanks in advance for any help
Attachments
meshfile.zip
msh file
(215.97 KiB) Downloaded 2127 times
groundwide.dat
bedrock dat file
(322.27 KiB) Downloaded 1963 times
meshstuff.zip
mesh files (elements boundary header nodes)
(225.24 KiB) Downloaded 2019 times
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by tzwinger »

The main issue is conected to the fact that your "Top Surface" keyword points to the initial data all the time, but rather should point to the result of your free-surface evolution equation. So, the good strategy would be to let the initial shape define the initial condition for the free-surface variable and link, as said before, the latter with "Top Surface"

Regards,

Thomas
marsice
Posts: 4
Joined: 08 Mar 2015, 06:25
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by marsice »

OK, thank you, that makes sense. I have changed my boundary condition for the top surface to just:

Boundary Condition 3
Name = "surface"
Target Boundaries = 2
Body ID = 2
Top Surface = Equals "Zs"
Depth = Real 0.0
End

To point to the result of the free-surface evolution equation. I think my problem is just that I am not putting in the initial condition correctly. If I do something like:

Initial Condition 2
Zs = Equals Variable Coordinate 2
Real cubic
include "icescarp70degground.dat"
End
RefZs = Equals Coordinate 2
End

I get an error for "Unknown Specifier" for property name "Real cubic." The same happens if I instead try and write a MATC expression in the initial condition. I think I must be missing something very easy - is this the incorrect syntax to input the initial shape? Or should it go in a different section, other than initial condition?
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by tzwinger »

Code: Select all

Initial Condition 2
Zs = Variable Coordinate 2
Real cubic
include "icescarp70degground.dat"
End
RefZs = Equals Coordinate 2
End
should do it. So no "Equals" for the Zs-initialization.

And then, of course, you have to take care that the solver that updates the mesh (StructuredMeshMapper) is executed before any other solver in the time-loop.

Regards,
Thomas
marsice
Posts: 4
Joined: 08 Mar 2015, 06:25
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by marsice »

OK, thank you! I am getting much closer. Elmer does run with that change, but the output does not have the glacier geometry, only the bedrock. I can get results with the correct glacier geometry using the following sif:

Code: Select all

Header
  Mesh DB "." "testglacier2"
  Include Path ""
  Results Directory ""
End


Simulation
  Max Output Level = 4
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = "Transient"
  Steady State Max Iterations = 1
  Timestepping Method = "BDF"
  BDF Order = 1
  Timestep Sizes = 1
  Timestep Intervals = 3
  Output Intervals = 1
  Output File = "scarp_thirdtest.result"
  Post File = "scarp_thirdtest.vtu" ! use .ep suffix for leagcy format
  Initialize Dirichlet Conditions = Logical False
End

Constants
  Gravity(4) = 0 -1 0 3.711
  Stefan Boltzmann = 5.67e-08
End

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

Body 2
  Name = "free surface"
  Equation = 2
  Material = 2
  Body Force = 2
  Initial Condition = 2
End

Equation 1
  Name = "Equation1"
  Convection = "computed"
  Flow Solution Name = String "Flow Solution"
  Active Solvers(5) = 1 2 3 4 6
End

Equation 2
   Active Solvers(1) = 5
   Flow Solution Name = String "Flow Solution"
   Convection = Computed
End

Initial Condition 1
  Velocity 1 = 0.0
  Velocity 2 = 0.0
  Pressure = 0.0
  Depth = Real 0.0
End

Initial Condition 2
   Zs = Variable Coordinate 2
   Real cubic
       include "icescarp70degground.dat"
   End
   ReferenceZs = Equals Coordinate 2
   Mesh Update 2 = Real 0.0
   Mesh Update 1 = Real 0.0
End


Solver 1
  Exec Solver = "before simulation"
  Equation = "MapCoordinate"
  Procedure = "StructuredMeshMapper" "StructuredMeshMapper"
  Active Coordinate = Integer 2 ! the mesh-update is y-direction
  Mesh Velocity Variable = String "Mesh Velocity 2"
  Mesh Velocity First Zero = Logical True
  Dot Product Tolerance = Real 0.01
End

! FlowDepth on generally unstructured mesh
!  (could be replaced by strucuterd version)
Solver 2
   Equation = "Flowdepth"
   Exec Solver = "Before Simulation"
   Procedure = File "ElmerIceSolvers" "FlowDepthSolver"
   Variable = String "Depth"
   Variable DOFs = 1
   Linear System Solver = "Direct"
   Linear System Direct Method = "UMFPACK"
   Linear System Max Iterations = 200
   Linear System Convergence Tolerance = 1.0E-09
   Linear System Abort Not Converged = False
   Linear System Preconditioning = "ILU0"
   Linear System Residual Output = 1
   Gradient = Real -1.0E00
   Exported Variable 1 = -dofs 3 "Mesh Velocity" 
   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 3
  Equation = "Navier-Stokes"
  Optimize Bandwidth = Logical True
  Linear System Solver = Direct
  Linear System Direct Method = "UMFPACK"   
  Linear System Max Iterations = 5000
  Linear System Convergence Tolerance = 1.0E-06
  Linear System Abort Not Converged = False
  Linear System Preconditioning = "ILU1"
  Linear System Residual Output = 1

  Steady State Convergence Tolerance = 1.0E-05
  Stabilization Method = Stabilized

  Nonlinear System Convergence Tolerance = 1.0E-04
  Nonlinear System Convergence Measure = Solution
  Nonlinear System Max Iterations = 50
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance =  1.0E-01
End

Solver 4
  Equation = Sij
  Procedure = "ElmerIceSolvers" "ComputeDevStress"
  Variable = -nooutput "Sij"
  Variable DOFs = 1
  Exported Variable 1 = -dofs 4 Stress
  Stress Variable Name = String "Stress"
  Flow Solver Name = String "Flow Solution"
  Linear System Solver = Direct
  Linear System Direct Method = umfpack
End

Solver 5
   Exec Solver = "After TimeStep"
   Equation = String "Free Surface Evolution"
   Variable = "Zs"
   Variable DOFs = 1
   Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"
   Apply Dirichlet = Logical True
   Linear System Solver = Iterative
   Linear System Iterative Method = BiCGStab
   Linear System Max Iterations = 100
   Linear System Preconditioning = ILU1
   Linear System Convergence Tolerance = 1.0e-8

   Nonlinear System Max Iterations = 100
   Nonlinear System Min Iterations = 2
   Nonlinear System Convergence Tolerance = 1.0e-6
  
   Steady State Convergence Tolerance = 1.0e-4

   Stabilization Method = Bubbles
   Flow Solution Name = String "Flow Solution"

   Exported Variable 1 = Zs Residual
   Exported Variable 1 DOFS = 1

   Exported Variable 2 = ReferenceZs
   Exported Variable 2 DOFS = 1
End

Solver 6
   Equation = "Mesh Update"
   Exec Solver = "After TimeStep"
   Linear System Solver = Iterative
   Linear System Iterative Method = BiCGStab
   Linear System Max Iterations = 500
   Linear System Preconditioning = ILU1
   Linear System Convergence Tolerance = 1.0e-6
   Nonlinear System Max Iterations = 1
   Nonlinear System Convergence Tolerance = 1.0e-06
End

! we use m-yr-MPa system 1 yr = 31556926.0 sec
Material 1
  Name = "ice-ice-baby"
  Density = Real $910.0*1.0E-06*(31556926.0)^(-2.0)
  Viscosity Model = String "Glen"
  Viscosity = Real 1.0 
  Glen Exponent = Real 3.0
  Critical Shear Rate = Real 1.0e-10
  Rate Factor 1 = Real 1.258e13  
  Rate Factor 2 = Real 6.046e28
  Activation Energy 1 = Real 60e3
  Activation Energy 2 = Real 139e3  
  Glen Enhancement Factor = Real 1.0
  Mesh Elastic Modulus = 1.0
  Mesh Poisson Ratio = 0.3
  ! 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 "Temp Homologous"
  ! the temperature to switch between the 
  ! two regimes in the flow law
  Limit Temperature = Real 0.0
  Constant Temperature = Real MATC "220-273.15"
  Cauchy = Logical False
End

Material 2
  Density = Real $910.0*1.0E-06*(31556926.0)^(-2.0)
!  Min Zs = Real 0.0
!  Max Zs = Real 900.0
End

Body Force 1
  Name = "BodyForce1"
  Heat Source = 1
  Flow BodyForce 1 = Real 0.0                          
  Flow BodyForce 2 = Real -3.711e15  !MPa - a - m
End
 
Body Force 2
  Zs Accumulation Flux 1 = Real 0.0
  Zs Accumulation Flux 2 = Real 0.0
End


Boundary Condition 1
  Name = "bedrock"
  Target Boundaries = 1
  Compute Normals = Logical True
  Bottom Surface = Variable Coordinate 1
  Real cubic
     include  "groundwide.dat" 
  End
  Velocity 1 = Real 0.0e0
  Velocity 2 = Real 0.0e0
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
End

Boundary Condition 2
  Name = "sides"
  Target Boundaries(2) = 3  4 ! combine left and right boundary
  Velocity 1 = 0.0
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
End

Boundary Condition 3
  Name = "surface"
  Target Boundaries = 2
  Body ID = 2
  Top Surface = Variable Coordinate 1  
  Real cubic
     include "icescarp70degground.dat" 
  End
! Top Surface = Equals Zs
  Depth = Real 0.0
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Variable Zs, ReferenceZs
      Real MATC "tx(0) - tx(1)"
End
However, I don't think this is quite correct, because as you pointed out before the surface Boundary Condition is pointing toward what I want to be the initial state. I have played around a lot with when the StructuredMeshMapper is executed like you suggested, and with a Mesh Update solver, but think I am probably using them incorrectly - is it necessary to have Mesh Update terms in the BC, and would it make sense to have the top surface term inside one of the solvers instead of the BC?
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: trouble going from diagnostic to prognostic

Post by tzwinger »

Hello,
I didn't realize your MeshUpdate in the first instance - sorry. Then you can do it such that you just initialize, imposing the data directly to Top Surface and then couple the FS to the MeshUpdate variable, sacrificing the advantage of having a vertically structured mesh. Or, you do as I suggested and link the Top Surface to FS - your choice. In any way, you should initialize the ReferenceFS (now set to Y-coords in the unaltered mesh) with the same values as the initial FS.

Best wishes,

Thomas
Post Reply