add face-centered velocities to `state_fc_new_`/`state_fc_old_`

[BenWibking]

Describe the proposal Add the face-centered velocities returned by the Riemann solver to the global state MultiFabs state_fc_new_ and state_fc_old_. (Note that these velocities do not have the same time centering as other variables. They are always at the midpoint of the timestep, since they reflect the average face-centered velocity over the timestep.)

This is necessary to implement tracer particles that are advected with the same velocities as the fluid. This requires using the face-centered velocities from the Riemann solver. Ghost cells (ghost faces?) that must be filled with AMR interpolation are required because particles may cross the coarse-fine boundary during a timestep.

Caveat: This causes a slight complication for users, because simulations must now initialize the face-centered velocity field as part of the initial conditions or the code will complain (or possibly crash). If tracer particles are used, boundary conditions for the face-centered velocity must be provided (or the code will crash).

Describe alternatives you've considered The alternative is to use cell-centered velocities for tracer particle advection. However, this is not exactly consistent with the mass flux of the fluid. While this difference may seem small, it can lead to order-of-magnitude errors in tracer particle advection in supersonic turbulence (see https://ui.adsabs.harvard.edu/abs/2013MNRAS.435.1426G/abstract and http://girichidis.de/ringberg2016/Genel.pdf).

Additional context This is implemented in https://github.com/quokka-astro/quokka/pull/490.

[markkrumholz]

Since this impacts users who will not use tracer particles, can we automate this so that users who do not care about face centered quantities can just invoke a single function? For example, can we provide an extrapFCVel function that automatically fill in the face-centered velocities from cell-centered data, which can be invoked by users who don't care about the FC data?

You should also coordinate with @AstroKriel on this, since he also uses FC velocities in the MHD update (which are not exactly the same as the ones you're proposing to use, since they are extrapolated from the cell-centered states, rather than the ones emerging from the Riemann solver).

[BenWibking]

Since this impacts users who will not use tracer particles, can we automate this so that users who do not care about face centered quantities can just invoke a single function? For example, can we provide an extrapFCVel function that automatically fill in the face-centered velocities from cell-centered data, which can be invoked by users who don't care about the FC data?

I'm not sure I understand. What would this function be used for?

You should also coordinate with @AstroKriel on this, since he also uses FC velocities in the MHD update (which are not exactly the same as the ones you're proposing to use, since they are extrapolated from the cell-centered states, rather than the ones emerging from the Riemann solver).

This would affect the indexing of the magnetic field variables, but that shouldn't require any code changes on the MHD side of things, since we already have Physics_Indices<problem_t>::mhdFirstIndex here: https://github.com/quokka-astro/quokka/blob/50ed0d4eb14f71ceda56bb72368addfd01e219b2/src/physics_info.hpp#L43 modified appropriately to allow the velocity index to go first: https://github.com/quokka-astro/quokka/blob/f0e5b481bfcc9b8a67aac84dac707dbd97e387e2/src/physics_info.hpp#L46

Off the top of my head, I can't think of any other way the MHD code might be affected. Maybe I'm missing something...

[markkrumholz]

I realize I may have misinterpreted your proposal: would you have the MF state include face-centered velocities if and only if tracer particles are turned on? In that case, I withdraw my proposal. I was thinking that you were suggesting adding FC velocities to the global state regardless, in which case I wanted to provide a function that would let users who weren't using tracer particles automatically fill in the required FC data, thereby saving them from having to manually set stuff that isn't going to be used. But if the FC variables exist only with tracer particles, then it's fine.

Regarding MHD, I don't think it affects the code directly, but I do think we need to comment the code -- both this new stuff and the MHD stuff, to explain that the FC velocities used by tracer particles and by MHD are not the same, and cannot be used interchangeably. I just want to minimize future confusion on this.

[BenWibking]

Since this impacts users who will not use tracer particles, can we automate this so that users who do not care about face centered quantities can just invoke a single function? For example, can we provide an extrapFCVel function that automatically fill in the face-centered velocities from cell-centered data, which can be invoked by users who don't care about the FC data?

I realize what you meant. The face-centered velocities do not need to have their initial conditions set. That step can be safely removed. The face-centered velocities are not used in the initial conditions, but rather only after the first hydro update has taken place, so their values are always overwritten.

The only reason that users might need to modify their code is to add boundary conditions for the face-centered velocities. But since these are only needed when using tracer particles, I can modify the code so these are only required when tracer particles are turned on.

[markkrumholz]

OK, I think we have converged in that case.

[BenWibking]

I think there are two residual complications:

1. In the current setup, if a user wants to run a simulation with magnetic fields but without tracer particles, boundary conditions must still be provided for the face-centered velocities. This would ideally be avoided, but the code will still work.
2. If a user wants to use both tracer particles and magnetic fields, separate AMR interpolators must be used for different components of the face-centered MultiFabs, since the divergence-free interpolator amrex::FaceDivFree must be used for magnetic fields (but cannot be used for non-divergence-free fields). The MHD code will not work unless we can use different interpolators per component. I will need to check with Weiqun about this.

[BenWibking]

2. If a user wants to use both tracer particles and magnetic fields, separate AMR interpolators must be used for different components of the face-centered MultiFabs, since the divergence-free interpolator amrex::FaceDivFree must be used for magnetic fields (but cannot be used for non-divergence-free fields). The MHD code will not work unless we can use different interpolators per component. I will need to check with Weiqun about this.

For future reference, this problem can be solved by creating an alias MultiFab that references a given range of components within an existing multifab.