"Valley" geometry not converging

[pohlan]

https://github.com/pohlan/SheetModel.jl/blob/bdba50356dffafef9a2f0e60e1d06afe25234e0c/examples/run_SHMIP.jl#L1-L10

Same solver and +/- same parameters as for "sqrt" geometry, but not converging. Tried a lot of different damping parameters/dtau scaling factors, as well as both error definitions, doesn't change much. Maybe something's wrong with the boundary conditions?

[luraess]

Had a quick try, seems indeed not to converge super well.

• Maybe worth checking the BCs and also initial condition.
• You could also plot the residual or update rates to see where things go wrong (i.e. lack of convergence in the entire domain, or in specific regions).
• If nothing helps, reverse-engineering could give more insights (removing temporarily nonlinearities, etc, until it works and figuring out the trouble-maker).

[pohlan]

Didn't find the issue so far. Now trying to change some features of "sqrt" towards "valley"

First try: Dirichlet b.c. (151c61c3b1d828bcbb60dedcb65530b285a836f4) Impose the Dirichlet b.c. only to a few x=0 grid points instead of the whole ϕ[1, :] column (for the "sqrt" geometry). Result: Even if I only leave away the corner points, the solver does not converge anymore.

https://github.com/pohlan/SheetModel.jl/blob/151c61c3b1d828bcbb60dedcb65530b285a836f4/src/SheetModel.jl#L223-L225

However, if I use the same procedure for the E1 case (pw=0 at all lower x-boundaries) it still does not converge.

[pohlan]

Playing around with dtau (cd4a834a94a35110b921f14493431f9849a4d795) I noticed that the jump in ϕ is always pretty large between ϕ[1, :] and ϕ[2, :]:

Now I randomly tried to decrease the dtau in the second column: https://github.com/pohlan/SheetModel.jl/blob/cd4a834a94a35110b921f14493431f9849a4d795/src/SheetModel.jl#L431

Leading to a ϕ distribution which looks more reasonable:

Unfortunately it does not have much impact on the error_ϕ but maybe this is something to think about? Should there be a dependency of dϕ_dx or dϕ_dy in the definition of dtau?

[pohlan]

Next try: modifying "sqrt" bed topography (ad18be325814880de63a118e029504714f45ea4d) I left everything the same for the A1 SHMIP case except the bed topography which I changed from zero everywhere to a parabolic shape in y-direction but keeping it zero at x=0: https://github.com/pohlan/SheetModel.jl/blob/ad18be325814880de63a118e029504714f45ea4d/examples/SHMIP_cases.jl#L84-L85 This change already caused the solver to not converge. Apparantly the upstream scheme was not able to handle this. If it already doesn't work here, then it's not surprising that it doesn't work for the "valley" geometry either. Cause and solution are still to be found.

[mauro3]

Are you sure the upstream scheme is correct?

[luraess]

I would also check the upstream scheme. Maybe give it a try on a simple 1D case ? Or add also and centred scheme such that you take the local average to get correct staggering of your diffusion coefficient values.

[pohlan]

Next try: setting β=2 Then the flux is a linear function of dϕ_dx or dϕ_dy: https://github.com/pohlan/SheetModel.jl/blob/bdba50356dffafef9a2f0e60e1d06afe25234e0c/src/SheetModel.jl#L238-L241 With this, the valley geometry cases converge and also the sqrt geometry with modified bed topography (see comment above).

[mauro3]

That's good news!

There is also a bug in the flow relation: abs(d\phi_du) should be sqrt(dphi_dx^2+dphi_dy^2). (although, what you solved should also converge)

[pohlan]

There is also a bug in the flow relation: abs(d\phi_du) should be sqrt(dphi_dx^2+dphi_dy^2). (although, what you solved should also converge)

Is fixed now (760433d61196342bb6e359ea9f9c277794fa51b8). E.g. for qx I took the average of the two upstream dphi_dy since they're not on the same grid. (In this piece of code i is already determined as the upstream index:)

https://github.com/pohlan/SheetModel.jl/blob/760433d61196342bb6e359ea9f9c277794fa51b8/src/SheetModel.jl#L386-L393

This fix does not change anything about the convergence though.

[pohlan]

Only solving for ϕ and commenting out the h solver (6459a631a6a848bc92f8a370218a4f8fd1b08e33) also only works for β=2

[pohlan]

Major fix (1dd6926e9482c712de195434962c2753d3723171) The effective diffusivity did not include the |∇ϕ|^(β-2) term, actually I already stumbled upon it somehow:

Should there be a dependency of dϕ_dx or dϕ_dy in the definition of dtau?

Now I define the effective diffusivity as https://github.com/pohlan/SheetModel.jl/blob/1dd6926e9482c712de195434962c2753d3723171/src/SheetModel.jl#L459 instead of k * h^α. Then the modified sqrt geometry (bed changing in y direction) converges, even though quite slowly (~10^4 iterations). For the valley geometry it's still too slow to converge in a reasonable time. Maybe it needs better definitions of the PT-parameters.

[pohlan]

e2ca6f7 Converging now (even though quite slowly) with different definition of d_eff, in particular |∇ϕ| is calculated by averaging the derivatives in x and y directions to the grid points where h is defined.

However, the solution is still not correct, which is probably due to the boundary conditions (see #7)