Output_dt in leapfrog time stepping and some adjustments for using Dates

[maximilian-gelbrecht]

This implements #416 and part of #417.

I made Δt_sec a Second as well, so it can be only full seconds. This also means that the adjusted Δt_at_T31 will be a bit smaller and the integration will take a few more steps. If we would want to implement it like Milan wrote in #416 then we would need Δt_sec as a Float or Millisecond. I am not opposed to that, but I thought in the initial draft I do everything in SI.

So with this I already implemented a part of #417 as well, but I did not adjust everything. After this is done, we definitely also need to go over the documentation. I already made some small corrections based on this PR, but no proper rewrite of the respective sections.

I am on vacation next week, so don't expect any replies from me. Feel free to modify this PR/branch.

[milankl]

Sorry I thought the following was a good idea (and therefore referenced it in the comments above) but it turns out it isn't. So forget about it, but for the record:

To avoid a*time_scale.value or a*value(time_scale) (and similar with /) we could also implement for SpeedyWeather the following promotion

julia> Base.:(*)(x::Second,y::Number) = x.value*y

julia> 5.5*a
ERROR: MethodError: *(::Second, ::Float64) is ambiguous.

Candidates:
  *(x::P, y::Real) where P<:Period
    @ Dates ~/.julia/juliaup/julia-1.9.4+0.x64.apple.darwin14/share/julia/stdlib/v1.9/Dates/src/periods.jl:90
  *(x::Second, y::Number)
    @ Main REPL[12]:1

which however hits ambiguity errors which don't seem to be easy to solve and therefore implement safely. I wanted to address this problem:

julia> a = Second(5)
5 seconds

julia> 5a
25 seconds

julia> 5.0*a
25 seconds

julia> 5.5*a
ERROR: InexactError: Int64(27.5)

[milankl]

@maximilian-gelbrecht I realised that this is actually a harder problem than I thought it would be. To illustrate:

Say you want output every 21601s (6 hours + 1 second) then the time step dt times number of time steps till output n needs to be n*dt = 21601. This doesn't work because 21601 is a prime number. So you need (a) either to output less than every 21601s (what we had before) or (b) you need a timestep more precise than seconds or (c) you need to restrict the possible output frequencies to something less precise than seconds.

So with the last commits I've used the chance to: (b) change the time integration time step to milliseconds precision. There's a Δt_millisec::Millisecond now that counts the clock up (=integer arithmetic no rounding errors yay!) as well as the corresponding Δt_sec::NF and Δt::NF which are the(unscaled/radius-scaled) float representations of the time step to be used in the dynamical core/parameterizations.

Some examples

julia> spectral_grid = SpectralGrid(trunc=42)
julia> time_stepping = Leapfrog(spectral_grid,adjust_with_output=true)
Leapfrog{Float32} <: TimeStepper
├ trunc::Int64 = 42
├ Δt_at_T31::Second = 1800 seconds
├ radius::Float32 = 6.371e6
├ adjust_with_output::Bool = true
├ robert_filter::Float32 = 0.05
├ williams_filter::Float32 = 0.53
├ Δt_millisec::Dates.Millisecond = 1350000 milliseconds
├ Δt_sec::Float32 = 1350.0
└ Δt::Float32 = 0.00021189766

which uses the DEFAULT_OUTPUT_DT (=21600s/6hrs), now lets be nasty and change to prime number of seconds output

julia> output = OutputWriter(spectral_grid,PrimitiveWet,output_dt=Second(21601))
julia> model = PrimitiveWetModel(;spectral_grid,time_stepping,output)
julia> simulation = initialize!(model)
[ Info: Adjust with output: Time step shortened to 1080050 milliseconds (-20%)

The output dt is in milliseconds 21601000, which is only about 16 time steps à 1350000 ms, given the divisors of 21601000

julia> Primes.factor(21601000)
2^3 * 5^3 * 21601

i.e. 1,2,4,5,8,10,20,25,40,... we would need to make it 20 time steps each 1080050ms to hit those 21601s exactly even with millisecond time step precision. That's a time step 20% shorter than before hence the [ Info: from above. I've also included in the initilisation of the output a check to throw an info as follows (here triggered with adjust_with_output=false)

julia> time_stepping = Leapfrog(spectral_grid,adjust_with_output=false)
julia> output = OutputWriter(spectral_grid,PrimitiveWet,output_dt=Second(21601));
julia> model = PrimitiveWetModel(;spectral_grid,time_stepping,output)
julia> simulation = initialize!(model)
[ Info: 16 steps of Δt = 1339535ms yield output every 21432560ms (=21432.56s), but output_dt = 21601s

[milankl]

okay next problem 😆

With timesteps in milliseconds but netcdf output in integer seconds the actual timestep written to file is rounded easily leading to an alternative roundup/down in the netcdf time axis (outputting every time step with OutputWriter(...,output_dt=Second(0)))

julia> ds["time"][:]
66-element Vector{DateTime}:
 2000-01-01T00:00:00
 2000-01-01T00:22:20
 2000-01-01T00:44:39
 2000-01-01T01:06:59
 2000-01-01T01:29:18
 ⋮
 
ulia> diff(ds["time"][:])
65-element Vector{Millisecond}:
 1340000 milliseconds
 1339000 milliseconds
 1340000 milliseconds
 1339000 milliseconds
 1340000 milliseconds
 1339000 milliseconds
 1340000 milliseconds
 ⋮

I believe we could just change it to hours since X and output in Float64?! Testing ...

[milankl]

Okay that looks good now!

julia> ds["time"][:]
66-element Vector{DateTime}:
 2000-01-02T00:11:09.775
 2000-01-02T00:33:29.310
 2000-01-02T00:55:48.845
 2000-01-02T01:18:08.380
 2000-01-02T01:40:27.915
 ⋮

julia> diff(ds["time"][:])
65-element Vector{Millisecond}:
 1339535 milliseconds
 1339535 milliseconds
 1339535 milliseconds
 1339535 milliseconds
 1339535 milliseconds
 1339535 milliseconds
 1339535 milliseconds
 ⋮

[maximilian-gelbrecht]

Sorry, for no responses so far, just got back from vacation yesterday. Got teaching and some other stuff to take care of first, but I can take a detailed look at it on Friday

[milankl]

I think it's all good now, just some last tweaks for the documentation (and the Documenter is being naughty) I'll merge beforehand, but if something is not right and you suggest something different we can fix that later. I'm quite confident that you like those changes ;)

[milankl]

Many thanks for initiating this in the first place Max!