Similar results from other solvers

[HannoSpreeuw]

I wanted to share that different solvers, i.e. different from "Euler" give similar results, at least wrt oscillations.
In rhythmite one can pick "Euler", "RK23", "RK45" or "DOP853".
The latter 3 call Scipy's solve_ivp and by default do not apply clamping of the rhs values.
This is with "DOP853" at 50 cm depth over 250 ka.

[HannoSpreeuw]

Similarly, in rhythmite one can choose an implicit solver, like "BDF" and you find a similar result.
This uses a numerically approximated Jacobian.

[HannoSpreeuw]

As @cedrict suggested at one of our meetings, the oscillations can possibly be associated with a change in sign of U.
Usually U is positive (downward) and backwards differencing adds stability to the integration.
When U changes in sign (becomes negative, i.e. upward motion) forward differencing should add stability to the integration.

This is accommodated for in rhythmite through this line:

                dAR_dt[i]= -u[i]*( (AR[i+1]-AR[i-1])/(2*h) - self.upwind_switch*np.sign(u[i])*h/2*(AR[i+1] - 2*AR[i] + AR[i-1])/h**2 )\
                    + self.R_AR(AR[i],CA[i],c_ca[i],c_co[i], x[i])    

and this line:

                dCA_dt[i]= -u[i]*( ( CA[i+1]-CA[i-1])/(2*h) - self.upwind_switch*np.sign(u[i])*h/2*(CA[i+1] - 2*CA[i] + CA[i-1])/h**2 )\
                    +self.R_CA(AR[i],CA[i],c_ca[i],c_co[i],x[i])

This works nicely, until U becomes negative (i.e. upward motion). Then you need a bottom boundary condition for aragonite and calcite, but those are not available. The compromise chosen in rhythmite is to still use backward differentiation at the deepest node despite a negative U through this line:

            elif (i==len(x)-1):
                dAR_dt[i]= -u[i]*(AR[i]-AR[i-1])/h + self.R_AR(AR[i],CA[i],c_ca[i],c_co[i],x[i])

and this line:

            elif (i==len(x)-1):
                dCA_dt[i]= -u[i]*( CA[i]-CA[i-1] ) / h + self.R_CA(AR[i], CA[i], c_ca[i], c_co[i], x[i])

If were to mimic marlpde where the backward derivative is always applied independent of the sign of U, for the reason described above i.e. because of missing boundary conditions we would have to adjust the rhythmite code like this, i.e. just remove sign(u[i]:

                           dAR_dt[i]= -u[i]*( (AR[i+1]-AR[i-1])/(2*h) - self.upwind_switch*h/2*(AR[i+1] - 2*AR[i] + AR[i-1])/h**2 )\
                    + self.R_AR(AR[i],CA[i],c_ca[i],c_co[i], x[i])
....
                dCA_dt[i]= -u[i]*( ( CA[i+1]-CA[i-1])/(2*h) - self.upwind_switch*h/2*(CA[i+1] - 2*CA[i] + CA[i-1])/h**2 )\
                    +self.R_CA(AR[i],CA[i],c_ca[i],c_co[i],x[i])
                    

Applying these two code modifications makes the oscillations disappear:

[HannoSpreeuw]

Oops. I thought I plotted for 250 ka, but it seems only half of that time, sorry.

Well yeah, taking out the sign(U) will only make a difference at a later stage in the evolution, when U turns negative (=upward).
Until that stage the upwind solution is still affected.

What do you mean by "without having the solution defined outside"?

[csummers25]

If the velocity is negative, we have material flowing in the boundary at the bottom, except we don't know what should be flowing in since we don't know what the solution is outside the bottom boundary. By doing a backwards derivative we avoid the problem of trying to use values from nodes that don't exist but we basically assume that the downwind derivative (the one between N-1 and N) is a reasonable approximation to the derivative at the final node (for upwind, N and N+1) when there's no guarantee that this is true.

[HannoSpreeuw]

Yes and the backward derivative can also give stability in the integration unless U changes sign, then you need to switch to a forward derivative - and fix the missing bottom boundary value problem e.g. in the way you implemented.

It turns out that, with soln_full being the output from solve_ivp, soln_full.t does not match all the values in t_eval although t_eval=t_eval was included in the solve_ivp arguments. I need to sort out why. This explains why the plot was not covering the full 250 ka.

[HannoSpreeuw]

Unfortunately, I cannot run rhythmite with an implicit solver and collect all the requested t_eval from the returned solution. No times beyond some 54% of the requested 250 ka are returned.

For an explicit "RK45" solver it crashes with ZeroDivisionError at that 54% of 250 ka.

To clarify again the perspective of this investigation: can I expand the plot above for the full 250 ka? That plot was produced with np.sign(u[i]) effectively replaced by 1, such that upwind solutions are applied until U changes sign and becomes negative. I wanted to find out if that kills oscillations.

Conclusion: it probably does not, but I cannot verify that. Sorry for the fuss.

[EmiliaJarochowska]

If the velocity is negative, we have material flowing in the boundary at the bottom, except we don't know what should be flowing in since we don't know what the solution is outside the bottom boundary. By doing a backwards derivative we avoid the problem of trying to use values from nodes that don't exist but we basically assume that the downwind derivative (the one between N-1 and N) is a reasonable approximation to the derivative at the final node (for upwind, N and N+1) when there's no guarantee that this is true.

We mostly use this (questionable) approach because it is used in Fortran and because we have nothing better to offer.

At the same time, it seems that this bottom boundary condition is essential for the oscillations. So in order to "rescue" this model we might revise the idea of adding a "storage" matrix to import solids back into the system. That is a bigger operation and the reward is uncertain, but it seems it would be the only "honest" way of handling this problem.

[cedrict]

thank you @HannoSpreeuw for carefully checking this code against yours and L'Heureux!

[EmiliaJarochowska]

I am working on a plan of a manuscript and working through a series of runs we should document. ⁠
@HannoSpreeuw, last Thursday we discussed adding a record of your runs to the table Charlotte and I have been using on Research Drive. I am thinking of re-doing a set of runs akin to what Charlotte has done with varying porosity diffusion coefficients and sizes of the grid.
For that purpose I would like to

• agree on what should be in the main branch of marlpde in terms of bottom boundary condition so we can carry out these runs.
• decide whether this should be shown only with the Euler scheme with a mention of the results when other solvers are used. Should we explicitly include the comparison with other solvers? I see it might not be always possible without further tweaks to rhytmite, but I am not convinced they are worth it.

[HannoSpreeuw]

Wrt the second bullet: I would disregard Eulerian schemes, since they cause problems for rhythmite: you essentially cannot run an Eulerian scheme in rhythmite without clamping. And we want to avoid clamping, I reckon.
All the other solvers, i.e. all solvers from solve_ivp can be invoked from both rhythmite and marlpde without clamping.
This includes the implicit solvers "BDF", "LSODA" and "Radau" with numerically approximated Jacobians.

I will comment on the first bullet tomorrow.

[csummers25]

I think the Eulerian scheme with clamping is fine, since it produces similar results to the other solvers (as you show above) and to the other codes. It is also stable over longer timescales, whereas the solve_ivp versions are not, so if we want to look at the long term behaviour of the larger grid size models I think this is the only way to do it.

[HannoSpreeuw]

I did not yet encounter the stability issue for other solvers over longer timescales.
Do you have an example for a derailed run at hand?

Any referee will probably ask questions about clamping, since it is a dubious measure.

[EmiliaJarochowska]

We are comparing with Fortran, which also uses clamping, so it makes sense to keep it similar between implementations. Perhaps we could report on the effects of clamping for one of the implementations, for example marlpde, to assess their importance. But they don’t seem essential for the research question, since the existence of oscillations doesn’t depend on clamping.

…

________________________________ From: Hanno Spreeuw ***@***.***> Sent: Tuesday, July 23, 2024 4:39:32 PM To: cedrict/rhythmite ***@***.***> Cc: Jarochowska, E.B. (Emilia) ***@***.***>; Comment ***@***.***> Subject: Re: [cedrict/rhythmite] Similar results from other solvers (Discussion #4) CAUTION: This email originated from outside of Utrecht University. Do not click links or open attachments unless you recognize the sender and know the content is safe. I did not yet encounter the stability issue for other solvers over longer timescales. Do you have an example for a derailed run at hand? Any referee will probably ask questions about clamping, since it is a dubious measure. — Reply to this email directly, view it on GitHub<#4 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AELMSVAM722EHL3GGOPVIQTZNZTKJAVCNFSM6AAAAABKYY4AHWVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTAMJSG43DENA>. You are receiving this because you commented.Message ID: ***@***.***>