Explicit and RK4 scripts to reproduce diffeq results

[luraess]

Updating the run_it_simple.jl script to reproduce the diffeq script results. Works - as long as one set e_v = 1.

Added run_it_simple_RK4.jl, the RK4 version of the run_it_simple.jl for comparison. Both seem as stable (or as unstable if run with stiff params ...)

[mauro3]

Also works with e_v=1e-3 but with dt=1s. Similar speed to DiffEq's DP5 solver (a RK 4-5 method). ROCK4 is about 2x faster.

[luraess]

Does it also run with run_it_simple.jl? I had the feeling RK4 does not help much here versus the simple forward Euler.

[mauro3]

Yes, that runs too and also faster (10s vs 60s (RK4) vs 20s (DiffEq ROCK4)). But it get oscillations at the top with dt=1.0s. With 0.5s it looks ok though (23s; i.e. at the same speed as ROCK4 but I suspect that ROCK4 is much more accurate).

[luraess]

Just pushed more accurate timing of the simple solvers (RK4 and the other version that actually has Crank-Nicolson stepper CN=0.5). Now the RK4 runs in 20sec which is the same as the diffeq ROC4. The CN code runs in only 5sec but yes, produces some wiggles with dt = 1s. See new codes here https://github.com/pohlan/SimpleSheet.jl/commit/6dfc088391439652a8d69acfa6fbfccad9e94fea .

Note: The RK4 could be boosted straight fwd with ParallelStencil and GPU 🙁

[mauro3]

Ok, but all in all the explicit methods do need a rather short time step. In cases where there is physics going on on short time scales, say diurnal melt input, this maybe ok (implict solvers there need timesteps of ~1h in the FEM GlaDS implementation of Elmer). For more steady system, say Antartica, this maybe less the case.

I think the only place where adaptive time stepping maybe useful for explicit methods is to iron-out initial transients (although that can be hand-coded easily too). Afterwards the CLF-constraint will be much lower than any physical constraint on dt.

I wonder whether we could use both a "physical e_v" and a "regularization e_v"; however, that would necessitate to solve the conservation equation for h as B&P do. In a similar vein, could a regularization be used within the pseudo-transient iteration? Maybe a large e_v at the beginning of PT which then shrinks with iteration count.

[mauro3]

Note: The RK4 could be boosted straight fwd with ParallelStencil and GPU slightly_frowning_face

That ROCK4 solver also looks interesting. It is kinda made for what you like to do: PDEs without Newton (and should thus work easily on the GPU). Not sure how tricky it would be to program https://doi.org/10.1137/S1064827500379549

[luraess]

I guess e_v will indeed be the key for making steady state converge fast. Note that in curent implementation, we have e_v * dϕ/dt = dϕdt. Thus, when dt->inf and dϕ/dt->0, it does no longer matter what e_v is. I guess one could use it straight-fwd as numerical param for steady state. For hybrid, I have to think. Relaxation or similar could do the job. Or adding a second e_v_num param...

[luraess]

Not sure how tricky it would be to program https://doi.org/10.1137/S1064827500379549

Thanks! However, I don't get what's in the paper. Anyway, the major issue with RK4 and other higher-order methods is the amount of info on has to save (RK4 needs 6 more arrays per DoF...) which is not optimal if memory bound...

[mauro3]

Note that when doing also R-channels, then having a large e_v can lead to instabilities (a jökulhlaup type instability). For ice sheets those start somewhere ~1e-2 to 1e-3. I'm not sure whether the same applies when just using it as regularization, maybe not.

But we shall cross that bridge when we get to it.

[mauro3]

Regarding the Crank-Nicolson: that method is an implicit method but what you coded is an explicit multi-step method, no?