diff --git a/.gitignore b/.gitignore index 3cc979bbb..533680177 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,9 @@ deps/src/ docs/build/ docs/site/ +# PythonCall added .CondaPkg to /docs +docs/.CondaPkg + # File generated by Pkg, the package manager, based on a corresponding Project.toml # It records a fixed state of all packages used by the project. As such, it should not be # committed for packages, but should be committed for applications that require a static @@ -32,4 +35,4 @@ Manifest.toml run_*/ # no video outputs -*.mp4 +*.mp4 \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml index eb86049e5..c6c3b8cef 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,11 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +PythonPlot = "274fc56d-3b97-40fa-a1cd-1b4a50311bf9" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] Documenter = "0.26, 0.27" -UnicodePlots = "^3.3" \ No newline at end of file +NCDatasets = "0.12" +UnicodePlots = "^3.3" +PythonPlot = "1" diff --git a/docs/make.jl b/docs/make.jl index 427a42577..68e47a329 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,8 +1,8 @@ using Documenter, SpeedyWeather makedocs( - format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing)=="true", - ansicolor=true), + format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing)=="true", + ansicolor=true), sitename="SpeedyWeather.jl", authors="M Klöwer and SpeedyWeather contributors", modules=[SpeedyWeather], @@ -20,8 +20,6 @@ makedocs( "Submodule: RingGrids"=>"ringgrids.md", "Submodule: LowerTriangularMatrices"=>"lowertriangularmatrices.md", "Submodule: SpeedyTransforms"=>"speedytransforms.md", - "Style and convention guide"=>"conventions.md", - "Development notes"=>"development.md", "Function and type index"=>"functions.md"] ) diff --git a/docs/src/conventions.md b/docs/src/conventions.md deleted file mode 100644 index 37a80f721..000000000 --- a/docs/src/conventions.md +++ /dev/null @@ -1,60 +0,0 @@ -# Style and convention guide - -In SpeedyWeather.jl we've been following the several conventions that are documented here. - -## Variable naming - -The prognostic variables in spectral space are called - -```julia - vor # Vorticity of horizontal wind field - div # Divergence of horizontal wind field - temp # Absolute temperature [K] - pres_surf # Logarithm of surface pressure [log(Pa)] - humid # Specific humidity [g/kg] -``` - -their transforms into grid-point space get a `_grid` suffix, their tendencies a `_tend` suffix. Further derived diagnostic dynamic variables are - -```julia - u - v - geopot - ... -``` - -## Preallocation - -All arrays representing variables are preallocated and grouped into two structs - -```julia - Prog::PrognosticVariables - Diag::DiagnosticVariables -``` - -The `Diag` struct contains further structs which represent the grid-point transformations of the prognostic variables and their tendencies. - -```julia - gridvars::GridVariables - tendencies::Tendencies - ... -``` - -Constant arrays are grouped into several structs - -```julia -Boundaries -``` - -## Julian conventions - -We follow Julia's [style guide](https://docs.julialang.org/en/v1/manual/style-guide/#Style-Guide) and highlight here some important aspects of it. - -- __Bang (!) convention__. A function `func` does not change its input arguments, however, `func!` does. -Hence, `func!` is often the in-place version of `func`, avoiding as much memory allocation as possible -and often changing its first argument, e.g. `func!(out,in)` so that argument `in` is used to calculate -`out` which has been preallocated before function call. -- __Number format flexibility__. Numeric literals such as `2.0` or `1/3` are only used in the model setup -but avoided throughout the code to obtain a fully number format-flexible package using the number format -`NF` as a compile-time variable throughout the code. This often leads to overly specific code whereas -a `Real` would generally suffice. However, this is done to avoid any implicit type conversions. \ No newline at end of file diff --git a/docs/src/development.md b/docs/src/development.md deleted file mode 100644 index fc14abb17..000000000 --- a/docs/src/development.md +++ /dev/null @@ -1,30 +0,0 @@ -# Development notes - -To run tests, from the path of your local clone of the repository do: - -``` -julia --project=. -e 'import Pkg; Pkg.test()' -``` - -To install dependencies: - -``` -julia --project -e 'import Pkg; Pkg.instantiate()` -``` - -then opening: - -``` -julia --project -``` - -you are able to `using SpeedyWeather`. - -To generate the docs: - -``` -julia --project=docs -e 'import Pkg; Pkg.develop(path="."); Pkg.instantiate()' -julia --project=docs docs/make.jl -``` - -If the docs are generated successfully, you view them by opening `docs/build/index.html` in your favorite browser. diff --git a/docs/src/how_to_run_speedy.md b/docs/src/how_to_run_speedy.md index fadba00f1..c9729ee06 100644 --- a/docs/src/how_to_run_speedy.md +++ b/docs/src/how_to_run_speedy.md @@ -14,59 +14,44 @@ but let's start with some examples first. We want to use the barotropic model to simulate some free-decaying 2D turbulence on the sphere without rotation. We start by defining the `SpectralGrid` object. -To have a resolution of about 100km, we choose a spectral resolution of -T127 (see [Available horizontal resolutions](@ref)) and `nlev=1` vertical levels. +To have a resolution of about 200km, we choose a spectral resolution of +T63 (see [Available horizontal resolutions](@ref)) and `nlev=1` vertical levels. The `SpectralGrid` object will provide us with some more information -```julia -julia> spectral_grid = SpectralGrid(trunc=127,nlev=1) -SpectralGrid: - Spectral: T127 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m - Grid: 40320-element, 192-ring OctahedralGaussianGrid{Float32} (quadratic) - Resolution: 112km (average) - Vertical: 1-level SigmaCoordinates +```@example howtorun +using SpeedyWeather +spectral_grid = SpectralGrid(trunc=63,nlev=1) ``` + We could have specified further options, but let's ignore that for now. Next step we create a planet that's like Earth but not rotating -```julia -julia> still_earth = Earth(rotation=0) -Main.SpeedyWeather.Earth - rotation: Float64 0.0 - gravity: Float64 9.81 - daily_cycle: Bool true - length_of_day: Float64 24.0 - seasonal_cycle: Bool true - length_of_year: Float64 365.25 - equinox: Dates.DateTime - axial_tilt: Float64 23.4 +```@example howtorun +still_earth = Earth(rotation=0) ``` There are other options to create a planet but they are irrelevant for the barotropic vorticity equations. We also want to specify the initial conditions, randomly distributed vorticity is already defined -```julia -julia> initial_conditions = StartWithRandomVorticity() -StartWithRandomVorticity - power_law: Float64 -3.0 - amplitude: Float64 1.0e-5 +```@example howtorun +using Random # hide +Random.seed!(1234) # hide +initial_conditions = StartWithRandomVorticity() ``` By default, the power of vorticity is spectrally distributed with ``k^{-3}``, ``k`` being the -horizontal wavenumber, and the amplitude is ``10^{-5}\text{ s}^{-1}``. +horizontal wavenumber, and the amplitude is ``10^{-5}\text{s}^{-1}``. Now we want to construct a `BarotropicModel` with these -```julia -julia> model = BarotropicModel(;spectral_grid, initial_conditions, planet=still_earth); +```@example howtorun +model = BarotropicModel(;spectral_grid, initial_conditions, planet=still_earth) +nothing # hide ``` The `model` contains all the parameters, but isn't initialized yet, which we can do -with and then run it. -```julia -julia> simulation = initialize!(model); -julia> run!(simulation,n_days=30) -``` -The `run!` command will always return the prognostic variables, which, by default, are +with and then run it. The `run!` command will always return the prognostic variables, which, by default, are plotted for surface relative vorticity with a unicode plot. The resolution of the plot is not necessarily representative but it lets us have a quick look at the result - -![Barotropic vorticity unicode plot](https://raw.githubusercontent.com/SpeedyWeather/SpeedyWeather.jl/main/docs/img/barotropic_vorticity.jpg) +```@example howtorun +simulation = initialize!(model) +run!(simulation,n_days=20) +``` Woohoo! I can see turbulence! You could pick up where this simulation stopped by simply doing `run!(simulation,n_days=50)` again. We didn't store any output, which @@ -78,108 +63,106 @@ with default settings. More options on output in [NetCDF output](@ref). As a second example, let's investigate the Galewsky et al.[^1] test case for the shallow water equations with and without mountains. As the shallow water system has also only one level, we can reuse the `SpectralGrid` from Example 1. -```julia -julia> spectral_grid = SpectralGrid(trunc=127,nlev=1) +```@example howtorun +spectral_grid = SpectralGrid(trunc=63,nlev=1) ``` Now as a first simulation, we want to disable any orography, so we create a `NoOrography` -```julia -julia> orography = NoOrography(spectral_grid) -NoOrography{Float32, OctahedralGaussianGrid{Float32}} +```@example howtorun +orography = NoOrography(spectral_grid) ``` Although the orography is zero, you have to pass on `spectral_grid` so that it can still initialize zero-arrays of the right size and element type. Awesome. This time the initial conditions should be set the the Galewsky et al.[^1] zonal jet, which is already defined as -```julia -julia> initial_conditions = ZonalJet() -ZonalJet - latitude: Float64 45.0 - width: Float64 19.28571428571429 - umax: Float64 80.0 - perturb_lat: Float64 45.0 - perturb_lon: Float64 270.0 - perturb_xwidth: Float64 19.098593171027442 - perturb_ywidth: Float64 3.819718634205488 - perturb_height: Float64 120.0 +```@example howtorun +initial_conditions = ZonalJet() ``` The jet sits at 45˚N with a maximum velocity of 80m/s and a perturbation as described in their paper. Now we construct a model, but this time a `ShallowWaterModel` +```@example howtorun +model = ShallowWaterModel(;spectral_grid, orography, initial_conditions) +simulation = initialize!(model) +run!(simulation,n_days=6) ``` -julia> model = ShallowWaterModel(;spectral_grid, orography, initial_conditions); -julia> simulation = initialize!(model); -``` -![Galewsky jet unicode plot](https://raw.githubusercontent.com/SpeedyWeather/SpeedyWeather.jl/main/docs/img/galewsky.jpg) - Oh yeah. That looks like the wobbly jet in their paper. Let's run it again for another 6 days but this time also store [NetCDF output](@ref). -``` -julia> run!(simulation,n_days=6,output=true) -Weather is speedy: run 0002 100%|███████████████████████| Time: 0:00:12 (115.37 years/day) -``` -The progress bar tells us that the simulation run got the identification "0002", meaning that -data is stored in the folder `/run_0002`, so let's plot that data properly (and not just using UnicodePlots). -```julia -julia> using PyPlot, NCDatasets -julia> ds = NCDataset("run_0002/output.nc"); -julia> ds["vor"] -vor (384 × 192 × 1 × 25) - Datatype: Float32 - Dimensions: lon × lat × lev × time - Attributes: - units = 1/s - missing_value = NaN - long_name = relative vorticity - _FillValue = NaN -``` -Vorticity `vor` is stored as a 384x192x1x25 array, we may want to look at the first time step, +```@example howtorun +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). +```@example howtorun +using PythonPlot, NCDatasets +ioff() # hide +ds = NCDataset("run_0001/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, which is the end of the previous simulation (time=6days) which we didn't store output for. -```julia -julia> vor = ds["vor"][:,:,1,1]; -julia> lat = ds["lat"][:]; -julia> lon = ds["lon"][:]; -julia> pcolormesh(lon,lat,vor') -``` -Which looks like - -![Galewsky jet pyplot](https://raw.githubusercontent.com/SpeedyWeather/SpeedyWeather.jl/main/docs/img/galewsky_nc_6days.png) - -You see that the unicode plot heavily coarse-grains the simulation, well it's unicode after all! -And now the last time step, that means time=12days is -```julia -julia> vor = ds["vor"][:,:,1,25]; -julia> pcolormesh(lon,lat,vor') -``` - -![Galewsky jet pyplot](https://raw.githubusercontent.com/SpeedyWeather/SpeedyWeather.jl/main/docs/img/galewsky_nc_12days.png) +```@example howtorun +time = 1 +vor = Matrix{Float32}(ds["vor"][:,:,1,time]) # convert from Matrix{Union{Missing,Float32}} to Matrix{Float32} +lat = ds["lat"][:] +lon = ds["lon"][:] +pcolormesh(lon,lat,vor') +xlabel("longitude") +ylabel("latitude") +tight_layout() # hide +savefig("galewsky1.png") # hide +nothing # hide +``` +![Galewsky jet pyplot](galewsky1.png) + +You see that in comparison the unicode plot heavily coarse-grains the simulation, well it's unicode after all! +And now the last time step, that means time = 12days is + +```@example howtorun +time = 25 +vor = Matrix{Float32}(ds["vor"][:,:,1,time]) +pcolormesh(lon,lat,vor') +xlabel("longitude") +ylabel("latitude") +tight_layout() # hide +savefig("galewsky2.png") # hide +nothing # hide +``` +![Galewsky jet pyplot](galewsky2.png) The jet broke up into many small eddies, but the turbulence is still confined to the northern hemisphere, cool! How this may change when we add mountains (we had `NoOrography` above!), say Earth's orography, you may ask? Let's try it out! We create an `EarthOrography` struct like so -```julia -julia> orography = EarthOrography(spectral_grid) -EarthOrography{Float32, OctahedralGaussianGrid{Float32}}: - path::String = SpeedyWeather.jl/input_data - file::String = orography_F512.nc - scale::Float64 = 1.0 - smoothing::Bool = true - smoothing_power::Float64 = 1.0 - smoothing_strength::Float64 = 0.1 - smoothing_truncation::Int64 = 85 +```@example howtorun +orography = EarthOrography(spectral_grid) ``` It will read the orography from file as shown, and there are some smoothing options too, but let's not change them. Same as before, create a model, initialize into a simulation, run. This time directly for 12 days so that we can compare with the last plot -```julia -julia> model = ShallowWaterModel(;spectral_grid, orography, initial_conditions); -julia> simulation = initialize!(model); -julia> run!(simulation,n_days=12,output=true) -Weather is speedy: run 0003 100%|███████████████████████| Time: 0:00:35 (79.16 years/day) + +```@example howtorun +model = ShallowWaterModel(;spectral_grid, orography, initial_conditions) +simulation = initialize!(model) +run!(simulation,n_days=12,output=true) ``` -This time the run got the id "0003", but otherwise we do as before. -![Galewsky jet pyplot](https://raw.githubusercontent.com/SpeedyWeather/SpeedyWeather.jl/main/docs/img/galewsky_nc_12days_mountains.png) +This time the run got a new run id, `0002` in our case, which you can 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 howtorun +ds = NCDataset("run_0002/output.nc") +time = 49 +vor = Matrix{Float32}(ds["vor"][:,:,1,time]) +pcolormesh(lon,lat,vor') +xlabel("longitude") +ylabel("latitude") +tight_layout() # hide +savefig("galewsky3.png") # hide +nothing # hide +``` +![Galewsky jet pyplot](galewsky3.png) Interesting! The initial conditions have zero velocity in the southern hemisphere, but still, one can see some imprint of the orography on vorticity. You can spot the coastline of Antarctica; the Andes and @@ -191,10 +174,71 @@ probably not surprising! The life of every SpeedyWeather.jl simulation starts with a `SpectralGrid` object. We have seen some examples above, now let's look into the details -```@docs -SpectralGrid +```julia +help?> SpectralGrid +search: SpectralGrid + + Defines the horizontal spectral resolution and corresponding grid and the + vertical coordinate for SpeedyWeather.jl. Options are + + • NF::Type{<:AbstractFloat}: number format used throughout the model + + • trunc::Int64: horizontal resolution as the maximum degree of + spherical harmonics + + • Grid::Type{<:SpeedyWeather.RingGrids.AbstractGrid}: horizontal + grid used for calculations in grid-point space + + • dealiasing::Float64: how to match spectral with grid resolution: + dealiasing factor, 1=linear, 2=quadratic, 3=cubic grid + + • radius::Float64: radius of the sphere [m] + + • nlat_half::Int64: number of latitude rings on one hemisphere + (Equator incl) + + • npoints::Int64: total number of grid points in the horizontal + + • nlev::Int64: number of vertical levels + + • vertical_coordinates::SpeedyWeather.VerticalCoordinates: + coordinates used to discretize the vertical + + nlat_half and npoints should not be chosen but are derived from trunc, Grid + and dealiasing. +``` + +Say we wanted double precision (`Float64`), a spectral resolution of T42 on +a regular longitude-latitude grid (`FullClenshawGrid`) with cubic truncation +(`dealiasing=3`) and 4 vertical levels, we would do this by + +```@example howtorun +spectral_grid = SpectralGrid(NF=Float64, trunc=42, Grid=FullClenshawGrid, dealiasing=3, nlev=4) +``` + +We don't specify the resolution of the grid (its `nlat_half` parameter) directly, +instead we chose a spectral truncation `trunc` +and through the `dealiasing` factor a grid resolution will be automatically chosen. +Here T42 will be a matched with a 192x95 regular longitude-latitude grid +that has 18240 grid points in total. For details see [Matching spectral and grid resolution](@ref). + +We could have also defined `SpectralGrid` on a smaller sphere than Earth, +or with a different vertical spacing +```@example howtorun +vertical_coordinates = SigmaCoordinates(0:0.2:1) +``` +These are regularly spaced [Sigma coordinates](@ref), defined through their half levels. +```@example howtorun +spectral_grid = SpectralGrid(;vertical_coordinates,radius=1e6) ``` +In the end, a `SpectralGrid` defines the physical domain of the simulation and its discretization. +This domain has to be a sphere because of the spherical harmonics, but it can have a different radius. +The discretization is for spectral, grid-point space and the vertical as this determines the size of many +arrays for preallocation, for which als the number format is essential to know. +That's why `SpectralGrid` is the beginning of every SpeedyWeather.jl simulation and that is why +it has to be passed on to (most) model components. + ## References [^1] Galewsky, Scott, Polvani, 2004. *An initial-value problem for testing numerical models of the global shallow-water equations*, Tellus A. diff --git a/docs/src/primitiveequation.md b/docs/src/primitiveequation.md index 90b04226e..aa6ff991c 100644 --- a/docs/src/primitiveequation.md +++ b/docs/src/primitiveequation.md @@ -368,7 +368,7 @@ equation above, then we can also write + \sigma_{k+\tfrac{1}{2}}(-\mathbf{\bar{u}} \cdot \nabla \ln p_s - \bar{\mathcal{D}}) ``` See also Hoskins and Simmons, 1975[^HS75]. These vertical averages are the same as required by the -[Surface pressure tendency](@ref) and in the [Temperature equaiton](@ref), they are therefore all calculated +[Surface pressure tendency](@ref) and in the [Temperature equation](@ref), they are therefore all calculated at once, storing the partial averages ``\overline{\mathbf{u}_k \cdot \nabla \ln p_s}`` and ``\bar{\mathcal{D}}_k`` on the fly. ## Pressure gradient diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 14dd3dff8..a9314b5d5 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -18,7 +18,7 @@ import Adapt: Adapt, adapt, adapt_structure # INPUT OUTPUT import TOML import Dates: Dates, DateTime -import Printf: @sprintf +import Printf: Printf, @sprintf import NetCDF: NetCDF, NcFile, NcDim, NcVar import JLD2: jldopen import CodecZlib diff --git a/src/dynamics/vertical_coordinates.jl b/src/dynamics/vertical_coordinates.jl index 7f5eca2ba..7c3ce76ac 100644 --- a/src/dynamics/vertical_coordinates.jl +++ b/src/dynamics/vertical_coordinates.jl @@ -12,6 +12,7 @@ end # obtain nlev from length of predefined σ_half levels SigmaCoordinates(σ_half::AbstractVector) = SigmaCoordinates(nlev=length(σ_half)-1;σ_half) +SigmaCoordinates(σ_half::AbstractRange) = SigmaCoordinates(collect(σ_half)) function Base.show(io::IO,σ::SigmaCoordinates) println("$(σ.nlev)-level SigmaCoordinates:") diff --git a/src/utility_functions.jl b/src/utility_functions.jl index 2e591af5c..ad2e48fbc 100644 --- a/src/utility_functions.jl +++ b/src/utility_functions.jl @@ -1,8 +1,7 @@ """ - true/false = isincreasing(v::Vector) - +$(TYPEDSIGNATURES) Check whether elements of a vector `v` are strictly increasing.""" -function isincreasing(x::Vector) +function isincreasing(x::AbstractVector) is_increasing = true for i in 2:length(x) is_increasing &= x[i-1] < x[i] ? true : false @@ -11,10 +10,10 @@ function isincreasing(x::Vector) end """ - true/false = isdecreasing(v::Vector) +$(TYPEDSIGNATURES) Check whether elements of a vector `v` are strictly decreasing.""" -function isdecreasing(x::Vector) +function isdecreasing(x::AbstractVector) is_decreasing = true for i in 2:length(x) is_decreasing &= x[i-1] > x[i] ? true : false