# Geometry related issues

Jump to: navigation, search

On this page one can find SIF inputs for geometry manipulations. The following properties are given below:

• Long/Lat as a function of projected Cartesian coordinates
• Implementation of a prognostic free surface problem

Long/Lat as a function of stereographic projected Cartesian coordinates

this is an example that is being used for Dome F region, EAIS, with input coordinated from the projection on -71 degrees latitude (with R=6371 km)

```\$ function longitude(X)  { _longitude = atan(X(0)/X(1))*360/(2*pi) }
\$ function latitude(X)  { _latitude = (-pi/2 + 2 * atan(sqrt(X(0)*X(0) + X(1)*X(1))/(2*6371225*0.97276)))*360/(2*pi) }```

one has to declare the variables Longitude and Latitude as exported variables in a solver

```Exported Variable 1 = String "Longitude"
Exported Variable 1 DOFs = 1
Exported Variable 2 = String "Latitude"
Exported Variable 2 DOFs = 1```

and then initialize it in the initial condition for the body on which this particular solver is run

```Longitude = Variable Coordinate
Real MATC "longitude(tx)"
Latitude = Variable Coordinate
Real MATC "latitude(tx)"```

Implementation of a prognostic free surface problem

In prognostic free surface problems, the kinematic boundary condition for the free surface is a partial differential equation (PDE) and hence has to be solved. The domain of the solution, though, is a "(probelm dimension - 1)" entity. I.e., in case of a 3D glacier, the domain is a surface, in case of a flow-line model being applied, the domain is a line.
To this end, Elmer has a feature that enables a boundary to be declared as an additional body. The syntax is the following:

```Boundary Condition 2
Name = "surf"
Target Boundaries = 2
! this is the assigned body number of the free surface domain
Body ID = 2
.
.
.
End```

This declares the free surface to act as the second body in the simulation (if only one 3D body is declared). For this body, separate sections (body force, initial conditions, ...) have to be declared in a separate body declaration:

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

the different sections attached to this body then may read as:

```! NB:  units in meters MPa years
!-------------------------------
Equation 2
Active Solvers(1) = 5 !(the Solver ID of the FreeSurfaceSolver)
Flow Solution Name = String "Flow Solution"
Convection = Computed
End

Initial Condition 2
! this sets the initial position of the free surface variable
FS = Equals Coordinate 3
! and this the reference to this initial condition
ReferenceFS = Equals Coordinate 3
! this sets a zero deviation from the input mesh
Mesh Update 3 = Real 0.0
Mesh Update 2 = Real 0.0
Mesh Update 1 = Real 0.0
End

Body Force 2
! the accumulation flux in the 3 Cartesian coordinates
FS Accumulation Flux 1 = Real 0.0e0
FS Accumulation Flux 2 = Real 0.0e0
! 1.0 m/a w.e. converted into ice equivalent)
FS Accumulation Flux 3 = Real \$ 1.0 * 910.0/1000.0
End

Material 2
Density =  Real MATC "910.0*1.0E-06*(31556926.0)^(-2.0)"
! the minimum value of the free surface variable
! keeping a minimum flow depth of 10 meters
Min FS = Variable Coordinate 3, Height
Real MATC "tx(0) - tx(1) + 10.0"
! allow a maximum thickening of 100 meters with
! respect to the initial position of the free surface
Max FS = Variable ReferenceFS
Real MATC "tx(0) + 100.0"
End```

There are at least two solvers needed, namely, the free surface solver itself and the mesh update solver, in order to shift the nodes of the mesh in order to comply with the altered geometry imposed by the change in the free surface:

```! the free surface solver to be solved on Body 2
!-----------------------------------------------
Solver 5
! usually, the solver is executed only after the thermo-mechanical
! problem has obtained a solution on the time-level
Exec Solver = "After TimeStep"

Equation =  String "Free Surface Evolution"
! the name of the variable
Variable = "FS"
Variable DOFs = 1

Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"

! this enables the limitation of the free surface
! by upper and/or lower limits (see material section above)
! using a variational inequality formulation
Apply Dirichlet = Logical True

Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations  = 1000
Linear System Preconditioning = ILU1
Linear System Convergence Tolerance = 1.0e-08

Nonlinear System Max Iterations = 100 ! variational inequality needs more than one round
Nonlinear System Min Iterations = 2
Nonlinear System Convergence Tolerance = 1.0e-06

Steady State Convergence Tolerance = 1.0e-4
! switch that off in parallel runs, as it may introduce
! partition dependent relaxation factors:
! Maximum Displacement = Real 10.0

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

! this is needed if the variational inequality method
! (Apply Dirichlet = Logical True) is applied
Exported Variable 1 =  FS Residual
Exported Variable 1 DOFS = 1

! this variable contains the free surface's initial
! values (set in initial condition above and needed
! for limiting maximum changes as well as in the boundary condition
! at the free surface)
Exported Variable 2 = ReferenceFS
Exported Variable 2 DOFS = 1
End

! the mesh update solver to be solved on body 1 (the main ice volume)
!--------------------------------------------------------------------
Solver 6
Equation = "Mesh Update"
! usually, the solver is executed only after the thermo-mechanical
! problem has obtained a solution on the time-level
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-06

Nonlinear System Max Iterations = 1
Nonlinear System Convergence Tolerance = 1.0e-06
End```

The material section entry for the mesh update solver then read as (have to be inserted in the glacial body's material section)

```Youngs Modulus = 1.0
Poisson Ratio = 0.3```

If the model consists of a bedrock boundary, a free surface and a side boundary, the boundary condition section read as follows:

```! the bedrock boundary
----------------------
Boundary Condition 3
Name = "bedrock"
Target Boundaries  = 1

! zero flow height at bedrock
Height = Real 0.0

! no-slip condition
Velocity 1 = 0.0
Velocity 2 = 0.0
Velocity 3 = 0.0

! keep the nodes at the bedrock, as they are
! NB.: iso-static rebound could  be linked in here
! by linking the change of the bedrock position
! to Mesh Update 3, similar as for the free surface
! at the upper surface boundary
Mesh Update 1 = Real 0.0
Mesh Update 2 = Real 0.0
Mesh Update 3 = Real 0.0
End

! the free surface boundary
!--------------------------
Boundary Condition 2
Name = "surf"
Target Boundaries = 2

! this is the assigned body number of the free surface domain
Body ID = 2

! zero flow-depth at the free surface
Depth = Real 0.0

Mesh Update 1 = Real 0.0
Mesh Update 2 = Real 0.0

! links the free surface to the mesh update variable
Mesh Update 3 = Variable FS, ReferenceFS
Real MATC "tx(0) - tx(1)"
End

! the side wall
!--------------
Boundary Condition 3
Name = "side"
Target Boundaries  = 3
! free slip, no penetration, hydrostatic pressure
! -----------------------------------------------
Normal-Tangential Velocity = True
Velocity 1 = 0
Velocity 3 = 0
Flow Force BC = True
Normal Pressure = Variable Depth
Real MATC "tx*910.0*1.0E-06" ! in mPa

! glacier boundary remains unchanged
Mesh Update 1 = Real 0.0
Mesh Update 2 = Real 0.0
Mesh Update 3 = Real 0.0
! freeze the free surface
FS = Equals ReferenceFS
End```

Additional two instances of the FlowDepthSolver might be needed to inquire the flow depth/height of the glacier (a non-trivial problem on an unstructured FEM-mesh)