## Rotated coordinate system for basal boundary condition

Extension of Elmer in computational glaciology

### Rotated coordinate system for basal boundary condition

Hi,

I am tracking down some issues in our code, and I am seeking some clarification on the way that Elmer handles the non-penetrating condition at the basal boundary condition (i.e. the ice cannot penetrate through the bedrock) in the Navier-Stokes solver.

As I understand it, Elmer rotates the coordinate system of the nodes that are on the bottom boundary so that the x-coordinate is normal to the base. If you set the parameter "Back Rotate N-T Solution = Logical True", it is supposed to rotate the coordinate system back to the original Cartesian coordinate system after DefaultSolve (and by extension SolveSystem) is finished executing. It is not clear, however, if the coordinate system of the input parameters prior calling DefaultSolve is supposed to be bed-normal coordinate system or in Cartesian coordinates (regardless of the value of Back Rotate N-T Solution). The reason I ask, is because I am setting the results of the Shallow Ice Approximation as an initial guess for the velocity of the ice, and this is stored in Cartesian coordinates in our program.

Another question is why it was set up like this. Wouldn't it be easier to rotate the non-penetrating condition to Cartesian coordinates, rather than rotating the nodes to a local bed-normal coordinate system? Then there would only have to be one rotation calculation at the start of the simulation and the coordinate system of all the nodes would remain consistent.

Cheers,
Evan
Evan

Posts: 12
Joined: 12 Aug 2014, 14:42

### Re: Rotated coordinate system for basal boundary condition

Just to follow up, after investigating the problem I am having with our code, it seems that the Back-Rotate routine is rotating not only the bottom layer of nodes, but also the layer directly above it. However, after exit from DefaultSolve, the second layer is not rotated back. Below is a stack of nodes at a single location (the mesh has been extruded to six layers) before and after calling DefaultSolve. As you can see, the velocity values in the second layer are incorrect.

Code: Select all
`vx (before)      vy(before)         vz(before)               vx(after)      vy(after)        vz(after)-669.15             -665.18              -5.46               -507.85        -628.50           -8.29-672.47             -668.10             -13.10               -64.700        -625.80           575.94-675.32             -670.22             -20.40               -574.88        -633.43          -61.43-677.63             -671.53             -26.88               -578.83        -634.89          -64.48-679.35             -672.06             -31.61               -582.10        -635.60          -67.81-680.40             -671.94             -33.28               -584.74        -635.68          -69.17`
Evan

Posts: 12
Joined: 12 Aug 2014, 14:42

### Re: Rotated coordinate system for basal boundary condition

Hi Even

I try to clarify the 2nd question. The N-T coordinate system does not rotate nodes. Instead it rotates the velocity vector related to the nodes on the boundary. The nodes on the N-T coordinate system should be marked so that the operation is performed for only them. When you rotate the contribution to the matrix equation to the N-T directions you can also set the Dirichlet conditions accordingly in the N-T directions. This is convenient since we know that normal component should be zero and the tangential component can remain unknown. If you would not do the rotation you cannot set either component using Dirichlet conditions since both include part if the unknown tangential velocity. If the tangential velocity would be known it would be easy to prescribe the velocity for all components in the cartesian system.

-Peter
raback

Posts: 3211
Joined: 22 Aug 2009, 11:57
Location: Espoo, Finland

### Re: Rotated coordinate system for basal boundary condition

Hi Peter,

Thank you for the explanation.

It turns out the problem I am having is not a problem with the rotation itself, but that we are using a reduced A matrix for DefaultSolve after the first time step. Somewhere in our code we have missed changing a pointer, so it is rotating too many nodes. Because the A matrix has dimensions that are about half what it was in the first time step, it looked like it was rotating two layers instead of one.
Evan

Posts: 12
Joined: 12 Aug 2014, 14:42