Perturbation formulation and variable-coefficient problems

[tristanmontoya]

As discussed at the first TrixiAtmo meeting two weeks ago, one possible strategy for implementing a perturbation formulation is to store a reference solution field in the cache, and give the flux function access to the cache, similarly to what I propose for the covariant form (which requires accessing geometric information from cache.elements) in #31.

This would be quite easy to add to my covariant form implementation, as the cache is already a parameter in the flux and source terms there. For more standard Cartesian formulations, this would require new kernels to be defined, which could be used more generally as a way of dealing with variable-coefficient problems. While passing the cache into the flux is a bit ugly, the "classical" way of treating variable coefficients as dummy variables (i.e. solution variables to which no flux is applied) is arguably even more of a hack, as is the alternative of subtracting the steady-state right-hand side at every residual evaluation.

I'd be happy to hear any ideas people may have about how to best implement a perturbation form (or variable coefficients more generally) and if they'd find this cache approach to be useful as an alternative to the "dummy variables" or "subtract the steady state" techniques.

[sloede]

Regarding the variable coefficient problem: There is no separate issue for it AFAIK, but there's been at least one lengthier discussion on the various options in trixi-framework/Trixi.jl#358 and this post gives a nice summary, while this one addresses the issue of AMR handling.

IMHO, the decision hinges to some extent on the question whether this feature is needed for (a) a solver variant that currently looks interesting to try out, or whether (b) it is clear that this will be needed for sure in the future no matter what.

If it is (a), I suggest to go for the current hack for now, since it causes the least amount of overhead for you to implement it. If, however, you are already quite sure it is (b), now would be a good time as ever to tackle this properly. In the latter case, I suggest to try to set a dedicated meeting with everyone involved to figure out the best way forward, since it is a likely to be a fundamental structural change.

On a personal note, I would try to refrain from passing cache into flux functions. It will forever merge solver implementation details with seemingly unrelated equation-specific functions, without even given access to everything a user might want to have (such as information about time and space). My gut tells me there is a better way to handle it (but of course I could be wrong)

[benegee]

Thanks @tristanmontoya for bringing this up and @sloede for hinting at what you already discussed.

In my opinion it is b) and I was already quiet concerned about how to do it without touching Trixi.jl.

I do not understand the details yet. To me trixi-framework/Trixi.jl#358 (comment) sounds good. Or the more general suggestion 4 in trixi-framework/Trixi.jl#358 (comment). But can this be done without passing cache around?

[tristanmontoya]

Thanks for pointing me to these previous discussions. I would probably classify this as (b) a feature that will definitely be needed. Implementing equations in perturbation form is a fairly standard approach in atmospheric modelling, and it would be good to have it implemented "natively". At the very least, this has already come up many times in the past few months.

I would prefer to take the approach of computing the nodal values of variable coefficients in initialization and then giving the flux access to these, rather than passing x and t to the flux function.

@sloede instead of passing the entire cache, would you be happy with passing the nodal values of the needed auxiliary variables or variable coefficients (in the covariant form, this would be the metric tensor components, and in the perturbation form, this would be the mean values) into the flux?

@amrueda this would affect the PR (#31) you are currently reviewing if I make this change, as that one currently passes around the cache. Either we could merge #31 first and then make changes, or we could hold off until this is figured out.

[benegee]

I would rather take the approach of computing the nodal values of variable coefficients in initialization, and then giving the flux access to these, rather than passing x and t to the flux function.

Just want to support this. For the examples involving the hydrostatic moist base state we have to solve an ODE globally to compute the initial condition.

[ranocha]

So we will go with some auxiliary variables. Could you please add this issue to the next Trixi.jl meeting agenda?

[tristanmontoya]

@ranocha I added it in. Is this something that you would want to have in Trixi.jl, or kept as a separate implementation in TrixiAtmo.jl and possibly added to the main package later?

[sloede]

IMHO a good solution that can live in Trixi.jl would be optimal. This would certainly maximize the number of potential applications that can benefit from this (e.g., "standard" variable coefficient problems, @benegee's applications, linearized equations with spatially varying coefficients etc.) and limit the amount of code duplication we have in TrixiAtmo.jl vs. Trixi.jl.

If we cannot find a good consensus for Trixi.jl base, only then would I implement it in TrixiAtmo.jl, gather some experience, and then decide later whether to upstream it or not.

IIRC, the two biggest conflicting priorities in terms of how to handle variable coefficients were these two:
a) Having (yet) another solver variant (curse of dimensionality) vs. everyone has to read/understand variable coefficients even if they don't use it.
b) AMR: interpolating coefficients vs. re-initializing them

Arguably, b) is easier to fix since it's not strictly either/or, while a) is the tough one that essentially prevented us from having first-class support in the past.

[tristanmontoya]

That makes sense. I'd be happy to work on the Trixi.jl implementation, but perhaps it makes sense if I did a minimal proof of concept where I use the same approach (auxiliary variables) for the covariant formulation in TrixiAtmo.jl first (replacing the current implementation in #31 where I pass around the cache), since that requires its own solver variant anyway?

Then, based on what works best there, we could start a PR in Trixi.jl. My (biased) preference would be to internally pass variable coefficients everywhere with an interface like flux(u, aux_vars, orientation_or_normal_direction, equations), but still allow users to just define flux(u, orientation_or_normal_direction, equations) if no variable coefficients are used, thus allowing for a simple interface if desired. Repeating the code for every part of rhs! for every solver type so as to have versions with and without variable coefficients is not ideal, in my opinion.

Edit: The minimal proof of concept is implemented now in #47.

[benegee]

Happy to help if possible @tristanmontoya

[tristanmontoya]

@benegee We have now merged #47, which adds the additional cache field for auxiliary variables. While we only use it for the covariant form in that PR, it could be used more generally by modifying rhs! to dispatch each kernel on the have_aux_node_vars trait, similarly to have_nonconservative_terms. The corresponding equation sets would need to be defined so as to take aux_vars as an additional argument in the fluxes, and init_auxiliary_node_variables! would have to be specialized for each equation to initialize the auxiliary variables appropriately.

[benegee]

Thanks @tristanmontoya and @amrueda for your work! This looks very promising!

I think it would be best to have a nice example using the aux variables when bringing this to Trixi.jl. I am thinking about the warm bubble with hydrostatic background state, but still need to wrap my head around the splitting.