Try the regularisation trick with the PT method

[pohlan]

Here's a file that runs the PT method for either of the two h equations: https://github.com/pohlan/SimpleSheet.jl/blob/891cf1bdf63ca9107fe975585e384faee0cf87fc/run_it_simple_PT.jl#L109-L112 The first one as previously done (with actual storage), the second one as in Büeler & van Pelt (storage just a regularisation), giving the following results, respectively:

To be investigated further tomorrow, maybe also with two different e_v values.

[mauro3]

Maybe you can have a look at doing both, as in PR #4? But first I'd just try with e_v large.

[mauro3]

Also, the h-timestep CFL condition can be gotten as in B&P. For the ODE it would be the same as for the PDE as they are equivalent.

[mauro3]

These are steady states, right? I would have thought that they should be indentical. Only transient states should differ.

[pohlan]

These are steady states, right? I would have thought that they should be indentical. Only transient states should differ.

Yes steady states. I have no idea why they are not the same.

.. first I'd just try with e_v large.

I've tried a range of magnitudes for e_v from 1e-3 to 1e3 but the result is the same.

Also, the h-timestep CFL condition can be gotten as in B&P. For the ODE it would be the same as for the PDE as they are equivalent.

The way the CFL condition is written now (CFL for both h and ϕ, but slightly different), it doesn't work with neither equation of h.

[pohlan]

However, coming back to these lines

https://github.com/pohlan/SimpleSheet.jl/blob/ab5ff1df935f3de5a753d64be111bc43fd90c6a4/run_it_simple_PT.jl#L113-L119

we see that div_q is cancelled out in any case if ev_num=0 and if we don't set dϕdt[1, :] .= 0. before. ev_num was mainly introduced for the explicit method and in the PT case we don't need it anyway so it shouldn't be a problem to leave it away.

[pohlan]

This works but then we're exactly where we were in the beginning. The idea here was to use regularisation also for the PT method but if we remove e_v_num and we want to regularise with e_v we have to remove dϕdt * e_v from the mass conservation equation so that we have:

dhdt .= .- div_q .+ Λ 

so we're again at the same point that we cannot calculate div_q at the Dirichlet boundaries.

[pohlan]

Last note on the regularisation: I now formulated it like this in all the scripts:

https://github.com/pohlan/SimpleSheet.jl/blob/66f4379a45f0154d6cfc6c57008dc910a6cb4c2b/diffeq_explicit.jl#L122-L125

So there is no use_masscons_for_h option anymore as it was equivalent to the other equation anyway. The only option to choose is whether there is regularisation (e_v_num>0) or not (e_v_num=0). To make it work, it is important that

1. Use the formulation for h that doesn't contain div_q as this one cannot be computed at the boundaries, i.e.

dhdt = Σ * vo - Γ * vc + e_v_num * dϕdt

instead of the previous formulation for use_masscons_for_h=true:

dhdt =  - div_q + Λ - e_v * dϕdt

Deriving those For ϕ we solve

(e_v + e_v_num) * dϕdt + div_q + Σ * vo - Γ * vc    = Λ
<=>
dϕdt = (- div_q - (Σ * vo - Γ * vc) + Λ) / (e_v + e_v_num)

and for h we solve the 'true' mass conservation equation

e_v * dϕdt + div_q + dhdt  = Λ

So to bring the equation of ϕ into the form for h we need to define dhdt as

dhdt = Σ * vo - Γ * vc + e_v_num * dϕdt

Solving the 'true' mass conservation equation for dhdt directly gives

dhdt =  - div_q + Λ - e_v * dϕdt

2. Another important point is that dϕdt is set to zero at the Dirichlet boundaries before setting dhdt += e_v_num * dϕdt. This way it is not required to set boundary conditions for h. This was wrong in a few scripts, which I fixed in 66f4379a45f0154d6cfc6c57008dc910a6cb4c2b.