Implement advanced steps for CMA

[psharda]

273 implements the basic step for consistent multi-fluid advection (CMA) where we conserve the sum of mass scalars to unity. However, more stuff needs to be done during interpolation of advection. Quoting @BenWibking's message from #273 :

There are three extra steps added to the reconstruction of the interface states. One is steepening the reconstruction at composition interfaces, another is flattening the reconstruction under certain conditions, and the last one is modifying the reconstructed partial density states so that they sum to the reconstructed total density. Those will all make it more accurate but require modifying a completely different part of the code, so I think that should go in another PR.

[BenWibking]

We may want to think about this in more detail.

The state-of-the-art method for scalar advection is the Piecewise Parabolic Boltzmann method (https://www.lcse.umn.edu/PPBdocs/PPB-exposition-9-19-06.pdf), which is significantly more complicated but also significantly more accurate than the extra CMA steps ("the accuracy of the PPB advection is roughly equivalent to that of the PPM advection scheme with a 3-fold grid refinement in each spatial dimension")

Update: on further thought, I think we should not implement PPB. However, I think it is worth carefully thinking about whether the contact steepening step in CMA makes sense for multi-D simulations.

[BenWibking]

Mike suggests we should probably not implement these steps. The one exception is that we might consider normalizing the reconstructed partial densities at interfaces, but the other stuff doesn't make sense to do.

[psharda]

Where does this reconstruction occur?

[BenWibking]

Where does this reconstruction occur?

Here: https://github.com/quokka-astro/quokka/blob/a14eabe13d50d025cf033855102ccab2331bf3be/src/hyperbolic_system.hpp#L210

[BenWibking]

@psharda How do the species abundances look in the sim you just ran? If everything is normalized correctly and there isn't anything else weird, can we close this issue?

[psharda]

@BenWibking in the cells where we actually do chemistry (i..e, cells with densities > min_density_allowed), we are conserving the sum of mass fractions to unity. In the cells where we don't do chemistry, the min and max sum of mass fractions is 0.99999999 and 1.00003192, respectively. So, mass fraction is not conserved in these background cells, which only undergo advection but not chemistry. Now, this is ofcourse not a problem for PopIII, but might become an issue for other setups? In any case, CMA is accurate enough for me for now, but we might want to come back to this later on if we encounter issues.

[BenWibking]

Could this be due to how the density floors are implemented?

[psharda]

Could be? The densityFloor is 3 orders of magnitude below the minimum density.

[BenWibking]

Did you enable a separate density floor for the species densities? If so, is there any renormalization of the other species when you do this?

[psharda]

Yes, species have a much lower density floor, because mass fractions of 1e-50 or so can occur. No, I didn't implement a renormalization if the floor is activated - I should do that.

[psharda]

Okay, this implemented this now as a part of #368, so closing this issue.

[psharda]

@BenWibking @markkrumholz I have a question about what we have done for CMA so far, in #273 and #368. The way we implemented CMA was inherited from what we did in FLASH: we conserve the mass fractions (or partial densities) of all the species. We have tested this only on the primordial chemistry network, and it works fine for this network. However, this network does not have any molecules that are not heterogeneous. In metal+ chemistry networks, we will have heterogeneous molecules (e.g., CO). While our scheme may be sufficient to conserve the sum of mass fractions of all species to unity, is it sufficient to conserve the amount of individual elements, like C and O? In such cases, don't we need something more advanced, like the higher order Plewa and Muller formalisms that go beyond a simple renormalization of the mass fractions? Relevant code in quokka/src/hydro_system.hpp here: https://github.com/quokka-astro/quokka/blob/d616a63e7a483785f954682ce7ef233ba71d3464/src/hydro_system.hpp#L1102

[BenWibking]

Reference: Appendix A2 of https://ui.adsabs.harvard.edu/abs/2010MNRAS.404....2G/abstract.

The renormalization step requires the additional step of inverting a matrix, rather than just dividing by the total species mass.

[psharda]

Ah yes, this is what we need. How do you do matrix inversion in amrex?

[BenWibking]

Ah yes, this is what we need. How do you do matrix inversion in amrex?

There is code in Microphysics to do this. Ask Mike ;)

[BenWibking]

Actually, it's here: https://github.com/AMReX-Astro/Microphysics/blob/development/util/linpack.H

Example of how to use it: https://github.com/AMReX-Astro/Microphysics/blob/0a96241b0f29b95347745ad74fd688dfd3d56660/unit_test/test_linear_algebra/test_linear_algebra.H#L130

It solves Ax = b, and overwrites the value of b with the solved-for value of x.

[BenWibking]

You might want to ask if that's the best way to do it, or if there's template magic in Microphysics that can construct the matrix from the network automatically.