Repository navigation
How to deal with source terms? #174
Description
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:
-
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_termsis 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
xto an index, needed to access the memory at the correct position
-
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_termsand make it available via a shared library - The function is then called each time
source_termsis called.
- The user has to provide a function similar to
Activity
sloede commentedon Mar 1, 2024
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):
Nsteps, whereNbeing variable? Or always withN == 1? Or do we need to take into account sub-step updates (i.e., per Runge-Kutta stage)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 commentedon Apr 5, 2024
Thanks for preparing this comparison and the code examples!
Right now I would say that from the application perspective the first approach fits better.
source_terms) will be called.t,xanduneeds to be shared; the agreement on common indexing (very important!)) are independent of the way libtrixi deals with the source terms.Regarding the other questions:
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.benegee commentedon Apr 8, 2024
Thanks for your feedback!
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.
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 commentedon Apr 8, 2024
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 commentedon Apr 8, 2024
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.