Merge pull request #418 from SpeedyWeather/mg/output-dt-dates

Output_dt in leapfrog time stepping and some adjustments for using Dates

docs/Project.toml

@@ -5,7 +5,7 @@ PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
-Documenter = "0.26, 0.27"
-NCDatasets = "0.12"
+Documenter = "0.26, 0.27, 1"
+NCDatasets = "0.12, 0.13"
 UnicodePlots = "^3.3"
 PythonPlot = "1"

docs/make.jl

@@ -1,12 +1,13 @@
 using Documenter, SpeedyWeather
 
 makedocs(
-     format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing)=="true",
-                            ansicolor=true),
-    sitename="SpeedyWeather.jl",
-    authors="M Klöwer and SpeedyWeather contributors",
-    modules=[SpeedyWeather],
-    pages=["Home"=>"index.md",
+    format = Documenter.HTML(prettyurls=get(ENV, "CI", nothing)=="true",
+                            ansicolor=true,
+                            size_threshold = 600_000),      # in bytes
+    sitename = "SpeedyWeather.jl",
+    authors = "M Klöwer and SpeedyWeather contributors",
+    modules = [SpeedyWeather],
+    pages = ["Home"=>"index.md",
             "Installation"=>"installation.md",
             "How to run SpeedyWeather.jl"=>"how_to_run_speedy.md",
             "Model setups"=>"setups.md",

docs/src/how_to_run_speedy.md

@@ -112,7 +112,7 @@ Meaning that if you want to have a shorter or longer time step you can create a
 `SpectralGrid` as first argument.
 ```@example howto
 spectral_grid = SpectralGrid(trunc=63,nlev=1)
-time_stepping = Leapfrog(spectral_grid,Δt_at_T31=15)
+time_stepping = Leapfrog(spectral_grid,Δt_at_T31=Minute(15))
 ```
 The actual time step at the given resolution (here T63) is then `Δt_sec`, there's
 also `Δt` which is a scaled time step used internally, because SpeedyWeather.jl
@@ -178,9 +178,9 @@ and we have initialized the `ShallowWaterModel` we have defined earlier.
 After this step you can continue to tweak your model setup but note that
 some model components are immutable, or that your changes may not be
 propagated to other model components that rely on it. But you can, for
-example, change the output time step (in hours) like so
+example, change the output time step like so
 ```@example howto
-simulation.model.output.output_dt = 1
+simulation.model.output.output_dt = Second(3600)
 ```
 Now, if there's output, it will be every hour. Furthermore the initial
 conditions can be set with the `initial_conditions` model component

docs/src/output.md

@@ -8,64 +8,83 @@ There are many options to this available.
 
 The output writer is a component of every Model, i.e. `BarotropicModel`, `ShallowWaterModel`, `PrimitiveDryModel` and `PrimitiveWetModel`, hence a non-default output writer can be passed on as a keyword argument to the model constructor
 
-```julia
-julia> using SpeedyWeather
-julia> spectral_grid = SpectralGrid()
-julia> my_output_writer = OutputWriter(spectral_grid, PrimitiveDry)
-julia> model = PrimitiveDryModel(;spectral_grid, output=my_output_writer)
+```@example netcdf
+using SpeedyWeather
+spectral_grid = SpectralGrid()
+output = OutputWriter(spectral_grid, ShallowWater)
+model = ShallowWaterModel(;spectral_grid, output=output)
+nothing # hide
 ```
 
 So after we have defined the grid through the `SpectralGrid` object we can use and change
-the implemented `OutputWriter` by passing on the following arguments
-```julia
-julia> my_output_writer = OutputWriter(spectral_grid, PrimitiveDry, kwargs...)
-```
-the `spectral_grid` has to be the first argument then the model type
-(`Barotropic`, `ShallowWater`, `PrimitiveDry`, `PrimitiveWet`)
-which helps the output writer to make default choices on which variables to output. However, we can
-also pass on further keyword arguments. So let's start with an example.
+the implemented `OutputWriter` by passing on additional arguments.
+The `spectral_grid` has to be the first argument then the model type
+(`Barotropic`, `ShallowWater`, `PrimitiveDry`, or `PrimitiveWet`)
+which helps the output writer to make default choices on which variables to output.
+Then we can also pass on further keyword arguments. So let's start with an example.
 
 ## Example 1: NetCDF output every hour
 
-If we want to increase the frequency of the output we can choose `output_dt` (default `=6` in hours) like so
-```julia
-julia> my_output_writer = OutputWriter(spectral_grid, PrimitiveDry, output_dt=1)
-julia> model = PrimitiveDryModel(;spectral_grid, output=my_output_writer)
+If we want to increase the frequency of the output we can choose `output_dt` (default `=Hour(6)`) like so
+```@example netcdf
+output = OutputWriter(spectral_grid, ShallowWater, output_dt=Hour(1))
+model = ShallowWaterModel(;spectral_grid, output=output)
+nothing # hide
 ```
-which will now output every hour. It is important to pass on the new output writer `my_output_writer` to the
+which will now output every hour. It is important to pass on the new output writer `output` to the
 model constructor, otherwise it will not be part of your model and the default is used instead.
-Note that `output_dt` has to be understood as the minimum frequency or maximum output time step.
-Example, we run the model at a resolution of T85 and the time step is going to be 670s
-```julia
-julia> spectral_grid = SpectralGrid(trunc=85)
-julia> time_stepper = Leapfrog(spectral_grid)
-Leapfrog{Float32}:
-...
- Δt_sec::Int64 = 670
-...
+Note that the choice of `output_dt` can affect the actual time step that is used for the model
+integration, which is explained in the following.
+Example, we run the model at a resolution of T42 and the time step is going to be
+```@example netcdf
+spectral_grid = SpectralGrid(trunc=42,nlev=1)
+time_stepping = Leapfrog(spectral_grid)
+time_stepping.Δt_sec
+```
+seconds. Depending on the output frequency (we chose `output_dt = Hour(1)` above)
+this will be slightly adjusted during model initialization:
+```@example netcdf
+output = OutputWriter(spectral_grid, ShallowWater, output_dt=Hour(1))
+model = ShallowWaterModel(;spectral_grid, time_stepping, output)
+simulation = initialize!(model)
+model.time_stepping.Δt_sec
+```
+The shorter the output time step the more the model time step needs to be adjusted
+to match the desired output time step exactly. This is important so that for daily output at
+noon this does not slowly shift towards night over years of model integration.
+One can always disable this adjustment with
+```@example netcdf
+time_stepping = Leapfrog(spectral_grid,adjust_with_output=false)
+time_stepping.Δt_sec
+```
+and a little info will be printed to explain that even though you wanted
+`output_dt = Hour(1)` you will not actually get this upon initialization:
+```@example netcdf
+model = ShallowWaterModel(;spectral_grid, time_stepping, output)
+simulation = initialize!(model)
 ```
-This means that after 32 time steps 5h 57min and 20s will have passed where output will happen as
-the next time step would be >6h. The time axis of the NetCDF output will look like
-```julia
-julia> using NCDatasets
-julia> ds = NCDataset("run_0001/output.nc");
-julia> ds["time"][:]
-5-element Vector{Dates.DateTime}:
- 2000-01-01T00:00:00
- 2000-01-01T05:57:20
- 2000-01-01T11:54:40
- 2000-01-01T17:52:00
- 2000-01-01T23:49:20
-
-julia> diff(ds["time"][:])
-4-element Vector{Dates.Millisecond}:
- 21440000 milliseconds
- 21440000 milliseconds
- 21440000 milliseconds
- 21440000 milliseconds
-```
-This is so that we don't interpolate in time during output to hit exactly every 6 hours, but at
-the same time have a constant spacing in time between output time steps.
+
+The time axis of the NetCDF output will now look like
+```@example netcdf
+using NCDatasets
+model.feedback.verbose = false # hide
+run!(simulation,n_days=1,output=true)
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+ds["time"][:]
+```
+which is a bit ugly, that's why `adjust_with_output=true` is the default. In that case we would have
+```@example netcdf
+time_stepping = Leapfrog(spectral_grid,adjust_with_output=true)
+output = OutputWriter(spectral_grid, ShallowWater, output_dt=Hour(1))
+model = ShallowWaterModel(;spectral_grid, time_stepping, output)
+simulation = initialize!(model)
+run!(simulation,n_days=1,output=true)
+id = model.output.id
+ds = NCDataset("run_$id/output.nc")
+ds["time"][:]
+```
+very neatly hourly output in the NetCDF file!
 
 ## Example 2: Output onto a higher/lower resolution grid
 
@@ -75,7 +94,7 @@ So for example `output_Grid=FullClenshawGrid` and `nlat_half=48` will always int
 regular 192x95 longitude-latitude grid of 1.875˚ resolution, regardless the grid and resolution used
 for the model integration.
 ```julia
-julia> my_output_writer = OutputWriter(spectral_grid, PrimitiveDry, output_Grid=FullClenshawGrid, nlat_half=48)
+julia> my_output_writer = OutputWriter(spectral_grid, ShallowWater, output_Grid=FullClenshawGrid, nlat_half=48)
 ```
 Note that by default the output is on the corresponding full of the grid used in the dynamical core
 so that interpolation only happens at most in the zonal direction as they share the location of the
@@ -163,12 +182,46 @@ search: OutputWriter
 
     •  pkg_version::VersionNumber
 
-    •  startdate::Dates.DateTime
-
-    •  output_dt::Float64: [OPTION] output frequency, time step [hrs]
+    •  startdate::DateTime
 
-    •  output_dt_sec::Int64: actual output time step [sec]
+    •  output_dt::Second: [OPTION] output frequency, time step
 
     •  output_vars::Vector{Symbol}: [OPTION] which variables to output,
        u, v, vor, div, pres, temp, humid
+
+    •  missing_value::Union{Float32, Float64}: [OPTION] missing value to
+       be used in netcdf output
+
+    •  compression_level::Int64: [OPTION] lossless compression level;
+       1=low but fast, 9=high but slow
+
+    •  shuffle::Bool: [OPTION] shuffle/bittranspose filter for
+       compression
+
+    •  keepbits::SpeedyWeather.Keepbits: [OPTION] mantissa bits to keep
+       for every variable
+
+    •  output_every_n_steps::Int64
+
+    •  timestep_counter::Int64
+
+    •  output_counter::Int64
+
+    •  netcdf_file::Union{Nothing, NCDatasets.NCDataset}
+
+    •  input_Grid::Type{<:AbstractGrid}
+
+    •  as_matrix::Bool: [OPTION] sort grid points into a matrix
+       (interpolation-free), for OctahedralClenshawGrid, OctaHEALPixGrid
+       only
+
+    •  quadrant_rotation::NTuple{4, Int64}
+
+    •  matrix_quadrant::NTuple{4, Tuple{Int64, Int64}}
+
+    •  output_Grid::Type{<:AbstractFullGrid}: [OPTION] the grid used for
+       output, full grids only
+
+    •  nlat_half::Int64: [OPTION] the resolution of the output grid,
+       default: same nlat_half as in the dynamical core
 ```

docs/src/ringgrids.md

@@ -72,8 +72,11 @@ based on [UnicodePlots](https://github.com/JuliaPlots/UnicodePlots.jl)' `heatmap
 ```@example ringgrids
 nlat_half = 24
 grid = randn(OctahedralGaussianGrid,nlat_half)
-plot(grid)
+RingGrids.plot(grid)
 ```
+(Note that to skip the `RingGrids.` in the last line you can do `import SpeedyWeather.RingGrids: plot`,
+`import SpeedyWeather: plot` or simply `using SpeedyWeather`. It's just that `LowerTriangularMatrices`
+also defines `plot` which otherwise causes naming conflicts.)
 
 ## Indexing RingGrids
 

docs/src/setups.md

@@ -106,11 +106,16 @@ run!(simulation,n_days=6,output=true)
 ```
 The progress bar tells us that the simulation run got the identification "0001"
 (which just counts up, so yours might be higher), meaning that
-data is stored in the folder `/run_0001`, so let's plot that data properly (and not just using UnicodePlots).
+data is stored in the folder `/run_0001`. In general we can check this also via
+```@example galewsky_setup
+id = model.output.id
+```
+So let's plot that data properly (and not just using UnicodePlots). `$id` in the following just means
+that the string is interpolated to `run_0001` if this is the first unnamed run in your folder.
 ```@example galewsky_setup
 using PythonPlot, NCDatasets
 ioff() # hide
-ds = NCDataset("run_0001/output.nc")
+ds = NCDataset("run_$id/output.nc")
 ds["vor"]
 ```
 Vorticity `vor` is stored as a lon x lat x vert x time array, we may want to look at the first time step,
@@ -127,7 +132,7 @@ ax.set_xlabel("longitude")
 ax.set_ylabel("latitude")
 ax.set_title("Relative vorticity")
 tight_layout() # hide
-savefig("galewsky1.png") # hide
+savefig("galewsky1.png", dpi=70) # hide
 nothing # hide
 ```
 ![Galewsky jet pyplot1](galewsky1.png)
@@ -139,7 +144,7 @@ And now the last time step, that means time = 12days is
 t = ds.dim["time"]
 vor = Matrix{Float32}(ds["vor"][:,:,1,t])
 ax.pcolormesh(lon,lat,vor')
-savefig("galewsky2.png") # hide
+savefig("galewsky2.png", dpi=70) # hide
 nothing # hide
 ```
 ![Galewsky jet pyplot2](galewsky2.png)
@@ -162,7 +167,7 @@ simulation = initialize!(model)
 run!(simulation,n_days=12,output=true)
 ```
 
-This time the run got a new run id, which you see in the progress bar, but can also always check
+This time the run got a new run id, which you see in the progress bar, but can again always check
 after the `run!` call (the automatic run id is only determined just before the main time loop starts)
 with `model.output.id`, but otherwise we do as before.
 ```@example galewsky_setup
@@ -180,7 +185,7 @@ ax.set_xlabel("longitude")
 ax.set_ylabel("latitude")
 ax.set_title("Relative vorticity")
 tight_layout() # hide
-savefig("galewsky3.png") # hide
+savefig("galewsky3.png", dpi=70) # hide
 nothing # hide
 ```
 ![Galewsky jet pyplot3](galewsky3.png)
@@ -198,7 +203,7 @@ using SpeedyWeather
 spectral_grid = SpectralGrid(trunc=63,nlev=1)
 forcing = JetStreamForcing(spectral_grid,latitude=60)
 drag = QuadraticDrag(spectral_grid)
-output = OutputWriter(spectral_grid,ShallowWater,output_dt=6,output_vars=[:u,:v,:pres,:orography])
+output = OutputWriter(spectral_grid,ShallowWater,output_dt=Hour(6),output_vars=[:u,:v,:pres,:orography])
 model = ShallowWaterModel(;spectral_grid,output,drag,forcing)
 simulation = initialize!(model)
 model.feedback.verbose = false # hide
@@ -233,7 +238,7 @@ ax.set_ylabel("latitude")
 ax.set_title("Zonal wind [m/s]")
 colorbar(q,ax=ax)
 tight_layout() # hide
-savefig("polar_jets.png") # hide
+savefig("polar_jets.png", dpi=70) # hide
 nothing # hide
 ```
 ![Polar jets pyplot](polar_jets.png)
@@ -247,11 +252,11 @@ using Random # hide
 Random.seed!(1234) # hide
 using SpeedyWeather
 spectral_grid = SpectralGrid(trunc=127,nlev=1)
-time_stepping = SpeedyWeather.Leapfrog(spectral_grid,Δt_at_T31=30)
+time_stepping = SpeedyWeather.Leapfrog(spectral_grid,Δt_at_T31=Minute(30))
 implicit = SpeedyWeather.ImplicitShallowWater(spectral_grid,α=0.5)
 orography = EarthOrography(spectral_grid,smoothing=false)
 initial_conditions = SpeedyWeather.RandomWaves()
-output = OutputWriter(spectral_grid,ShallowWater,output_dt=12,output_vars=[:u,:pres,:div,:orography])
+output = OutputWriter(spectral_grid,ShallowWater,output_dt=Hour(12),output_vars=[:u,:pres,:div,:orography])
 model = ShallowWaterModel(;spectral_grid,orography,output,initial_conditions,implicit,time_stepping)
 simulation = initialize!(model)
 model.feedback.verbose = false # hide
@@ -302,7 +307,7 @@ ax.set_xlabel("longitude")
 ax.set_ylabel("latitude")
 ax.set_title("Divergence")
 tight_layout() # hide
-savefig("gravity_waves.png") # hide
+savefig("gravity_waves.png", dpi=70) # hide
 nothing # hide
 ```
 ![Gravity waves pyplot](gravity_waves.png)
@@ -349,7 +354,7 @@ ax.set_xlabel("longitude")
 ax.set_ylabel("latitude")
 ax.set_title("Surface relative vorticity")
 tight_layout() # hide
-savefig("jablonowski.png") # hide
+savefig("jablonowski.png", dpi=70) # hide
 nothing # hide
 ```
 ![Jablonowski pyplot](jablonowski.png)

docs/src/speedytransforms.md

@@ -53,6 +53,7 @@ By default, the `gridded` transforms onto a [`FullGaussianGrid`](@ref FullGaussi
 into a vector west to east, starting at the prime meridian, then north to south, see [RingGrids](@ref).
 We can visualize `map` quickly with a UnicodePlot via `plot` (see [Visualising RingGrid data](@ref))
 ```@example speedytransforms
+import SpeedyWeather.RingGrids: plot    # not necessary when `using SpeedyWeather`
 plot(map)
 ```
 Yay! This is the what the ``l=m=1`` spherical harmonic is supposed to look like!

src/SpeedyWeather.jl

@@ -4,6 +4,7 @@ module SpeedyWeather
 using DocStringExtensions
 
 # NUMERICS
+import Primes
 import Random
 import FastGaussQuadrature
 import LinearAlgebra: LinearAlgebra, Diagonal
@@ -17,7 +18,7 @@ import Adapt: Adapt, adapt, adapt_structure
 
 # INPUT OUTPUT
 import TOML
-import Dates: Dates, DateTime
+import Dates: Dates, DateTime, Millisecond, Second, Minute, Hour, Day
 import Printf: Printf, @sprintf
 import NCDatasets: NCDatasets, NCDataset, defDim, defVar
 import JLD2: jldopen
@@ -27,7 +28,7 @@ import UnicodePlots
 import ProgressMeter
 
 # to avoid a `using Dates` to pass on DateTime arguments
-export DateTime
+export DateTime, Second, Minute, Hour, Day
 
 # EXPORT MONOLITHIC INTERFACE TO SPEEDY
 export  run_speedy,