How to deal with source terms?

[benegee]

For simulations where source terms (aka right hand side, forcing, ...) change over time, we need a way for an external application to provide the respective values.

Trixi.jl expects a method source_terms(u, x, t, ...), which is called in every time step for every DG node, i.e.

• time t
• coordinates of the DG node x
• current values of conservative variables u(x, t)

are given and the values of the source term S(u(x,t), x, t) need to be returned.

The method source_terms(u, x, t, ...) can simply be defined in a libelixir.
At his point there is a lot of flexibility. One could think about caching, interpolation onto DG nodes, ...
But in the end we need a mechanism therein that allows to retrieve whatever values there are from the calling external application.

We implemented two prototypes for demonstration:

1. Create a global data storage which can be written to via API functions: RFC: source terms via global vector storage #173

   • The external application allocates memory and libtrixi stores pointers to it
   • The external application computes source terms (whenever it wants) and stores the results in this memory
   • libtrixi reads values from this memory every time source_terms is called
   • To be able to compute the source term, the external application possibly needs to separately get t, x, u
   • A mechanism is required to map coordinates x to an index, needed to access the memory at the correct position

2. Let the user call her own external functions from a libelixir: RFC: Source terms via callback #172

   • The user has to provide a function similar to source_terms and make it available via a shared library
   • The function is then called each time source_terms is called.

[sloede]

I think creating PRs to evaluate different ideas is a great idea, as it is (usually) much easier to reason over concrete code than to just talk about concepts in an abstract manner!

Yet, I think it would be helpful (for the audience, but maybe for you too) to go into more detail in your description above (maybe directly amend the original post):

• Temporal properties: How often will the source terms change? Once every N steps, where N being variable? Or always with N == 1? Or do we need to take into account sub-step updates (i.e., per Runge-Kutta stage)
• Spatial properties: What is the original set of spatial coordinates at which the source terms are provided? Is it always the cell/element centers? Is it the DG nodes? Is it some cell-local distribution (but not the DG nodes)? Is the position of the source terms fixed with respect to the mesh?

I feel some additional information would make it easier to decide - by the process of elimination - which approach could be viable and which one isn't.

[KerstinHartung]

Thanks for preparing this comparison and the code examples!

Right now I would say that from the application perspective the first approach fits better.

• The sharing of memory can also be used to transfer information on the initial state but wouldn't be covered by the second approach.
• The source terms will be computed after a new state was received from libtrixi. After this is completed the next iteration of libtrixi (reading the source_terms) will be called.
• The last two points of the first approach (information on t, x and u needs to be shared; the agreement on common indexing (very important!)) are independent of the way libtrixi deals with the source terms.
• However, if the first approach has benefits from the Trixi.jl side (at least there seemed to be fewer code modifications) we could also work with that.

Regarding the other questions:

• I think we will mostly use N==1. Maybe not all processes contribute to the source terms each time step but once the state is updated MESSy will respond to that. Right now I think that sub-stepping is not necessary.
• So far in GCMs, most variables are defined at the cell centers but wind components are often defined at the cell edges. However, if this is not required by the dynamical core (and the tracer transport) then MESSy doesn't require this either. We can discuss further if sharing of data on the DG nodes would be beneficial. Whatever the approach the position of the source terms would be fixed with respect to the mesh.

[benegee]

Thanks for your feedback!

The sharing of memory can also be used to transfer information on the initial state but wouldn't be covered by the second approach.

Correct! Indeed, Trixi.jl uses a similar callback for the initial condition and this could be fed from the global database in an analogous manner. In case of the second approach another external function would be required.

The last two points of the first approach (information on t, x and u needs to be shared; the agreement on common indexing (very important!)) are independent of the way libtrixi deals with the source terms.

That's right, but to clarify: In the second approach you get all the information you need as parameters in every call to your external function. This function would be called once for every single point where source terms need to be computed. In the first approach you would have to get the coordinates and the current state (as vectors!) in advance and then do a batch-process to fill the source terms vector.

So, in a way the data storage approach is more generic, but it comes at the cost of organizing the required information in advance.

[benegee]

Regarding cell indexing: For the cell -> index mapping, I think it would be best to follow t8code's indexing (based on the Morton curve).

Indexing points withing cells, e.g. DG nodes, would be another thing and require a convention. However, I am not sure yet where we will ultimately need this. In this simple examples here the source term is given as a mathematical function which can compute values for arbitrary x. In practice I assume that we will compute source terms from the state point by point. So the indexing just needs to be consistent.

[sloede]

Thanks for the feedback @KerstinHartung! I think while it is not the most elegant approach, having a dedicated and shareable "data container" for source terms, IC data etc. is probably the most straightforward approach. The main downside here is that it will have to be adjusted accordingly with the number of elements upon AMR. But if that turns out to be a major pain point, we could also switch to a more sophisticated approach later.