Mask to restrict sliding in ice free areas

Extension of Elmer in computational glaciology
ChristianB
Posts: 10
Joined: 28 Apr 2014, 22:02
Antispam: Yes

Mask to restrict sliding in ice free areas

Post by ChristianB »

Hi Everyone,

I am trying to get columb-friction sliding running on a model of a small'ish glacier I am currently working on.
I am using the ParStokes solver but it keeps diverging when I activate the CF law. From reading on the forum I have been trying to write a mask to restrict sliding in areas with less then 10 meters of ice. The model contains some areas that are ice free where I set the ice to 1 m,

Code: Select all

  Top Surface = Variable ZsDEM, bedrockDEM                       
    Real MATC "if (tx(0)>tx(1)+1.0) {tx(0)} else {tx(1)+1.0}"  
If I understand this post correct: viewtopic.php?f=3&t=281#p908
This syntax should set Flow Solution 1 equal zero if the ice is less then 10 meters.

Code: Select all

  
  Flow Solution 1 = Real 0.0 
  Flow Solution 1 Condition = Variable ZsDEM, bedrockDEM 
    Real MATC "-1*(tx(0)-tx(1)-10.0)" 
Dirichlet is active as default. But the solver is still not converging. Anyone know the best way to do CF sliding with ice-free areas? I am still new to Elmer/ICE so it may just be a beginner problem.

My full bedrock BC is,

Code: Select all

!! bedrock:                                                           
Boundary Condition 3                                                  
  Name = "bedrock" 
  Height = Real 0.0                                                  
  Bottom Surface = Equals bedrockDEM                               

  Normal-Tangential Velocity = Logical True                           
  Flow Force BC = Logical True                                        
                                                                      
  Flow Solution 1 = Real 0.0                                          
  Flow Solution 1 Condition = Variable ZsDEM, bedrockDEM              
    Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                      
  Velocity 1 = REal 0.0                                               
  Velocity 1 Condition = Variable ZsDEM, bedrockDEM                   
    Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                      
  Velocity 2 = REal 0.0                                               
  Velocity 2 Condition = Variable ZsDEM, bedrockDEM                   
    Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                     
  Flow Solution 2 = REal 0.0                                          
  Flow Solution 2 Condition = Variable ZsDEM, bedrockDEM             
    Real MATC "-1*(tx(0)-tx(1)-10.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 
  Friction Law Maximum Value = Real 1.0  
  Friction Law PowerLaw Exponent = Real 3.0
  Friction Law Linear Velocity = Real 0.01  
End 
Thanks,
Christian
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by tzwinger »

Hi Christian,
I made an approach with a mountain glacier using Weertman law. I didn't switch between a Dirichlet condition and the sliding law, but simply, depending on the depth, switched from one to another sliding coefficient. Something like that:

Code: Select all

$ function depthwise(D) {\
   if (D > 10.0) {\
               _depthwise = 2.5E-02;\
   } else {\
         _depthwise = 10.0E00;\
  }\
}
With the call in the boundary condition:

Code: Select all

  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 = Variable Depth
    Real MATC "depthwise(tx)" 
So a high enough coefficient basically keeps the ice from sliding. For me there were convergence issues with this (but I did not try the Dirichlet-Robin switching approach)

Cheers,

Thomas
ChristianB
Posts: 10
Joined: 28 Apr 2014, 22:02
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by ChristianB »

Thanks for the suggestion thomas. Unfortunately it did not help.

To reduce the complexity of my model I changed to a Weertman sliding law and forced a constant sliding coefficient of 10.0E00. That should in practis result in a no-sliding condition everywhere. But also in this "simplere" case the velocity does not converge.
I should note that the residual in the new model is smaller (~1 - 100) whereas it before was much higher (~1E20 - 1E40). But no matter what I try it does not come close to converge (1E-6).
Without the sliding boundary condition the model converges with no problems.

Perhaps the sliding boundary is causing problems with my solver setup? Are there any recommended settings for using the ParStokes with a sliding law? I tried switching from a CG algorithm to a BiCGStab which helped a bit.
My current solver for the velocity system is a follows,

Code: Select all

Solver 3                                                              
  Equation = "Velocity Preconditioning"                               
  Procedure = "DummyRoutine" "DummyRoutine"                           
  Variable = -dofs 3 "V"                                              
  Variable Output = False                                             
  Exec Solver = "before simulation"                                   
  Element = "p:1 b:4"                                                 
  Bubbles in Global System = False                                    
                                                                      
  Linear System Symmetric = True                                      
  Linear System Scaling = True                                        
  Linear System Row Equilibration = Logical True                     
  Linear System Solver = Iterative                                    
  Linear System Iterative Method = BiCGStab                           
  Linear System Max Iterations = 250                                  
  Linear System Preconditioning = ILU0                                                                           
  Linear System Convergence Tolerance = 1.0e-6                        
  Linear System Abort Not Converged = False                           
  Skip Compute Nonlinear Change = Logical False                       
  Back Rotate N-T Solution = Logical False                            
  Linear System Timing = True                                         
End 
Cheers,
Christian
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by tzwinger »

Christian,
I didn't realise that you use Block Pre-Conditioner (BPC). As you apparently were able to run with no-slip condition (Dirichlet) the issue is not directly connected to your Mesh or Solver rather than the boundary condition. Solver settings for ParStokes are well documented here: http://elmerice.elmerfem.org/wiki/lib/e ... stokes.pdf. Unfortunately, there is only Dirichlet at bedrock in this document. Check also, whether your bubble-degrees are consistent with your element type(s) (but I guess that should not be the issue).
If you haven't already done so, could you simply switch to linear sliding law (so insert a number or fixed distribution for the sliding coefficient) and see what you get?
I will try to dig out a case where sliding is used in combination with BPC.

Best,

Thomas
tzwinger
Site Admin
Posts: 99
Joined: 24 Aug 2009, 12:20
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by tzwinger »

Hi Christian,
I apparently missed a few potential issues in your setup:

Code: Select all

!! bedrock:                                                           
    Boundary Condition 3                                                 
      Name = "bedrock"
      Height = Real 0.0                                                 
      Bottom Surface = Equals bedrockDEM                               

      Normal-Tangential Velocity = Logical True                           
      Flow Force BC = Logical True                                       
                                                                         
      Flow Solution 1 = Real 0.0                                         
      Flow Solution 1 Condition = Variable ZsDEM, bedrockDEM             
        Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                         
      Velocity 1 = REal 0.0                                               
      Velocity 1 Condition = Variable ZsDEM, bedrockDEM                   
        Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                         
      Velocity 2 = REal 0.0                                               
      Velocity 2 Condition = Variable ZsDEM, bedrockDEM                   
        Real MATC "-1*(tx(0)-tx(1)-10.0)"                                 
                                                                         
      Flow Solution 2 = REal 0.0                                         
      Flow Solution 2 Condition = Variable ZsDEM, bedrockDEM             
        Real MATC "-1*(tx(0)-tx(1)-10.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
      Friction Law Maximum Value = Real 1.0 
      Friction Law PowerLaw Exponent = Real 3.0
      Friction Law Linear Velocity = Real 0.01 
    End
First of all, you should be sure that the Variable name in your ParStokes solver is "Flow Solution" and has Variable DOFs = 4
Then you should change your "Normal-Tangential Velocity = Logical True" to a "Normal-Tangential Flow Solution = Logical True"
Then, if the switch above is set, your normal velocity (in direction of the normal, hence pointing into the bedrock) is "Flow Solution 1" - that one you want to (apart from prescribing melting/refreezing) have set to zero, i.e.,

Code: Select all

Flow Solution 1 = Real 0.0
without any condition to it (else you let it free falling out of your domain for thick enough glacier). You want to set the same condition (Normal-Tangential and 1st component zero also for the variable in your pre-conditioning routine (I guess it was V).

Components 2 and 3 then can be set as they are - at least in the example I found, where I use linear sliding coefficients, that worked. As mentioned before, start with simple sliding law and - once that works - use the fancy stuff.

Hope that helps,

Thomas
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by Martina »

Hi Christian,

I am using ParStokes together with the normal siding law and it works. I attach you below some lines of my sif-file.
However, in that case I didn't have convergence problems so I did not need to use a mask.
I had similar problems in another case (but without using ParStokes), so I was wondering if you could try with the normal Navier-Stokes Solver to check if your problem is really linked to ParStokes or to your setup as a such? I had a case where it is not possible to use a sliding law with a the same constant parameter everywhere.

Martina

P.S. Just saw that Thomas was posting while I was writing, so sorry for partly double posting....

Code: Select all

Solver 3
  Equation = "Stokes equations"
  Procedure = "./ParStokes" "StokesSolver"
  Variable = FlowVar
  Variable Dofs = 4
  Element = "p:1 b:4"  
  Bubbles in Global System = False
  Nonlinear System Convergence Tolerance = 1.0e-5
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Tolerance = 1.0e-2
  Linear System Row Equilibration = Logical True
  Linear System Solver = "Iterative"
  Linear System Iterative Method = "GCR"
  Linear System Max Iterations = 200
  Linear System GCR Restart = Integer 50
  Linear System Convergence Tolerance = 1.0e-6
  Block Preconditioning = Logical True
  Linear System Adaptive Tolerance = Logical True
  Linear System Base Tolerance = Real 1.0e-3
  Linear System Relative Tolerance = Real 1.0e-2
 End


Solver 8
  Equation = "Velocity Preconditioning"
  Procedure = "VelocityPrecond" "VelocityPrecond"
  Exec Solver = "before simulation"
  Variable = -dofs 3 "V"
  Variable Output = False
  Element = "p:1  b:4"
  Bubbles in Global System = False
  Linear System Row Equilibration = Logical True
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStabL
  Linear System Max Iterations = 1000
  Linear System Convergence Tolerance = 1.0e-3
  Linear System Preconditioning = ILU0
  Skip Compute Nonlinear Change = Logical True
  Back Rotate N-T Solution = Logical False
End


Boundary Condition 3
  Name = "Bottom"
  ....
  Normal-Tangential FlowVar = Logical True
  Normal-Tangential V = Logical True
  V 1 = Real 0.0e0
  FlowVar 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
ChristianB
Posts: 10
Joined: 28 Apr 2014, 22:02
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by ChristianB »

Sorry for the late reply. Have been without internet for a couple of weeks.

Thanks for all the suggestions. I have all week to play with Elmer so I will try the suggestions and write back about the outcome.

Martina: I will try the normal Navier-Stokes solver at some point and see what happens.
ChristianB
Posts: 10
Joined: 28 Apr 2014, 22:02
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by ChristianB »

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.

Image

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
Martina
Posts: 39
Joined: 17 Apr 2013, 14:06
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by Martina »

Hi Christian,

I have the feeling that something is wrong with your mesh extrusion. Did you try to run some very simple solver with an obvious result (for example FlowDepth) to check?
On the wiki it says "Note: “Target Boundaries” must not be given in the bottom and top boundary sections in the sif-file, however they need to be Boundary Condition number i+1 and i+2 (i being the last lateral boundary of the 2D-footprint mesh)." http://elmerice.elmerfem.org/wiki/doku. ... cturedmesh If this comment is still valid, you have to delete these lines in your sif-file.
And finally you should check for example with ElmerPost if your lateral boundaries are really the boundaries 1 2 3 4, there is no general rule which boundary has which number!
Also I'm a bit puzzled that your DOFs of dS, dSdt and so on are 1 in your "exported variables..." lines, aren't you in 3D? And you have these export commands twice (in 2 solvers).

Good luck,
Martina
ChristianB
Posts: 10
Joined: 28 Apr 2014, 22:02
Antispam: Yes

Re: Mask to restrict sliding in ice free areas

Post by ChristianB »

Thanks for the suggestion Martina.

I have tried playing around with the mesh but it is quite difficult debugging mesh problems. The flowdepth solver does not work at the moment (had it working previously) so there must be something with the mesh. I have tried using ElmerPOST but the mesh is quite big and with all the numbers it is completly impossible to move around and see which numbers are where. In my output when running ElmerSolver I get

Code: Select all

WriteVtuFile: Writing body and BC indexes
Anyone know how I can see boundary numbering in Paraview?
Post Reply