In particular, in a test problem with no physics except hydrodynamics, and a maximum refinement level of 0, an initially uniform metallicity field developed features near strong shocks. The changes to the two shock Riemann solver fix this issue by multiplying the color density fields by the same factors by which density itself is corrected, so that metal density and total density are treated identically.
In an AMR version of the same test problem, features again appeared even with the changes in the Riemann solver. These were eliminated with a slight modification to the flux correction code, namely including the metal field in the list of color fields which had to be divided by density before flux correction, then multiplied by density after flux correction.
Here are a number of caveats which maybe should be fixed before this is merged:
I'm not sure this is the correct thing to do physically (maybe there's some drawback I haven't considered)
The density floor near line 138 of euler.F should probably be passed in as a user-specified parameter; not sure if there are other hard-coded density floors elsewhere in the code.
These fixes will almost certainly be non-effective if certain PPM options (e.g. diffusion, steepening) are turned on. It may be straightforward to apply a similar set of fixes that cover these cases.
Other solvers may have similar problems which these exact fixes won't solve.
There are other color fields which will have the issue mentioned above which these changes won't fix, e.g. SNIa, or the optional extra metal fields.
The changes in TestProblemData.h have nothing to do with these issues - they were just necessary to get the AgoraRestart problem type working.
Is the same issue present in the other Riemann solvers?
I'll check. Some version of it was certainly present with HLLC, but I haven't checked since these changes were made, and I haven't checked any of the solvers besides these two.
Alright, it looks like under the these changes HLL, HLLC, and TwoShock (1, 4, and 5) all leave the metallicity featureless (again with the caveats regarding other PPM parameters).
@gbryan and @yipihey should probably take a look. I remember looking at the euler.F options for fraction versus mass conservation, but I can't speak to these particular changes.
So if I understand this correctly, any additional color fields which may be added in the future will also suffer from this problem unless specifically modified in the Riemann solver as you have done here? If this is true, it might be worth mentioning somewhere in the docs so others don't get this poor behavior when adding new fields. Or perhaps this solves the problem for all other color fields and thus doesn't need to be mentioned explicitly in the docs?
Also, great detective work on this, @jforbes ! Thanks for persevering on a very challenging problem and coming up with a solution that will aid us all!
John Wise and I looked at this and think it makes sense. The Flux correction issue is a clear bug and needs to be fixed (it would be good if this was handled more generally and I think John Forbes is working on that but will come as a separate pull request). We also believe the changes to twoshock_flux and inteuler are correct (although at some point it would be good to double-check with some simple tests and cosmological tests). I'm merging.