Skip to content
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

Vorticity initial conditions #153

Closed
wants to merge 5 commits into from
Closed

Vorticity initial conditions #153

wants to merge 5 commits into from

Conversation

maximilian-gelbrecht
Copy link
Member

@maximilian-gelbrecht maximilian-gelbrecht commented Sep 26, 2022

I first wanted to do it by requiring to hand over an instance of PrognosticVariables, but for that you have the problem again that PrognosticVariables are not defined at the time the inputs are handles. Changing that would need larger changes, so I am thinking handing over initial conditions as Vector{LowerTriangluarMatrix} is fine

@maximilian-gelbrecht
Copy link
Member Author

maximilian-gelbrecht commented Sep 26, 2022

There's a Julia subtlety that I don't get, but noticed while doing this PR:

(LowerTriangularMatrix{Complex{NF}} <: LowerTriangularMatrix) == true

but

(Vector{LowerTriangularMatrix{Complex{NF}}} <: Vector{LowerTriangularMatrix}) == false

I know that for parametric types the type without the parameter is automatically the supertype (and abstract) to all discrete types, but this doesn't work as I expected for the Vector anymore.
Any idea why?

One can do

(Vector{LowerTriangularMatrix{Complex{NF}}} <: Vector{<:LowerTriangularMatrix}) == false

though

@white-alistair
Copy link
Member

Isn't this precisely the invariance of Julia types discussed here?

@maximilian-gelbrecht
Copy link
Member Author

Ah yes, you are right. I also did encounter this before, but somehow wasn't thinking about it in that way.

@milankl
Copy link
Member

milankl commented Sep 26, 2022

I very much like the idea to have the initial conditions easily choosable with arrays that are defined in the REPL before run_speedy is called. However, I'm not sure about the current idea how this PR suggests to implement it. The reason is that conceptually Parameters should contain parameters and not entire arrays that contain data. This may sound pendatic, but given the discussion we had about also including the parameters into the restart file, suddenly including a potentially massive array doesn't sound neat. Therefore first question, what's wrong with

p,d,m = initialize_speedy(kwargs...)
p.layers[1].leapfrog[1].vor .= ...   # change vorticity, or any other prognostic variable
time_stepping!(p,d,m)

where the second line could also be replaced by something like

set_vorticity!(p,some_array...)
set_divergence!(p,some_other_data...)
...

which could take LowerTriangularMatrix or some gridded data (which is then transformed to spectral space), scalars, etc. The reason I'm suggesting this is also that it feels more modular, as the function how to get the data into the prognostic variables struct doesn't have to be buried in initialize_speedy!, and people can easily define their own extensions.

lmax,mmax = size(vor_new) .- (2,1)

vor_ic_trunc = spectral_truncation(layer_ic, lmax+1, mmax)
vor_new .= convert(LowerTriangularMatrix{Complex{NF}}, vor_ic_trunc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that we still have #135 open, meaning that vor_new .= some_array may not do what you expect it too. To avoid this broadcasting issue in the first place I've already defined copyto! for :LowerTriangularMatrix which can handle the type conversion in convert implicitly.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is just copied from the initialize_from_file! function. Then that should also be looked at, but we can also replace that function with set_vorticity! etc as well

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, good point!

progn, diagn, model = initialize_speedy(initial_conditions=ic,model=:barotropic)
for layers in progn.layers
for leapfrog in layers.leapfrog
@test all(leapfrog.vor .== 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given the open broadcasting issue related to the comment above, note that this is a redundant test: Replacing zeros with zeros would pass even if nothing is actually done.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here as well, that's essentially the same test as already done for the other initialisations

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Haha, boooh past me setting up shitty tests 😆

Comment on lines +83 to +88
@test all(leapfrog.vor .== 0)
@test all(leapfrog.div .== 0)
end
end
@test all(progn.pres.leapfrog[1] .== 0)
@test all(progn.pres.leapfrog[2] .== 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.

@@ -109,7 +109,8 @@ The default values of the keywords define the default model setup.

# INITIAL CONDITIONS
seed::Int = abs(rand(Int)) # a random seed that's used in initialize_speedy for the global RNG
initial_conditions::Symbol=:barotropic_vorticity # :rest, :barotropic_vorticity or :restart
initial_conditions::Union{Symbol, Vector}=:barotropic_vorticity # :rest, :barotropic_vorticity, :restart
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another problem I see here is that with output=true the parameters struct is printed to parameters.txt such that one can easily trace back what parameters were used in a given simulation. I'm not sure what happens if a big array is part of the parameters struct. This also applies to setting the vertical levels manually, which I haven't checkedn. Those ones, ideally, should be written into that .txt file (as just a vector) but entire vectors or matrices obviously not...

@maximilian-gelbrecht
Copy link
Member Author

Mh, initial conditions do belong to a well defined differential equation problem, so I am not fundamentally opposed to including them in the Parameters (as this is also currently the only way input arguments are handled).

However, I also like your idea of having flexible set_vorticity!(p,some_array...) etc. functions. This does look like a clean solution.

@maximilian-gelbrecht
Copy link
Member Author

I am working on a new draft for this

@milankl
Copy link
Member

milankl commented Sep 27, 2022

Thanks! I'm busy this week, but happy to review code.

I absolutely agree with you that technically the initial conditions should be part of the problem. Having said that, at the moment we do set up a fully defined "speedy problem" with initialize_speedy, but split into 3 groups P,D,M: prognostic variables/initial conditions P, diagnostic variables/work arrays (whose initial state is irrelevant) D, and the model struct M that contains the parameters. So the idea behind is that one should be able to mix and match prognostic variables and model structs without running into something being inconsistently defined. For example, when using initial_conditions=:restart then there's that and the restart id stored in model, and the initial conditions in prognostic variables. If we now include the initial conditions also in M then they are defined twice, in P and in M.

In order to be consistent I'd probably suggest to allow any Symbol for initial_conditions, such that when someone does initialize_speedy(initial_conditions=:hurricance_sandy) it will just fall back to setting up things from rest, one can then use set_vorticity! to set things to hurricane sandy, such that we have the initial conditions stored in P and a matching description :hurricane_sandy in M. How does that sound? That also means if we export both P and M then the problem is fully defined again, as one can recreate D from M.

@maximilian-gelbrecht
Copy link
Member Author

This is closed in favour of #154

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants