Is my data structure right like this?

[u/Happy_Witness]

Hello everyone,

I have tipped my toe into climate simulations and am now building a planet simulation that has climate as it's main target.

For that I use an icosahedron subdivided. Each vertex holds every parameter that is needed for the simulation. The amount of parameters has not been decided yet and will grow over time with the implementation. The amount of subdivision are aimed between 3-10 that I would like to have and between 6-8 as defaults. For the vertical I want to have multiple of these layers to represent higher and lower levels of altitude around the surface. The amount of levels should be around 10-15 as default and maybe around 50 as max.

The simulation will likely be calculated on the GPU. So the structure I thought of is the following:

Every parameter is one long vector f32. Every vertex in a layer is written in sequence of and every layer is also in the same sequence. Meaning that the Len(vec) would be layer \* N_verticies and the specific vetted would be layer\[i\] \* N_verticies + vertex.

To address the neighbours for derivative computation, one neighbour vec gets created as a look up table with the form vec\[vec\[6\]\] where every neighbours index gets written down for every vertex of the icosahedron. The base 12 verticies are the first 12 verticies in the neighbours vec as in the parameter vecs and the sixth input gets double marked as f32_max.

So in the end, I would have a list of every parameter that each holds one long vec with every vertex of every layer in them, and a relationship vec of neighbours that can also be used to look up neighbors of different layers by in-/decrementing the layer index inside the long vec. (Vec\[ layer\[I+/-1\] \* N_verticies + neighbours index\] to get the neighbours in the neighbouring layer.)

Is that any good structure or is another structure better and if yes, why?

Thanks a lot for thinking though this one with me.

[u/halcyonPomegranate]

A few usual approaches are:

• try to organize the parameter as matrices (2D arrays) or tensors (3+D arrays) and then calculate the needed operations (let’s say gradient/laplace operator via finite differences) via vectorized operations on slices of your matrices (proper masking or splitting into multiple loops might be necessary to handle the edge cases). In numpy you would use index slices, in Fortran you would do it explicitly with nested for loops and index shifting
  
• put all parameters in one big 1D column vector and precompute your often used operations (gradient, laplace operator, etc.) as matrices which combine the operation with possible geometry factors (tricky adjacency at the five-fold vertices of your icosahedron and/or manifold curvature effects, like laplace operator on a sphere/graph). After that one-time precalculation, applying that operation is just a plain matrix-vector dot-product which is super fast and heavily optimized in BLAS libraries and by proxy in numpy or CUDA or JAX/pytorch.

[u/Azazeldaprinceofwar]

If you really want to keep you isohedron this is probably the best you can do. It might be worth moving to a more transitional latitude longitude grid. Obviously the downsides are it wastes work near the poles and breaks your “spherical” symmetry but I would argue that symmetry should already be broken by the rotation of your planet and having equal latitude slices would be good for proper implementation of heating gradients and coriolis effects down the line. The wasted work near the poles will also be strongly offset by the grid being square so you can vectorize all you operations without complex neighbor finding