-
Notifications
You must be signed in to change notification settings - Fork 29
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Conserved quantity verifying #400
Comments
If I read your y axis correctly then you have julia> h1 = 6200+6.3875e9
6.3875062e9
julia> h0 = 5500+6.3875e9
6.3875055e9
julia> h1/h0
1.0000001095889468 It looks like it's growing "fast" but actually it's conserved to 1e-7 which I find very reasonable given there are still several sources of error
Do you just want to be sure that volume is conserved? Then that's ☝️ probably good enough, if you want to try and close the mass budget as best as possible I'd start by those points. |
Oh! I see! Thank you very much! I made a mistake, because I just look at the trend. If I want to
How can I adjust, especially with the grid? maybe I should the even size lon-lat grid. In addition, at #391, when making the initial conditions:
How can I determine this number 24 which is the size of the array. I think this have something with the
So, how can I decide the array size by trunc-number? In other words, how do I know the the array size of initial conditions in diffrernt trunc-number? |
I am interested in the energy budget, so I also want to know the default viscosity coefficient and how to modify it. I read the output's parameters.txt which seem don't mention the viscosity. |
What is the mathematical formula of the conservations you want to use? Once you multiply h with u I'd definitely include orography. And what happened to hv? There's no forcing by default, you can check that with
There's no viscosity in the simulation. There's hyper diffusion which removes energy at the smallest scales to keep to model stable. Read the documentation https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#diffusion for more information. While the operator is dimensionless, the associated time scale is by default julia> model.horizontal_diffusion.time_scale
2.4 hours. You can translate that back into a diffusion coefficient. But note that by default we use a power 4 Laplacian so it'll have units of m^8/s! |
It's all in the documentation! |
Orography is included by default in the shallow water model, which you can check with |
Can you provide a reference for this? I need to check whether the assumptions in that derivation are also given for our |
Quick note, I'd recommend to use directly the values from the model julia> model.planet.rotation
7.29e-5
julia> model.spectral_grid.radius
6.371e6 The omega you use is slightly off because the planet actually rotates in 24 hours more than around itself. This is because in 24 hours it also travels 1/365.25 around the sun. Check https://en.wikipedia.org/wiki/Sidereal_time Normally I wouldn't bother but because you are trying to close some budget you should use what's actually used inside the model. Bottom friction/drag by default is julia> model.drag
SpeedyWeather.NoDrag{Float32} <: AbstractDrag which you should always check as default setups may change from version to version. |
For a normal Laplacian do julia> horizontal_diffusion = HyperDiffusion(spectral_grid,power=1,time_scale=24)
HyperDiffusion{Float32}
├ trunc::Int64 = 31
├ nlev::Int64 = 8
├ power::Float64 = 1.0
├ time_scale::Float64 = 24.0
├ resolution_scaling::Float64 = 0.5
├ power_stratosphere::Float64 = 2.0
├ tapering_σ::Float64 = 0.2
├ adaptive::Bool = true
├ vor_max::Float64 = 0.0001
├ adaptive_strength::Float64 = 2.0
└── arrays: ∇²ⁿ_2D, ∇²ⁿ_2D_implicit, ∇²ⁿ, ∇²ⁿ_implicit
julia> model = ShallowWaterModel(;spectral_grid,horizontal_diffusion) You will also need to increase the |
Oh, thanks, I will use the the values from the model directly. And I check the bottom drag: model.drag
SpeedyWeather.NoDrag{Float32} <: AbstractDrag I think this means no drag. So why is the absolute angular momentum decreasing in the first two months? |
Hi @milankl, I am not very familiar with spectral method but interested here. Standard viscous term, I and @mini-DONG just want to close the budget as best as we could, so we now choose the simplest form of viscous dissipation. Another thing is that if we choose to output at daily frequency, does the default model output is instantansous or daily-averaged one? I guess this may also affect the budget. |
On time scale, see https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#Normalization-of-diffusion
|
Yes, but it doesn't really make a difference whether you use time scale times power 1,2,3 or 4 Laplacian on your side, does it? I don't know how you want to calculate that term, but to be consistent with the model you'll need to transform into spectral space and apply it there anyway, in which case you also run into the problem that the diffusion inside the model is applied implicitly, check https://speedyweather.github.io/SpeedyWeather.jl/dev/barotropic/#diffusion
Output is always instantaneous, the only exception is precipiation which is accumulated - but you do not use that anyway. Non-instantaneous values would give you a hard time anyway as any non-linear metric wouldn't close I guess: For example in your formulation you have h times u, which is not the same as h averaged over some time times u averaged over some time. |
I don't know. For these kind of budgets I would always try to replicate as best as possible what's done in the model. With angular momentum I do not have enough experience to say what's affecting that budget and what is not. We would need to play around with it. There's many numerics inside the model that make this a non-trivial exercise
|
Hi, I ran a shallow-water simulation earlier, and now I want to add forcing on that. If I think of the shallow-water as an ocean, and I want to put wind stress on the upper surface, what should I do? |
Hey @mini-DONG, sorry I never responded to that. We have a help?> JetStreamForcing
search: JetStreamForcing
Forcing term for the Barotropic or ShallowWaterModel with an idealised jet
stream similar to the initial conditions from Galewsky, 2004, but mirrored
for both hemispheres.
• nlat::Int64: Number of latitude rings
• latitude::Any: jet latitude [˚N]
• width::Any: jet width [˚], default ≈ 19.29˚
• speed::Any: jet speed scale [m/s]
• time_scale::Second: time scale [days]
• amplitude::Vector: precomputed amplitude vector [m/s²]
julia> forcing = JetStreamForcing(spectral_grid)
JetStreamForcing{Float32} <: Main.SpeedyWeather.AbstractForcing
├ nlat::Int64 = 256
├ latitude::Float32 = 45.0
├ width::Float32 = 19.285715
├ speed::Float32 = 85.0
├ time_scale::Second = 2592000 seconds
└── arrays: amplitude This part in the documentation gives you an idea of how to set this up https://speedyweather.github.io/SpeedyWeather.jl/stable/setups/#Polar-jet-streams-in-shallow-water |
Hi @milankl, I just move back to this issue, where previously I have tried some schemes to see several conserved quantities like:
I and @mini-DONG found that the evaluation of these integrations depend on the grids and schemes. So I would like to know:
We need to close the energy budget using speedy simulation. So we want to know how well speedy is doing this (not perfect for every metric). Considering previous discussions at here, there could be several other things affecting the accuracy. Please remind us if you feel it is necessary. |
On 1: The linear quantities you can calculate in spectral space, e.g. total mass. Linear here means that any variable that depends on space can be added but not multiplied with another variable that also depends on space. E.g. where layer thickness spectral_grid = SpectralGrid(trunc=31, nlev=1)
model = ShallowWaterModel(;spectral_grid)
simulation = initialize!(model) Now the julia> simulation.prognostic_variables.surface.timesteps[1].pres[1]
4634.5415f0 + 0.0f0im It's imaginary part is always zero, so you can julia> a = model.spectral_transform.norm_sphere
3.5449078f0
julia> real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a
1307.38f0 So with the chosen initial conditions julia> run!(simulation, period=Day(10));
Weather is speedy: 100%|████████████████████████████████████████████| Time: 0:00:00 (12.73 millenia/day)
julia> real(simulation.prognostic_variables.surface.timesteps[1].pres[1]) / a
1307.38f0 And in fact it cannot, because |
On 2: The Gaussian and Clenshaw-Curtis grids are exact, meaning that the transform has only floating-point errors but no transform error. For Gaussian grids that happens at a dealiasing factor of spectral_grid = SpectralGrid(Grid=FullClenshawGrid, dealiasing=3) The FullClenshawGrid is a regular longitude-latitude grid, which likely makes your analysis easier. With Gaussian grids, the latitudes are not equally spaced because they are Gaussian latitudes, meaning you'll get another complication when calculating global integrals. Instead of calculating the global integral by hand you can also just transform a grid back into spectral space and then use the So for julia> u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid
julia> v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid
julia> KE = zero(u) # allocate grid of same size
julia> @. KE = 1/2 * (u^2 + v^2) # with @. everything element-wise
julia> real(spectral(KE)[1])/a # to spectral space, get l=m=0 and normalise
160.22296f0 |
Wrapped into a function julia> function total_energy(u,v,η,model)
E = zero(u)
H = model.atmosphere.layer_thickness
Hb = model.orography.orography
g = model.planet.gravity
@. E = 1/2 * ((η + H - Hb)*(u^2 + v^2) + g*(η + H - Hb)^2)
real(spectral(E)[1])/model.spectral_transform.norm_sphere
end You'd have initially julia> simulation = initialize!(model)
julia> run!(simulation, period=Day(0)); # run for 0 days to propagate spectral initial conditions to grid
julia> u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid; # flat copy
julia> v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid;
julia> η = simulation.diagnostic_variables.surface.pres_grid;
julia> total_energy(u, v, η, model)
4.5406352f8 then after 100days julia> run!(simulation, period=Day(100));
Weather is speedy: 100%|████████████████████████████████████████████| Time: 0:00:01 (14.91 millenia/day)
julia> total_energy(u, v, η, model) # u,v, eta still point to the same mutable array
4.5268032f8 The model has lost about 0.3% in total energy. I find that acceptable, and I wasn't expecting anything more than that to be honest. |
That's awesome! I am trying now but a tech problem is how to calculate absolute vorticiy. I have:
it is not possible to sum them up |
Great, that works for me.
What's the difference between reduced and full grids? Does the output always use interpolation so that there are unavoidable errors (e.g., cannot get exact total energy as in the spectral grid) even the Gaussian and Clenshaw-Curtis grids you suggested? I also found that the outputs do not have grid points at poles (+90 or -90). Does that affect the calculation of total (potential) energy? |
Read on grids here https://speedyweather.github.io/SpeedyWeather.jl/dev/grids/#Implemented-grids Output is by default on the corresponding full grid of whatever reduced grid the simulation uses and using the same resolution parameter julia> RingGrids.full_grid(OctahedralGaussianGrid)
FullGaussianGrid So when using reduced grids, the output is interpolated, for full grids it's not (unless you specify a different resolution). Read here https://speedyweather.github.io/SpeedyWeather.jl/dev/output/#Example-2:-Output-onto-a-higher/lower-resolution-grid No grid in SpeedyWeather (and in any other spectral model afaik) has a grid point on the pole. It's not needed for an exact transform and, in fact, its FFT would be just the identity. If you do need a point on the pole you can average the pole-most ring, e.g. # create a random grid of given type (octahedral Gaussian) and resolution parameter nlat_half = 4
grid = rand(OctahedralGaussianGrid, 4)
rings = RingGrids.eachring(grid) # index range per ring
first_ring = rings[1] # ring around north pole
last_ring = rings[end] # ring around south pole
using Statistics
north_pole = mean(grid[first_ring])
south_pole = mean(grid[last_ring]) for full grids you may intuitively think of doing To answer your question: If you calculate the global integral via the spectral transform as I did it, no it doesn't matter that we don't have a point on the pole. If you calculate the integral based on grid points by simply adding them up with area weights then you'll have a much larger error from numerical integration anyway (because that's not a good quadrature rule https://en.wikipedia.org/wiki/Numerical_integration). For the Gaussian grids we use the Gaussian quadrature (hence the name) https://en.wikipedia.org/wiki/Gaussian_quadrature which is exact (rounding errors only) if you use at leas the default On that note, I'm imagining a
as we discussed them here. If you'd like to contribute I'm happy to co-create a pull request to the documentation based on everything you've learned through this issue. We have already a similar section on Geostrophy which shows how to use the gradient operators directly from the model on the model grids https://speedyweather.github.io/SpeedyWeather.jl/dev/speedytransforms/#Example:-Geostrophy |
Great. I just learn a lot on grids here. So even if I prescribe the full grid for I see your concern on the summation over grid points. But the offline diagnoses would be more popular in practice. I and @mini-DONG are trying to see how well the offline calculation is close to the online calculation. We can put all these into a doc section once we figure these out clearly. You can point us to an outline or a template (probably the Geostrophy example page), so that we can fill in step by step. |
No no. You get what you ask for. julia> spectral_grid = SpectralGrid(trunc=31, nlev=5)
SpectralGrid:
├ Spectral: T31 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid: 48-ring OctahedralGaussianGrid{Float32}, 3168 grid points
├ Resolution: 401km (average)
└ Vertical: 5-level SigmaCoordinates
julia> output = OutputWriter(spectral_grid, ShallowWater)
julia> output.output_Grid
FullGaussianGrid
julia> spectral_grid.Grid
OctahedralGaussianGrid if you try outputting directly on the julia> output = OutputWriter(spectral_grid, ShallowWater, output_Grid=OctahedralGaussianGrid);
ERROR: MethodError: no method matching get_nlon(::Type{OctahedralGaussianGrid}, ::Int64)
Closest candidates are:
get_nlon(::Type{<:SpeedyWeather.RingGrids.AbstractFullGrid}, ::Integer) suggesting you to use a full grid instead. (Maybe that error message could be clearer).
I know because people don't expect that you can do these things interactively with a model! I think that's a massive problem in the climate modelling community that we have such a strong divide between the model and its analysis. The model is too often a stiff blackbox with code that you cannot reuse, and so everyone invents their own workflow based on the data the model outputs. What a missed opportunity! Models should be libraries that allow you to do all these things directly! That's at least what I'm working towards with SpeedyWeather!!
I would create a pull request but then you cannot add to it. Can you create a dummy pull request? Then I'll add the outline in the docs and you can fill in the rest. That would be very much appreciated! |
Hi, thanks a lot for your help at #391 . Now, I have run a shallow water model like:
with initial conditions:
I have run 1 year, and the first output (i.e. initial conditions) and last output is like:
Like #323, I also want to verify the conserved quantity. I calculate the total thickness with python codes:
The plot of H:
H is growing fast, which doesn't seem reasonable. Or maybe I should use other way to calculate?
The text was updated successfully, but these errors were encountered: