Boundary conditions

[pohlan]

The solution that the solver converges to is not reasonable. I think it has something to do with the boundary conditions. Setting https://github.com/pohlan/SheetModel.jl/blob/e2ca6f7e0bdc5920675160d96885ee175d8d9363/src/SheetModel.jl#L454-L456 does not seem to be enough. Furthermore, for the "A1" case the solution is only reasonable when the ϕ values are explicitly set outside the boundary, even when qx and qy are already set to zero across the boundary as above. When uncommenting https://github.com/pohlan/SheetModel.jl/blob/e2ca6f7e0bdc5920675160d96885ee175d8d9363/src/SheetModel.jl#L210-L212 the solutions looks ok.

(Here the algorithm adds a row/column on each boundary with zero ice thickness.)

[pohlan]

High ϕ outside of glacier (902417c)

Instead of applying the boundary conditions manually, I tried setting ϕ very high outside of the boundary to prevent outflow. There shouldn't be any inflow since h is zero where the ice thickness H is zero. https://github.com/pohlan/SheetModel.jl/blob/902417c92b79cc2dff770da3d7deb3f7a24d3836/src/SheetModel.jl#L506 ϕ_outbound is a mask for the grid points where H=0 but which are neighbours of cells where H>0. The convergence is faster (~10^4 iterations instead of ~10^5) but the result is not very meaningful, there is quite some outflow and inflow to/from the ϕ_outbound cells. If I only apply the above to ϕ[H .== 0.0] then there is no convergence at all.

[pohlan]

Yet another way to define d_eff (9601ec95590e38cf2e05c1717484b86f2704d9dc) I tried the equivalent of the iceflow code in Ludovic's EGU21 git repo https://github.com/luraess/julia-parallel-course-EGU21/blob/b700a9ad0d1d14f26e0590ce41a8e63380a1df15/scripts/iceflow.jl#L66-L72 where gradϕ and d_eff are defined on points centered between ϕ/h points. No upstream but a lot more points are considered for averaging. https://github.com/pohlan/SheetModel.jl/blob/9601ec95590e38cf2e05c1717484b86f2704d9dc/src/SheetModel.jl#L400-L403

But nothing changed in terms of convergence or meaningfulness of the solution.

Maybe d_eff needs some special treatment at the boundaries?

[pohlan]

The solution that the solver converges to is not reasonable. I think it has something to do with the boundary conditions. Setting

https://github.com/pohlan/SheetModel.jl/blob/e2ca6f7e0bdc5920675160d96885ee175d8d9363/src/SheetModel.jl#L454-L456

does not seem to be enough. Furthermore, for the "A1" case the solution is only reasonable when the ϕ values are explicitly set outside the boundary, even when qx and qy are already set to zero across the boundary as above.

Fix (df7b186): The flux boundary condition needs to be defined before d_eff is calculated, e.g. through dϕ_dx and dϕ_dy directly: https://github.com/pohlan/SheetModel.jl/blob/df7b18620b41a92d4efd9a8cfda68bc403f26851/src/SheetModel.jl#L377-L386

Then the ϕ values don't need to be set explicitly outside the boundary and qx values are zero automatically across the boundaries.

However, the valley glacier still has a problem which may or may not have something to do with the boundary, not sure.

[pohlan]

I think the fix of the previous comment actually solved the boundary issue.