Flux functions that depend explicitly on space

[ranocha]

This came up in #356: It would be nice to allow flux functions that depend explicitly on space x, e.g. for linear advection with variable coefficients. There are at least two possible approaches:

1. Just provide the space coordinate explicitly as argument to all flux functions (calcflux and all numerical fluxes)
2. Create a trait for this property, e.g. has_variable_coefficients(equations) = Val(false) as default value and use that for dispatch. This would be similar to our approach for nonconservative terms used in the MHD equations.

If we choose the first approach, we should benchmark carefully whether there are any performance issues with this.

[sloede]

I think this is a good suggestion in general. However, there are two questions that we might also want to consider when addressing this issue:

i) Should we also pass t in that case to be consistent with other methods, which IIRC always accept both x and t and never just one of them? ii) How to handle variable coefficients that are not a (simple) function of space or time, but that should be prescribed as part of the initialization? I am thinking about linearized Euler equations or acoustic perturbation equations, where you need to specify, e.g., mean density or speed of sound etc. in often a non-trivial manner. Thus re-calculation of the coefficients is typically either not feasible (for performance reasons) or not possible (if they are from an input file) - in that case, you need permanent storage for these additional coefficient variables.

Especially ii) is a lot tougher to answer than the original issue, I think, but we should have at least the scenario in mind before deciding which way to go.

P.S.: @ranocha I took the liberty of adding numbers to your suggestions to reduce ambiguity.

[ranocha]

i) Should we also pass t in that case to be consistent with other methods, which IIRC always accept both x and t and never just one of them?

Good point, we should definitely allow time-dependent coefficients.

ii) How to handle variable coefficients that are not a (simple) function of space or time, but that should be prescribed as part of the initialization? I am thinking about linearized Euler equations or acoustic perturbation equations, where you need to specify, e.g., mean density or speed of sound etc. in often a non-trivial manner. Thus re-calculation of the coefficients is typically either not feasible (for performance reasons) or not possible (if they are from an input file) - in that case, you need permanent storage for these additional coefficient variables.

That's also a good point. A solution might be to introduce a possible storage for auxiliary variables, which could also be used to store the primitive variables (see #54). Then, we can decide whether we want to pass the auxiliary variables as additional argument (might be cleaner) or whether we want to augment the vector of conserved variables when calling get_node_vars. In that case, we could use functions like update_coefficients!(semi, t) that will update the auxiliary variables if necessary and we won't need to pass the time explicitly to flux functions etc. These kind of auxiliary variables will also provide the necessary cache for parabolic problems when we want to store the gradients separately. The default choice would be to use nothing as auxiliary variables and functions like get_node_vars (or similarly named variants, get_node_aux_vars???) would need to be specialized to return nothing accordingly. We would also need a mechanism to decide how AMR should handle auxiliary variables

• do nothing (when storing things that need to be re-evaluated anyway, e.g. gradients, primitive variables)
• use AMR transfer operators as usual (for variable coefficients given as initial data)
• use something like update_coefficients!(semi, t) for coefficients given as functions

[ranocha]

On the other hand, it's probably just okay to use a number of auxiliary variables which can also be zero.

julia> Array{Float64, 3}(undef, 0, 4, 5)
0×4×5 Array{Float64,3}

julia> Array{Float64, 3}(undef, 0, 4, 5)[:, 1, 1]
Float64[]

julia> ntuple(n -> Array{Float64, 3}(undef, 0, 4, 5)[n, 1, 1], 0)
()

julia> using StaticArrays

julia> SVector{0,Float64}(ntuple(n -> Array{Float64, 3}(undef, 0, 4, 5)[n, 1, 1], 0))
0-element SArray{Tuple{0},Float64,1,0} with indices SOneTo(0)

This kind of auxiliary variables could live for example in semi.cache.elements.

[sloede]

Having such auxiliary variables together with semi.cache.elements would make sense to me as well, since this is exactly where I put such variables when I implemented this features in my previous code ZFS :-) However, one challenge to consider is what to do when you have several orthogonal features that all require theses auxiliary variables. E.g., you want to store primitive variables because you need them in the source terms where they would be expensive to calculate (for whatever reason), while at the same time you need mean values for your linearized system. In that case:

• Do you use a single data structure? But how to properly access the different categories of data then, e.g., in get_node_aux_vars?
• Do you use one data structure per category? But how to pass them to different functions then, given that we want to be able to keep using node-based functions?

Regarding the AMR mechanism: I think these are exactly the three types of update mechanisms we need to support.

[ranocha]

What kind of applications of these auxiliary variables will we have in the (near) future? Right now, I don't think we will have, say, hundreds of auxiliary variables. In that case, it shouldn't be a performance issue to store them in one place as a single vector. The order and kind of auxiliary variables will depend on the equations, so they will be specified and used correctly using multiple dispatch. On the other hand, we could also distinguish the three different kinds of auxiliary variables determined by their AMR behavior and store them in three different vectors with a convenience function returning all of them together.

[sloede]

Hm, you are right. Let's KISS: one storage location for all auxiliary nodal variables, and everything else is trait-based. This looks to me like the least intrusive approach for a full set of features, and if we need something more complex in the future, we'll re-evaluate. This way, each system of equations can support these variables, and it just depends on the way the equation is created. As long as the total number of auxiliary variables will become part of the equation type, we can dispatch/trait anywhere it's required.

[sloede]

We had a discussion at this week's Trixi meeting about how to implement such variable coefficients for the acoustic perturbation equations (see related PR #450). In this case, the coefficients represent time-averaged quantities of the mean flow field, such as mean density, mean velocity, or mean speed of sound. Several possible approaches were considered:

1. Treat variable coefficients as if they were state variables and store them in u, after the actual conservative state. Then set their fluxes to zero such that they are never updated. This is probably the least intrusive approach, but it is not very clean and has a performance overhead. Also, it does not allow to easily treat the coefficients specially where it might be required, e.g., for AMR.
2. A modification of the first approach, where we still store variable coefficients in u but modify the volume integral, surface integral, and other relevant routines to just ignore the additional coefficients to reduce computational overhead.
3. Modify the existing solver such that it treats all problems as variable coefficient problems. That is, variable coefficient data is passed to each function that potentially might need it. For equations that don't need it, the passed data is just an empty tuple/SVector. The downside here is that we add code to many routines that every user has to take into account but only a subset actually uses.
4. Create a new solver for such variable coefficient problems that allows to store the coefficients separately and pass them explicitly into every function that may require them. This has the advantage that we have a lot of flexibility regarding special treatment of such variable coefficient problems, with the downside that it creates another implementation variant.

There was no final conclusion on what the best approach would be for us in Trixi. Instead, we decided to approach this iteratively by starting out with a simple and easy-to-implement approach, and upgrade it in case we face serious issues or once variable coefficient problems start playing a more important role. In case of #450, we specifically decided to use the first approach for now (store and treat the coefficients as additional conservative quantities with zero-valued fluxes).

[DanielDoehring]

Has anything changed in the past 1,5 years regarding this issue? For instance, what about equation dependent function signatures, i.e., for flux? Then one would only read and pass $t$ and $\boldsymbol x$ if actually needed.

Reason I am bringing this up again is that for instance the advection equation in form $$u_t + a(x) u_x = 0$$ with $$a(x) = \begin{cases} 1 & x \leq x_\max/2 \ 1000 & x_\max/2 < x \leq x_\max \end{cases} $$

is a convenient way to set up "locally stiff" problems.

[jlchan]

I believe that the approach @andrewwinters5000 took for the shallow water equations with variable bathymetry could be adapted here. If a(x) doesn't depend on time, it can be represented in the same basis as the solution, stored as a "dummy" field in the ODE, and extracted within the flux function to compute a(x) * u(x).

A similar approach would be to simply store a vector of values for a(x) within the equations struct and extract it each time to compute the flux. The downside to this is that it forces equations to depend on the discretization, which isn't as clean.

Your suggestion about creating equations which explicitly depend on space/time should also work, but would require a few more changes to the solver. If you'd like to add this, we can help you through the process.

[ranocha]

The workaround @jlchan mentioned is described in one of our tutorials: https://trixi-framework.github.io/Trixi.jl/stable/tutorials/adding_nonconservative_equation

Sadly, we did not have the time to design a proper solution to this issue.