Hi,
I've been working with Elmer this year to develop a 2D flowline model of a glacier in Greenland. I've hit a bit of a stumbling block with regards to modelling the additional mass input at the point where a tributary glacier joins the main channel.
Is there a way to artificially create a 'mass source' by adding mass at individual nodes in my 2D mesh to account for the confluence of this tributary?
I thought about simply adding an accumulation flux at the surface BC, but I don't think this would produce realistic results.
Any help is greatly appreciated,
Thanks,
Joe
Modelling ice tributary in 2D flowline case

 Site Admin
 Posts: 3264
 Joined: 22 Aug 2009, 11:57
 Antispam: Yes
 Location: Espoo, Finland
 Contact:
Re: Modelling ice tributary in 2D flowline case
Hi Joe
Well, if you want a nodal load on node then basically you can do this by the following construct in 2D. It sets a nodal load at the node closest to origin for the continuity equation.
Peter
Well, if you want a nodal load on node then basically you can do this by the following construct in 2D. It sets a nodal load at the node closest to origin for the continuity equation.
Code: Select all
Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 0.0 0.0
Flow Solution Loads 3 = Real 1.23
End
Re: Modelling ice tributary in 2D flowline case
Hi Peter,
Thanks for your quick reply. I'm having some trouble implementing your proposed solution, unfortunately. Would you mind having a quick look over my SIF to see if you can spot any obvious mistakes?
I've been trying to pinpoint the problem. I'm able to implement a nodal temperature load using
but, for some reason, it seems to ignore "Flow Solution Loads 3"
Thanks again for your help. Much appreciated!
Joe
Thanks for your quick reply. I'm having some trouble implementing your proposed solution, unfortunately. Would you mind having a quick look over my SIF to see if you can spot any obvious mistakes?
I've been trying to pinpoint the problem. I'm able to implement a nodal temperature load using
Code: Select all
Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 70000.0 0.0
Temp Load = Real 1.23
End
Code: Select all
$Rho = 9.03760022565934E019
!$RhoW = 1.00417780285104E018
$RhoW = 1.03430313693657E018
$g = 9.85098424596868E015
!For Glens
$GlenN = 3.0
!For Width
$N = 3.0
$Bforce = 6.0E9
!For basal slip
$As = 4.0E000
$q = 2.0
$C = 0.42
$ut0 = 50.0
$m = 3.0
!Piezometric stuff
$L = 109550.0
$Piez = 0.015
Header
Mesh DB "." "mesh"
Include Path ""
Results Directory "." "Results"
End
Constants
End
Simulation
Coordinate System = "Cartesian 2D"
Simulation Type = SteadyState
Coordinate Mapping(3) = 1 2 3
Steady State Min Iterations = 1
Steady State Max Iterations = 5
Set Dirichlet BCs By BC Numbering = Logical True
max output level = 3
End
Body 1
name = All
Equation = 1
Material = 1
Body force = 1
Initial Condition = 1
End
Material 1
Density = Real $ Rho
Viscosity Model = Power Law
Viscosity Exponent = $ 1.0/3.0
Critical Shear Rate = $ 3.0e09*31556926.0
Viscosity = Variable Temp
Real MATC "glen(tx) * 31556926.0^(1.0/3.0)*1.0E06"
! the heat capacity as a MATC function of temperature itself
!
Temp Heat Capacity = Variable Temp
Real MATC "capacity(tx)*(31556926.0)^(2.0)"
! the heat conductivity as a MATC function of temperature itself
!
Temp Heat Conductivity = Variable Temp
Real MATC "conductivity(tx) * (31556926) * 1.0E06"
! Upper limit  pressure melting point
! as a MATC function of the pressure (what else?)
!
Temp Upper Limit = Variable Pressure
Real MATC "pressuremeltingpoint(tx)"
! lower limit (to be save) as 0 K
!
Temp Lower Limit = Real 0.0
Cauchy = Logical False
End
!!!!! SOLVERS
Solver 1
Equation = "Depth"
Exec Solver = "Before TimeStep"
Procedure = File "src/DepthSolver" "DepthSolver"
Variable = String "Depth"
Variable DOFs = 1
Gradient = Real 1.0E00
Linear System Direct Method = "Mumps"
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End
Solver 2
Equation = "Elevation"
Exec Solver = "Before TimeStep"
Procedure = File "src/DepthSolver" "DepthSolver"
Variable = String "Elevation"
Variable DOFs = 1
Gradient = Real 1.0E00
! Linear System Solver = "Direct"
Linear System Direct Method = "Mumps"
! Linear System Direct Method = umfpack
Linear System Solver = Iterative
Linear System Max Iterations = 300
Linear System Iterative Method = "BiCGStab"
Linear System Preconditioning = ILU0
Linear System Convergence Tolerance = Real 1.0E6
Linear System Abort Not Converged = False
Linear System Residual Output = 0
End
Solver 3
Exec Solver = "Before Simulation"
Equation = "Normal vector"
Variable = "Normal Vector"
! in 3dimensional simulations we have 3 entries
Variable DOFs = 2
!NB: does not need to actually solve a matrix
! hence no BW optimization needed
Optimize Bandwidth = Logical False
Procedure = "src/ComputeNormal" "ComputeNormalSolver"
! if set to True, all boundary normals would be computed by default
ComputeAll = Logical True
End
Solver 4
Equation = "NavierStokes"
Procedure = "FlowSolve" "FlowSolver"
Variable = Flow Solution[Velocity:2 Pressure:1]
Exec Solver = Always
Flow model = String "Stokes"
Stabilize = True
! Bubbles = True
Linear System Direct Method = "MUMPS"
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStabl"
BiCGStabl Polynomial Degree = 4
Linear System Max Iterations = 500
Linear System Convergence Tolerance = Real 1.0E6
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU1"
Linear System Residual Output = 0
Nonlinear System Max Iterations = 50
Nonlinear System Convergence Tolerance = 1.0E5
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e4
Nonlinear System Relaxation Factor = 0.999999 !0.999999
Steady State Relaxation Factor = Real 0.9
Steady State Convergence Tolerance = 1.0d3
Exported Variable 1 = Width
Exported Variable 2 = KCoef
Exported Variable 3 = Flow BodyForce 1
Exported Variable 4 = Draggy
Exported Variable 5 = Slip Coefficient 2
Exported Variable 6 = Overburden
Update Exported Variables = Logical True
Nonlinear Update Exported Variables = Logical True
End
Solver 5
Equation = String "StressSolver"
Procedure = File "src/ComputeDevStressNS" "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 = 4 ! 4 in 2D, 6 in 3D
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 300
Linear System Convergence Tolerance = 1.0E09
Linear System Abort Not Converged = True
Linear System Preconditioning = "ILU0"
Linear System Residual Output = 0
Cauchy = Logical False
End
Solver 6
! Exec Solver = "After TimeStep"
Exec Solver = "After Saving"
! Exec Solver = "After All"
Procedure = "SaveData" "SaveLine"
Filename = "surface.dat"
Optimize Node Ordering = Logical False
File Append = Logical False
Variable 1 = Coordinate 1
Variable 2 = Velocity 1
Variable 3 = Slip Coefficient 2
End
Solver 7
Equation = "Result Output"
! Exec Solver = "After Saving"
Exec Solver = "After All"
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = "results"
Output Format = String "vtu"
Vector Field 1 = String "Velocity"
Vector Field 2 = String "Stress"
Scalar Field 1 = String "Depth"
Scalar Field 2 = String "Elevation"
Scalar Field 3 = String "Pressure"
Scalar Field 4 = String "Temp"
Scalar Field 5 = String "Temp Homologous"
Scalar Field 6 = String "Temp Residual"
Scalar Field 7 = String "Velocity 1"
Scalar Field 8 = String "Width"
Scalar Field 9 = String "KCoef"
Scalar Field 10 = String "Flow BodyForce 1"
Scalar Field 11 = String "Draggy"
End
Solver 8
Equation = String "Homologous Temperature Equation"
Procedure = File "src/TemperateIce" "TemperateIceSolver"
! Before Linsolve = "EliminateDirichlet" "EliminateDirichlet"
Variable = String "Temp"
Variable DOFs = 1
Linear System Solver = "Iterative"
Linear System Iterative Method = "BiCGStab"
Linear System Max Iterations = 10000
Linear System Convergence Tolerance = 1.0E06
Linear System Abort Not Converged = False
Linear System Preconditioning = "ILU4"
Linear System Residual Output = 1
Steady State Convergence Tolerance = 1.0E04
Nonlinear System Convergence Tolerance = 1.0E05
Nonlinear System Max Iterations = 50
Nonlinear System Relaxation Factor = Real 0.999 !0.9999
Steady State Relaxation Factor = Real 0.9
Apply Dirichlet = Logical True
Stabilize = True
Exported Variable 1 = String "Temp Homologous"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "Temp Residual"
Exported Variable 2 DOFs = 1
End
!!!!! EQUATION
Equation 1
Active Solvers (8) = 1 2 3 8 4 5 6 7
Flow Solution Name = String "Flow Solution"
Convection = Computed
NS Convect = True
End
!!!!! BOUNDARY CONDITIONS
Boundary Condition 1
name = Base
Target Boundaries = 2
Flow Force BC = Logical True
NormalTangential Velocity = Logical True
Velocity 1 = 0.0
External Pressure = Variable Coordinate 1, Coordinate 2
Real MATC "getPw(tx(0),tx(1))"
Slip Coefficient 2 = Variable Coordinate 1
Real Procedure "./src/USF_Sliding" "Friction_jgr2007"
Overburden = Variable Depth
Real MATC "tx*Rho*g"
Friction Law Sliding Coefficient = Real $As
Friction Law PostPeak Exponent = Real $q
Friction Law Maximum Value = Real $C
Friction Law Linear Velocity = Real $ut0
Friction Law PowerLaw Exponent = Real $m
Temp Flux BC = Logical True
Temp Heat Flux = Real 55.0E03
Elevation = Real 0.0
Save Line = Logical True
End
Boundary Condition 2
name = Surface
Target Boundaries = 3
Temp = Real 253.15
Depth = Real 0.0
Save Line = Logical True
End
Boundary Condition 3
name = Right
Target Boundaries = 1
Flow Force BC = Logical True
External Pressure = Variable Coordinate 2
Real MATC "if(tx>10.0){0.0}else{if(tx<0.0){RhoW*g*tx}else{0.0}}"
! Real MATC "if(tx>0.0){0.0}else{RhoW*g*tx}"
! External Pressure = Variable Coordinate 2, Depth, Elevation
! Real MATC "Rho*g*(tx(1))+0.5*Rho*g*(1.0Rho/RhoW)*(tx(1)+tx(2))"
Save Line = Logical True
End
Boundary Condition 4
name = Left
Target Boundaries = 4
Flow Force BC = Logical True
External Pressure = Variable Depth, Elevation
Real MATC "0.9635 * Rho*g*tx(0)"
!Real MATC "1.0E5 * tx(1)/(tx(0)+tx(1))"
Save Line = Logical True
End
Boundary Condition 5
Name = "Nodal Mass Source"
Target Coordinates(1,2) = 70000.0 0.0
Flow Solution Loads 3 = Real 5.0E10
End
Body Force 1
Width = Variable Coordinate 1
Real
include Halfwidth.dat
End
KCoef = Variable Width, Temp
Real MATC "((N+1)^(1/N))/((tx(0)^((N+1)/N))*(2*getArrheniusFactor(tx(1))^(1/N)))"
Draggy = Variable Velocity 1
Real MATC "if(tx==0){0.0}else{tx * abs(tx)^((1.0/N)1.0)}"
Flow BodyForce 1 = Variable KCoef, Draggy
Real MATC "Bforce*tx(0)*tx(1)"
!!2.8E9 works.
Flow BodyForce 2 = Real $ g
Temp Volume Source = Real 1.0E00
End
Initial Condition 1
Temp = Real 263.15
End
!Functions
!! pressuremeltingpoint (in SI units)
$ function pressuremeltingpoint(PIN) {\
P = PIN;\
if (P<0.0) P=0.0;\
beta=9.8E08*1.0E06;\
_pressuremeltingpoint=273.15(beta*P);\
}
!! in SI units, input in Kelvin
$ function conductivity(T) { _conductivity=9.828*exp(5.7E03*T)}
!! in SI units, input in Kelvin
$ function capacity(T) { _capacity=146.3+(7.253*T)}
!! Glen's flow law
!
$ function glen(Th) \
import GlenN {\
EF = 1.0;\
AF = getArrheniusFactor(Th);\
_glen = (2.0*EF*AF)^(1.0/GlenN);\
}
!! Arrhenius factor needed by glen
!! (in SI units)
!
$ function getArrheniusFactor(Th){ \
Th1 = Th  273.15;\
if (Th1<10) {_getArrheniusFactor=3.985E13 * exp( 60.0E03/(8.314 * (273.15 + Th1)));}\
else {\
if (Th1>0) _getArrheniusFactor=1.916E03 * exp( 139.0E03/(8.314 * (273.15)));\
else _getArrheniusFactor=1.916E03 * exp( 139.0E03/(8.314 * (273.15 + Th1)));}\
}
!! Basal Water Pressure  calculates based on a given piezometric gradient.
!
!$ function getPw(C1,Depth) \
! import L,Base,RhoW,g,Divide { \
! x = (L  C1)/L;\
!y = Depth * (Base  ((Base  Divide)*x));\
!if (y>0) {_getPw = RhoW*g*y;}\
! else {_getPw = 0.0;}\
!}
$ function getPw(C1,C2) \
import L,Piez,RhoW,g { \
x = L  C1;\
y = C2 + (x*Piez);\
if (y>0) {_getPw = RhoW*g*y;}\
else {_getPw = 0.0;}\
}
Joe

 Site Admin
 Posts: 3264
 Joined: 22 Aug 2009, 11:57
 Antispam: Yes
 Location: Espoo, Finland
 Contact:
Re: Modelling ice tributary in 2D flowline case
Hi
Ok, my fault. I mixed two different conventions: the name of computed nodal loads and the name of given nodal loads. These should rather be same, but unfortunately they are not Both features are rather rarely used so it is not a major problem but still a slightly disturbing. So try
Peter
Ok, my fault. I mixed two different conventions: the name of computed nodal loads and the name of given nodal loads. These should rather be same, but unfortunately they are not Both features are rather rarely used so it is not a major problem but still a slightly disturbing. So try
Code: Select all
Flow Solution 3 Load = Real 5.0e10
Re: Modelling ice tributary in 2D flowline case
Hi Peter,
Still no luck I'm afraid; am I maybe missing a keyword in the NS Solver section or something like that?
Cheers,
Joe
Still no luck I'm afraid; am I maybe missing a keyword in the NS Solver section or something like that?
Cheers,
Joe
Re: Modelling ice tributary in 2D flowline case
Interestingly, Displacement 1 Load doesn't seem to do anything either. Is there something wrong with my NS Equation settings?

 Site Admin
 Posts: 3264
 Joined: 22 Aug 2009, 11:57
 Antispam: Yes
 Location: Espoo, Finland
 Contact:
Re: Modelling ice tributary in 2D flowline case
Hi
Ok, if the vector variable has components that have been named differently, as is the case with FlowSolve, then the name of the Load is derived from that name. I.e. use
There is a new consistency test for this as "pointload2". Sorry for the confusing instructions.
Peter
Ok, if the vector variable has components that have been named differently, as is the case with FlowSolve, then the name of the Load is derived from that name. I.e. use
Code: Select all
Pressure Load = Real 5.0e10
Peter
Re: Modelling ice tributary in 2D flowline case
Yep that did it! Thanks very much for your help, Peter.
Joe
Joe