Towards standalone packages + ColumnFields

[maximilian-gelbrecht]

Ok, I really quickly got the basics, but somehow this broke GPU broadcasting. So I'll have to finish this another time.

[milankl]

Looks great already!! I like it how you tried to reuse some of the field functionality but I guess in the end there's just some duplication we can't avoid or it would be too cumbersome to do. Big question for me is just though how we handle arithmetic operations between a field and a column field, like should +, *, etc just be allowed? With A+B one expects to allocate another array so +50% in memory what A,B already use. For field + column_field naively this would be +100% as the tranpose scratch memory also has to be allocated. I say "naively" because technically you can allocate the scratch and then accumulate field into it, so +50% also here. But then would we say field + column_field returns a Field so has precedence in type promotion?

[maximilian-gelbrecht]

I fixed the broadcasting issue.

Wrt arithmetics: What do we need for the model? I would probably keep this quite minimal now, and just see what we need later. E.g. with multiplication it's also not clear to me immediately how we should define it.

One thing I didn't adapt is the interpolation, I think that's just not necessary. But this means, I'll probably have to restrict the interpolation to only accept Field and not AbstractField anymore.

[milankl]

Wrt arithmetics: What do we need for the model? I would probably keep this quite minimal now, and just see what we need later. E.g. with multiplication it's also not clear to me immediately how we should define it.

I don't think we need anything for the model because the idea isn't to mix field and column fields but use the latter only after explicit tranposes. We could completely skip the idea of +, * etc of field and column field to be honest. Makes our live easier and it would only be for convenient interactive playing around anyway

One thing I didn't adapt is the interpolation, I think that's just not necessary. But this means, I'll probably have to restrict the interpolation to only accept Field and not AbstractField anymore.

Yes, I would do that too!

[milankl]

Having said this, we could check whether another tranpose ... tranpose block could be put around the vertical advection. WENO is currently quite costly and it might be cheaper if the memory was reordered.

[maximilian-gelbrecht]

Tests are passing now on my laptop, so this is almost ready. I wouldn't add a lot of arithmetics, at least not now.

I am just not sure about the exact way to do the transpose_unsafe!, I took over your sketch from the issue, but I had to add a reshape, which I think made the whole idea with the scratch memory pointless? This is also related to your comment about reshaping in the constructor.

I don't have further time today though, I'll look at it again later.

[milankl]

I am just not sure about the exact way to do the transpose_unsafe!, I took over your sketch from the issue, but I had to add a reshape, which I think made the whole idea with the scratch memory pointless? This is also related to your comment about reshaping in the constructor.

Did this in f0125f1 I hope I didn't break anything! reshape is a view on the same data just with different dimensions as long as the length stays the same! So regardless the size of scratch, reshape(scratch, size...) will reshape the scratch memory how we need it -- as long as the length fits!

[milankl]

Okay I understand that transpose(transpose(vector)) is defined as unity, so a 3-element vector would become a 1x3 matrix but the singleton dimension is dropped in the backtranspose

julia> transpose(transpose([1,2,3]))
3-element Vector{Int64}:
 1
 2
 3

julia> transpose([1,2,3])
1×3 transpose(::Vector{Int64}) with eltype Int64:
 1  2  3

but the original size here is somehow stored, see how it says Vector despite being 1x3 ... we don't do that. In that sense maybe what we currently do, not removing the singleton trailing dimensions is best

julia> field
24-element, 4-ring OctaminimalGaussianField{Float32, 1} on Array on Main.RingGrids.Architectures.CPU{KernelAbstractions.CPU}
 0.0f0
...
 0.0f0

julia> transpose(field)
1×24, 4-ring ColumnField{Float32, 2, Matrix{Float32}, OctaminimalGaussianGrid{Main.RingGrids.Architectures.CPU{KernelAbstractions.CPU}, Vector{UnitRange{Int64}}, Vector{Int64}}}
 0.0f0  0.0f0  0.0f0  0.0f0  0.0f0  …  0.0f0  0.0f0  0.0f0  0.0f0  0.0f0

julia> transpose(transpose(field))
24×1, 4-ring OctaminimalGaussianField{Float32, 2} on Array on Main.RingGrids.Architectures.CPU{KernelAbstractions.CPU}
 0.0f0
...
 0.0f0

seee how 24-element transposes to 1x24, transposes to 24x1. Hence we have

julia> field == transpose(transpose(field))
false

which is consequent but counter intuitive. But it only arises with Field2D which transpose isn't made for anyway, so I'd be happy to just ignore that and note it here. One can always transpose! in-place and reuse the old reference again without the trailing singleton dimension.

[maximilian-gelbrecht]

Ok... SpeedyInternals is there.

So, are we happy with this now? Then, the docs would need to be adjusted to reflect this new structure. Then merge, then register.

[milankl]

I can do a review later?

[milankl]

Just checked, I can develop SpeedyWeather and any submodules independently in VS code on this branch with no package or scope issues. Fantastic!

[maximilian-gelbrecht]

I revised the readme and the leading paragraph of the docs of the sub-packages. Anything else that we to adjust in the docs?

[milankl]

Do we want to add a README.md at src/LowerTriangularArrays/README.md and similar for other packages too? Maybe they are even required by the registrator?

[maximilian-gelbrecht]

Yeah, we could do that. I am also not sure what exactly the registrator even requires for packages within other repos, the manual doesn't even mention that it's possible to register submodules at all.

[milankl]

Sorry not done yet with the readmes ...

[maximilian-gelbrecht]

Ok, I finished the readmes now.

[maximilian-gelbrecht]

I think what's still missing for registering the sub packages is a license file in each subfolder. Our current license is including e.g. the copyright notice for the fortran speedy parametrization. I don't want to copy that to every submodule, but we have to keep the MIT licence obviously.

Shall we just use a standard MIT licence with
"The SpeedyWeather.jl Contributors for SpeedyWeather.jl" in each of the submodules?

[milankl]

Or is that the moment where we think about using EUPL? I would like to have some copyleft aspect in our license and even if it's a weak one

[maximilian-gelbrecht]

We would need the agreement of all contributors to the submodules though, right?

[bgroenks96]

Yes, unless you get the consent of all contributors, you need to keep it under the same license. I would recommend keeping MIT on everything in the Speedy monorepo for now.

[milankl]

Okay happy with that

[milankl]

Merge now? 😄

[bgroenks96]

I think you should do the honors @milankl!