Climatemachine.jl: 'Float32' issue in AtmosDycore

[jkozdon]

It's possible there is an issue when the AtmosDycore is run in Float32 mode. Results for rising_thermal_bubble.jl give weird results (visualize the velocity field)

edit: thinking there may be an issue with the background state being included in the solution. Plan is to try to subtract it off to see if that helps.

h/t @smarras79

[jkozdon]

This may be related to the large background state, so we may want to have reference and perturbation state.

[smarras79]

The issue is the same with Float64. I am testing the later revision: 5f1226ac5b054fee5dc4cee0e2ecfc5997fb171d

[smarras79]

This issue should now be closed

The issue was simply related to the fact that viscosity is zero by default

[jkozdon]

I think that there is still an issue we should look at with the background field. Without viscosity the code will eventually blow up, but the stuff you pointed to was before that happens. (For me, Float64 stuff looks smooth at early time.)

I'm going to leave this open until we explore the removing of the background field.

[smarras79]

perfect.

[ali-ramadhan]

Out of curiosity, is the issue with the AtmosDycore or is Float32 just not enough precision?

I know for a finite volume ocean, Float32 works fine for the time stepping (which is all local operations). But for the global nonhydrostatic pressure solve/correction step which involves FFTs and DCTs, Float32 is not precise enough and we get weird artifacts (https://github.com/ali-ramadhan/Oceananigans.jl/issues/55) that build up over 1000's of iterations and things blow up some time around time step number 1,000-10,000+.

It may just be because we're solving the incompressible Boussinesq equations: since the continuity equation is satisfied only to eps(Float32) = 1.1920929f-7, mass eventually accumulates quickly enough to blow up the model. But with eps(Float64) = 2.220446049250313e-16 this isn't a numerical issue anymore.

[jkozdon]

Float32 has been shown to be precise enough for many PDE simulation, so we do not know. We have perturbations on top of background states that are several orders of magnitude smaller than the background state, so we are hoping that we can separate out the two states and still run Float32. The advantage being the code will be faster and we can run bigger problems on the same hardware.

There may certainly be places that we can and cannot get away with Float32, but it will certainly take a bit of care to make sure we don't get major cancellation errors.

When there are linear solves involved, even more care will certainly be involved.

Properly scaling the variables may certainly be key.

[jkozdon]

old dycore is going away...