Salinity drifting when there is no salinity forcing

[fabien-roquet]

When there is no salinity forcing, salinity is still drifting.
We believe it happens with the implicit diffusion scheme which produces noise which accumulates over time, but we are not sure exactly why. Because salinity is used in the equation of state, that creates disturbances on the dynamics that will increase over time.
It would be best to enforce/control salinity conservation. Some similar bug might be affecting temperature as well.

[dionhaefner]

Can you share more details about your setup? How can I reproduce this?

[Titouan-Moulin]

neverworld2.zip
We use this neverworld2 setup inspired from the acc setup and we ran it for 1d before plotting a salinity profile

ds_snap = xr.open_dataset(neverworld2.snapshot.nc)
SA = (ds_snap["salt"].isel(Time=0)-35)
SA.isel(xt=30).plot()

[dionhaefner]

Some numerical noise is certainly expected. Your reproducer outputs salinity fluctuations on the order of 10^-14 which is well below the noise threshold and should not lead to significant model drift.

[fabien-roquet]

Normal maybe but problematic for salt and heat conservation...
Is there a way to keep it under control?

[dionhaefner]

Say more how that is problematic? Surely the model is running for << than 1e14 steps, so that numerical noise should never reach a scale where it influences results.

[Titouan-Moulin]

We try to use a no forcing scheme : no focing salinity nor temperature and no wind and we try to investigate for the issue.
It seems to come from the solver in the function vertmix_tempsalt called in thermodynamics.py : utilities.solve_implicit. It is a call to 'core.operators.solve_tridiagonal' which use lapack.dgtsv. Maybe the issue reside elsewhere in vertmix_tempsalt than in the solver, but vertmix_tempsalt introduce noise that seem to amplify over the time.

We try to comment the lines :

    with state.timers["vmix"]:
        vs.update(vertmix_tempsalt(state))

to compare with our run. Here is a plot of the drifting, in salinity and temperature for ~2000 timesteps ~80 days (timestep are 1h), we see a clear difference if vmix is computed or not :

Here is a second plot showing salinity profile at the equator, the noise appear after the first iteration and continue to amplify :

[dionhaefner]

Can you please share a standalone script that reproduces the issue? In that case I'm happy to take a look at the numerical noise introduced by vertmix_tempsalt specifically.

[Titouan-Moulin]

What do you mean by a standalone script ? Is the setup file sufficient or do you ask for a script with the plots as well ?

[dionhaefner]

Ideally both (so I can verify if the issue is fixed or not).

[Titouan-Moulin]

no_forcing_share.zip
I ran run_veros_manually.py in debug mode which uses the neverworld2_noforcing.py setup file and I pickle data at several timesteps. (I can't give you all the extracted files because it is too heavy.)
I converted the data to xrarray (dat.nc and data_wo_vmix.nc)
The plots can be found in the .ipynb file

[dionhaefner]

I see, thanks for sharing.

It's not obvious to me whether this is because of finite numerical precision in the tridiagonal solve or because the linear system isn't set up in a way that is conservative or because your cell volume computation is off (either of which is plausible to me).

Unfortunately I don't have cycles to do the math on this right now. If someone wants to dive into this, here are some things I would check:

• Compute residuals after each linear solve and compare to observed drift.
• Write down the discretized conditions for tracer conservation (i.e., we want to preserve the integral T dV) and compare to what's implemented.

If you want a workaround, I'd disable TKE so this particular mixing scheme becomes a no-op (which doesn't mean that there won't be other sources of slight leakage due to finite precision once you add dynamics to the model).

I also still challenge the notion that this is a problem that needs fixing, a drift of O(1e-12) is totally acceptable (perhaps even inevitable) in my book.

[Titouan-Moulin]

It turns out that I was running my long_time runs in simple precision, float32 and the drifting was then several order of magnitudes higher (~1e-6) but we didn't thought to investigate that sooner
Running with double-precision with float64 seems to remove the problem. As you said, the 1e-14 drift is acceptable.
Thanks for your help