[JOSS review] Nat's Review of SpeedyWeather.jl

[natgeo-wong]

Review of SpeedyWeather.jl by Nat (openjournals/joss-reviews#6323 (comment))

See the pdf for detailed comments.
reviewv1.pdf

Current outstanding issues on the checklist:

• Functionality: Have the functional claims of the software been confirmed?
  • We need a more detailed explanation of the different modules / AbstractTypes available for the different models
  • which of theseAbstractType modules are relevant for each of the different types of models (e.g. ShallowWaterModel, PrimitiveDryModel) available in SpeedyWeather.jl
• Functionality documentation: Is the core functionality of the software documented to a satisfactory level (e.g., API method documentation)?
  • A reorganization of the documentation into a more hierarchical form
• Community guidelines: Are there clear guidelines for third parties wishing to 1) Contribute to the software 2) Report issues or problems with the software 3) Seek support

[natgeo-wong]

openjournals/joss-reviews#6323

[milankl]

Nat also suggests looking into Paulius and Garner, 2006 for a simple radiation scheme

[milankl]

Nat also suggests looking into Paulius and Garner, 2006 for a simple radiation scheme

implemented with #444

[maximilian-gelbrecht]

@milankl Let me know if you want me to do something

[milankl]

I can create a type tree for the documentation and outline more clearly what new types within that tree need to fulfill in order to be used properly, e.g. what would a MyConvection <: AbstractConvection need to have/do in order to serve as an alternative to our current SimplifiedBettsMiller convection?

However, @maximilian-gelbrecht what do you think our community guidelines are or should be? I think "open an issue" is the answer to all of them, but maybe we can just make that clearer in the readme or documentation?

[maximilian-gelbrecht]

I do think this is also a problem many Julia library have, so much so that I think that Julia should have some kind of built-in system that handles that. But that's probably besides the point here. Personally, I'd put it in the doc string of each abstract type.

However, @maximilian-gelbrecht what do you think our community guidelines are or should be? I think "open an issue" is the answer to all of them, but maybe we can just make that clearer in the readme or documentation?

We could think about enabling the GitHub Discussions as well, but currently I think it's fine to just use GitHub Issues for all of it. But yes, we should add a paragraph prominently to the documentation that encourages people to just open an issue if there are any questions, bugs, ideas, etc. We should also add an CONTRIBUTING.md to the repo

[navidcy]

But yes, we should add a paragraph prominently to the documentation that encourages people to just open an issue if there are any questions, bugs, ideas, etc.

Agreed

We should also add an CONTRIBUTING.md to the repo

Definitely

[navidcy]

Regarding @natgeo-wong comment:

For some reason, for the Shallow Water w/ Mountains example, the model crashed when I changed the resolution from T63 to T85. Is there a reason? Do I need to specify a smaller timestep for stability? If so, perhaps it would be useful to put this somewhere, maybe we should have a section/subsection on changing the timestep to ensure stability.

@natgeo-wong I understand you are talking about https://speedyweather.github.io/SpeedyWeather.jl/dev/setups/#Shallow-water-with-mountains

Personally I can change the resolution without any issue. Look at results below using #main.

using SpeedyWeather

for resolution in [63, 85, 127]
    @info "resolution is T$(resolution)"
    @show spectral_grid = SpectralGrid(trunc=resolution,nlev=1)
    orography = NoOrography(spectral_grid);
    initial_conditions = ZonalJet();
    @show model.time_stepping
    model = ShallowWaterModel(;spectral_grid, orography, initial_conditions);
    simulation = initialize!(model);
    run!(simulation,period=Day(6));
end

[ Info: resolution is T63
spectral_grid = SpectralGrid(trunc = resolution, nlev = 1) = SpectralGrid:
├ Spectral:   T63 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid:       96-ring OctahedralGaussianGrid{Float32}, 10944 grid points
├ Resolution: 216km (average), 192km × 207km (Equator)
└ Vertical:   1-level SigmaCoordinates
model.time_stepping = Leapfrog{Float32} <: TimeStepper
├ trunc::Int64 = 63
├ Δ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 = 900000 milliseconds
├ Δt_sec::Float32 = 900.0
└ Δt::Float32 = 0.0001412651
Weather is speedy: 100%|██████████████████████████████████████████████████████████| Time: 0:00:00 ( 3.37 millenia/day)
[ Info: resolution is T85
spectral_grid = SpectralGrid(trunc = resolution, nlev = 1) = SpectralGrid:
├ Spectral:   T85 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid:       128-ring OctahedralGaussianGrid{Float32}, 18688 grid points
├ Resolution: 165km (average), 147km × 156km (Equator)
└ Vertical:   1-level SigmaCoordinates
model.time_stepping = Leapfrog{Float32} <: TimeStepper
├ trunc::Int64 = 85
├ Δ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 = 675000 milliseconds
├ Δt_sec::Float32 = 675.0
└ Δt::Float32 = 0.00010594883
Weather is speedy: 100%|██████████████████████████████████████████████████████████| Time: 0:00:01 (1279.78 years/day)
[ Info: resolution is T127
spectral_grid = SpectralGrid(trunc = resolution, nlev = 1) = SpectralGrid:
├ Spectral:   T127 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid:       192-ring OctahedralGaussianGrid{Float32}, 40320 grid points
├ Resolution: 112km (average), 100km × 104km (Equator)
└ Vertical:   1-level SigmaCoordinates
model.time_stepping = Leapfrog{Float32} <: TimeStepper
├ trunc::Int64 = 127
├ Δ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 = 450000 milliseconds
├ Δt_sec::Float32 = 450.0
└ Δt::Float32 = 7.063255e-5
Weather is speedy: 100%|██████████████████████████████████████████████████████████| Time: 0:00:04 (303.66 years/day)

Notice that when we construct the model, the tilmestep is different for different resolutions. I believe the CFL-type of instability you are hinting is taken care of by this automatic timestep change. Isn't it @milankl?

My machine details

julia> versioninfo()
Julia Version 1.10.1
Commit 7790d6f064 (2024-02-13 20:41 UTC)
Build Info:

    Note: This is an unofficial build, please report bugs to the project
    responsible for this build and not to the Julia project unless you can
    reproduce the issue using official builds available at https://julialang.org/downloads

Platform Info:
  OS: macOS (arm64-apple-darwin23.3.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 8 default, 0 interactive, 4 GC (on 8 virtual cores)
Environment:
  JULIA_EDITOR = code

julia> using Pkg; Pkg.status()
Project SpeedyWeather v0.8.0
Status `~/Research/SpeedyWeather.jl/Project.toml`
  [621f4979] AbstractFFTs v1.5.0
  [79e6a3ab] Adapt v4.0.1
  [2119f1ac] AssociatedLegendrePolynomials v1.0.1
  [de688a37] BitInformation v0.6.1
  [052768ef] CUDA v5.2.0
⌃ [944b1d66] CodecZlib v0.7.3
  [ffbed154] DocStringExtensions v0.9.3
⌃ [7a1cc6ca] FFTW v1.7.2
  [cc61a311] FLoops v0.2.1
  [442a2c76] FastGaussQuadrature v1.0.2
  [a8297547] GenericFFT v0.1.6
⌃ [033835bb] JLD2 v0.4.44
  [63c18a36] KernelAbstractions v0.9.16
⌃ [85f8d34a] NCDatasets v0.14.0
  [27ebfcd6] Primes v0.5.5
  [92933f4c] ProgressMeter v1.9.0
  [b8865327] UnicodePlots v3.6.3
  [ade2ca70] Dates
  [37e2e46d] LinearAlgebra
  [de0858da] Printf
  [9a3f8284] Random
  [10745b16] Statistics v1.10.0
  [fa267f1f] TOML v1.0.3
Info Packages marked with ⌃ have new versions available and may be upgradable.

[milankl]

Yes, the timestep is

1. automatically chosen given the resolution. At T31 it currently uses 30min which is the default set in model.time_stepping.Δt_at_T31. It then scales ~linearly to your resolution
2. Given the output frequency, the timestep might be slightly adjusted to have e.g. output exactly every hour without interpolation. This rarely influences the time step by more than +- a few percent. model.time_stepping.adjust_with_output = false can disable that, also a warning is thrown if that's more than +5% or less than -5%

Overall the time step is part of the numerics that SpeedyWeather should hide from the general user. A general user shouldn't worry about time stepping, scaling, transforms, etc. Sometimes that's inevitable, but that should be our mantra!

[milankl]

@natgeo-wong As you have may seen in the recent pull requests the documentation is growing to address your points (feel free to check it out and comment back, appreciate any feedback!), currently missing are still your points from above, which I plan to adress next week, then you'll get a proper point-by-point response. Many thanks!!!

[kthyng]

@milankl How is this work coming along?

[milankl]

@kthyng Thanks for the reminder, it's just a lot of documentation to write... I'm still working on @natgeo-wong's points (any feedback welcome to our recent progress on this issue!) who rightly requested to have more documentation, in particular we've added now

• Vertical diffusion Vertical diffusion documentation #532
• Surface fluxes Surface flux documentation #530
• Orography Orography documented #526
• Land-sea mask and aquaplanet simulations ConstantOceanClimatology, AquaPlanet #496
• Held-Suarez simulations Held-Suarez example, NoTypes for surface fluxes #514
• Radiation Radiation documentation #513
• Convection Convection documentation #512
• Large-scale condensation Large-scale condensation documented #510
• Using gradient operators from SpeedyWeather Exports more gradient operators from SpeedyTransforms #523
• Online model analysis with operators in SpeedyWeather and callbacks Add to analysis doc (from new branch) #505, Document typical online analysis steps #497
• Model and Simulation structure Field trees #492
• Ocean model Ocean model documentation #533
• An overview of how to define new parameterizations Parameterization structure documented #534

[milankl]

Responses to individual points @natgeo-wong

1. In light of the above, the page Model Setups in the documentation should be re-
   named as Examples. The page Model Setups should more follow the style of the Oceananigans.jl documentation, in that it specifies and expands upon explaining the
   different AbstractTypes that are available for the different model types.

We now have 2D and 3D examples

• https://speedyweather.github.io/SpeedyWeather.jl/dev/examples_2D/
• https://speedyweather.github.io/SpeedyWeather.jl/dev/examples_3D/

introducing some of the types/components users can change and play around with. Others are discussed in other parts of the documentation, e.g. https://speedyweather.github.io/SpeedyWeather.jl/dev/land_sea_mask/#Time-dependent-land-sea-mask documents a time-dependent land-sea mask, or https://speedyweather.github.io/SpeedyWeather.jl/dev/analysis/ explains how to check for conservation of mass, energy etc on the fly using gradient operators and the spectral transform for spatial integrations: https://speedyweather.github.io/SpeedyWeather.jl/dev/gradients/

2. There is a lot of explanation on the theoretical physics that underlies this model in the documentation (i.e., a page dedicated to the physics behind say, the Primitive Equation model). I recommend grouping these all under a single category in the documentation named Model Physics.

We did this. Thanks for the suggestion.

3. While it is always useful to keep a reference list of API of Functions and Types at the end, some of these Functions and Types would be better placed at where the documentation either explains the Model Physics or somewhere within the Model Setup.

Thanks for the suggestion but in the end we decided against this because we wanted to separate the theoretical/physical/mathematical aspects from the implementation details, API, function signatures. We now usually have some example outlining how a given function/type is supposed to be used or extended and in that sense one can always request the docstring in the REPL like ?model or search the documentation if the usage of a function is not clear. I personally prefer to have most parts of the documentation example-driven instead of hitting a user with all optional keywords to a function which the docstring, however, definitely should provide. But please let us know if there's anything in the documentation where any of this is not clear.

4. There needs to be a section somewhere in the Documentation or README on how to officially contribute to SpeedyWeather.jl, even if it’s just a note to make an issue/pull-request

We have added this to the readme and the documentation. There's many ways to contribute but we encourage people being interested to open an issue or reach out to us individually.

5. what specific AbstractType modules are planned to be made available in Speedy-
   Weather.jl for the different physical parameterizations

See https://speedyweather.github.io/SpeedyWeather.jl/dev/parameterizations/ and also more in answers below.

6. what are the standard fields required for each of these AbstractType modules

Only few have required fields, like orography, see https://speedyweather.github.io/SpeedyWeather.jl/dev/orography/, others like a custom coriolis struct would need to have a field f::Vector but we don't really encourage people to change them although it's certainly possible to do (e.g. you can mathematically think of a planet where northern and southern hemisphere rotate in opposite directions, meaning you have low pressure system rotating anti-clockwise on both hemispheres). Whenever there is a required field we specify this in the respective part of the documentation.

7. which of these AbstractType modules are relevant for each of the different types of
   models (e.g. ShallowWaterModel, PrimitiveDryModel) available in SpeedyWeather.jl

For the parameterizations we document the abstract types now in https://speedyweather.github.io/SpeedyWeather.jl/dev/parameterizations/

Technically you can always do

julia> subtypes(SpeedyWeather.AbstractModelComponent)
16-element Vector{Any}:
 SpeedyWeather.AbstractAdiabaticConversion
 SpeedyWeather.AbstractAlbedo
 SpeedyWeather.AbstractAtmosphere
 SpeedyWeather.AbstractCoriolis
 SpeedyWeather.AbstractDrag
 SpeedyWeather.AbstractForcing
 SpeedyWeather.AbstractGeopotential
 SpeedyWeather.AbstractHorizontalDiffusion
 SpeedyWeather.AbstractImplicit
 SpeedyWeather.AbstractInitialConditions
 SpeedyWeather.AbstractOrography
 SpeedyWeather.AbstractParameterization
 SpeedyWeather.AbstractParticleAdvection
 SpeedyWeather.AbstractPlanet
 SpeedyWeather.AbstractSchedule
 SpeedyWeather.AbstractTimeStepper

julia> subtypes(SpeedyWeather.AbstractParameterization)
12-element Vector{Any}:
 SpeedyWeather.AbstractBoundaryLayer
 SpeedyWeather.AbstractClausiusClapeyron
 SpeedyWeather.AbstractCondensation
 SpeedyWeather.AbstractConvection
 SpeedyWeather.AbstractRadiation
 SpeedyWeather.AbstractSurfaceEvaporation
 SpeedyWeather.AbstractSurfaceHeatFlux
 SpeedyWeather.AbstractSurfaceThermodynamics
 SpeedyWeather.AbstractSurfaceWind
 SpeedyWeather.AbstractTemperatureRelaxation
 SpeedyWeather.AbstractVegetation
 SpeedyWeather.AbstractVerticalDiffusion

But we are hesitating to put this list into the documentation as it sounds like we encourage users to just go ahead and define new concrete types for any of them. While some are straightforward to subtype, others are more complicated and so we'd rather explain those that we encourage users to experiment with and document them well but slowly growing the documentation over time.

8. My preference would be for the creation of detailed documentation for at least two or
   three different modules before this paper is accepted for publication. The documentation
   for the AbstractPlanet and AbstractOrography superTypes, for example, could be pretty
   helpful.

Documented now are:

• Forcing and drag, https://speedyweather.github.io/SpeedyWeather.jl/dev/forcing_drag/
• Parameterizations https://speedyweather.github.io/SpeedyWeather.jl/dev/parameterizations/
• Orography: https://speedyweather.github.io/SpeedyWeather.jl/dev/orography/
• Land-sea mask: https://speedyweather.github.io/SpeedyWeather.jl/dev/land_sea_mask/
• Ocean: https://speedyweather.github.io/SpeedyWeather.jl/dev/ocean/
• Callbacks: https://speedyweather.github.io/SpeedyWeather.jl/dev/callbacks/

On AbstractPlanet: We don't document or expect someone to define another planet like Jupiter <: AbstractPlanet, while this is mathematically straight forward, we haven't tested any of this (see #462) and so don't want to document something that we don't know works. When we can guarantee that a planet with another radius, rotation, gravity does not cause numerical stabilities (it likely will) then we'll document this. Otherwise we're afraid that too many users will run into instabilities. In other words, while easy to implement, other planets are currently not supported from our side and entirely experimental.

9. And also, a question I had is, what is the criteria that determines whether something
   needs to become modularized? For example, I believe that SpeedyWeather does have the
   AbstractOrography superType, but within this superType after some exploring I have seen
   that you can point towards different files for orography. And if you change the file for orog-
   raphy file, technically, you can run drastically non-Earthlike orographies in your experiment.
   So why do you need an AbstractOrography superType?

From a software engineering perspective it's easier and clearer to implement a very similar structure for (almost) all components. That's why we have abstract VerticalCoordinates, AbstractTimeStepper, AbstractImplicit although implementing a new concrete type is anything but straight forward. Some components are easy to abstract and make modular, others like the timestepping are so deeply ingrained into the model structure that it's more work to abstract them (and hence provide an API for custom) than we currently believe is worth it. So these "criteria" you are requesting are somewhat a grey zone: Some components feel very natural to abstract and modularise others require much more effort or are components that we believe a user would not want to experiment with. But so depending on future demands this line can shift: At the moment we don't document or expect someone to define another planet like Jupiter <: AbstractPlanet, while this is mathematically straight forward, we haven't tested any of this (see #462) and so don't want to document something that we don't know works.

On Orography: We have EarthOrography, ZonalRidge and NoOrography. You are right than one can just use another file within EarthOrography and get, e.g. Mars' orography. However, maybe someone wants to define some orography as a function of lat, lon or maybe also orography with some randomness (rough ground etc). Orography is now documented here https://speedyweather.github.io/SpeedyWeather.jl/dev/orography/

10. In the Primitive Equations, there is a note that radiation is subsumed into a forcing
    term Q/cp. Is there any plan to add a radiative scheme, and which scheme? And right
    now, for clarification purposes, the model is entirely driven by the surface temperatures?

We have now several radiation schemes, see #502 #444 even though they are still very simple. More complicated ones on the todo list are #461, the Frierson 2006 2-band radiation (grey) and Fortran SPEEDY's 4-band radiation scheme (humidity and clouds impact radiation, +CO2 forcing). A simple TransparentShortwave radiation can force the model, but if disabled the forcing comes from sea surface temperatures indeed. The current ones are also documented: https://speedyweather.github.io/SpeedyWeather.jl/dev/radiation/. While there is much more todo on this, we agree, we believe it's beyond the scope of a JOSS paper and more something towards a model evaluation paper. So we're in favour of not implementing more complicated radiation schemes at this stage.

[natgeo-wong]

Very nice work and detailed documentation! I have one passing comment, is that certain parameterizations and physics may only be relevant to PrimitiveModel Types, and not ShallowWater or Barotropic types. Maybe some demarcation of this would help? Other than that, all things look good to me.

[milankl]

All parameterizations act on vertical columns and are therefore only available for PrimitiveEquation as the others are 2D models where any additional terms have to be formulated as forcing/drag. Hmmm maybe that's not clear?

[natgeo-wong]

Maybe it should be mentioned somewhere?

[milankl]

@natgeo-wong just added a prominent info block in the Parameterization docs for this! ☝🏼 Hope this is what you imagined and we can close this issue?

[milankl]

Sorry was automatically closed with the pr merged. You should close it

[natgeo-wong]

Okay, I had a look. Perhaps you should split the "Forcing and Drag" part into 2: (1) Introduction to Extensions and (2) Forcing and Drag.

Also, I suggest a reordering. "Running SpeedyWeather" and "Extending SpeedyWeather" should probably come before the documentation on the physics and dynamics. But I will go ahead and close this issue.

[milankl]

Good suggestion, I'll address those and references the pull requests here!