CurrentDensityMax is updated in Grid_SolveHydroEquation only.
Hence using a RK based solver, the default value was never changed and OutputOnDensity did not work.
This change updates CurrentDensityMax after the RK solvers were executed on each grid.
Hi @Philipp Grete ! Thanks for adding this -- I missed it when I added OutputOnDensity initially. I'm a bit surprised that this works as is, as I think the BaryonField would need a prefix for Grids[grid1] in front of it. Does it in fact work and compile?