Geometry related issues
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)