Skip to content

feature: Struct-of-Array like iteration interface

Florian Atteneder requested to merge feature-soa-iteration into develop

The Problem

Currently, computing fluxes or other quantities is done by writing functions that compute stuff inplace of the cache variable. Within each function one has to manually iterate the grid points and do the computations. This can be easily done using broadcasting operations (https://git.tpi.uni-jena.de/dg/dg1d.jl/-/blob/develop/src/projects/ScalarEquation/equation.jl#L46).

There is at least two issues with this approach:

  1. Broadcasting does not work well with inplace methods or unpacking multiple return values (the latter could perhaps be solved with StructArrays.jl, but I am not sure).
  2. Functions are not easily reusable, e.g. we always have to call them with a full Cache object and the results are always computed over the full mesh.

Both of these drawbacks appear when computing numerical fluxes. E.g. the Lax-Friedrich formula requires the computation of the flux for only two points on the mesh (for each interface):

f^\ast = \frac{1}{2} (f(u^+) + f(u^-)) + \frac{1}{2}(u^+ - u^-)

Furthermore, for systems of conservation laws the flux is a vector, so it should actually return multiple values for each call. Note: Right now we compute the fluxes needed for the LLF along with the normal fluxes, because we have chosen to collocate the interface boundaries with our grid. But this might not work anymore if we chose to use other numerical flux methods (Riemann solvers) or non-conforming grids. And it already fails to work when using boundary conditions.

The Solution

This MR tries to solve the problem by adapting an Struct-of-Array like approach. To this end we redefine the interface for flux computations, etc. to compute the quantity in question for a single point only. By making the functions accept and return only tuples of numbers, we can utilize StructArrays.jl to perform 'broadcasting inplace' over (possibly any range of) the mesh without allocations.

One advantage of the previous implementation was that it established a single source of truth for the variables that are used for the computation. To establish that also for the new function interface we add a macro to automate the broadcasting.

Lastly, I should mention that this approach is inspired by Trixi.jl's equation interface, although, we broaden the interfaces to allow for more flexibility.

Merge request reports