Consolidate pressure routines, improve VL+CT & Grackle compatability

[mabruzzo]

This PR builds on changes introduced in PR #151. EDIT: This is now ready for review

Overview

This PR had 2 primary objectives:

1. Refactor the EnzoEOSIdeal::pressure_from_integration (which is used in EnzoMethodMHDVlct) to call the code in EnzoComputePressure. Previously, we basically had parallel implementations of the same code.
2. Removes the restrictions on using EnzoMethodMHDVlct with EnzoMethodGrackle (when primordial_chemistry>1).

Detailed Description

1st Objective

The 1st objective required the following steps:

• introduced the EnzoFieldAdaptor class.

  • The goal of this class is to provide a light-weight wrapper with a single interface for accessing data from instances of the Block and EFltArrayMap classes.
  • The idea is to only initialize an instance of EnzoFieldAdaptor temporarily.
  • For that reason, I took some care to avoid heap allocations within the class's implementation.
  • I think that this class can broadly be useful in other circumstances as well.
  • The implementation of this class is a little quirky to avoid creating problems for one of the solutions for the EnzoEFltArrayMap problem described in issue Unsound const semantics in EnzoEFltArrayMap #172.

• I refactored EnzoComputePressure so that the main calculation is now performed in a static method called EnzoComputePressure::compute_pressure that is compatible with the EnzoFieldAdaptor class (I choose to introduce a new method to avoid breaking the old interfaces).

• To complete the previous step I needed to also refactor EnzoMethodGrackle's setup_grackle_fields, setup_grackle_units, and calculate_pressure methods to use of EnzoFieldAdaptor. (It also made sense to tweak EnzoMethodGrackle's calculate_cooling_time and calculate_temperature methods because they are built on the same machinery as EnzoMethodGrackle::calculate_pressure).

• modified EnzoEOSIdeal::pressure_from_integration to call EnzoComputePressure::compute_pressure.

2nd Objective

To complete the second objective, I mostly just removed some checks from EnzoEOSIdeal that aborted the simulation if it detected primordial_chemistry>1. I also added some documentation to the website (see doc/source/user/problem_method.rst) explaining some of the inconsistencies that are present when using a hydro solver and Grackle with primordial_chemistry>1. To briefly summarize:

• when primordial_chemistry>1 Grackle routines modify the fixed adiabatic index (specified in the Field:gamma - in this scenario, it should be 5/3) based on the abundance of molecular hydrogen and gas temperature.

• Both the Ppm and VL+CT hydro/mhd solvers almost entirely rely the fixed adiabatic index (which overestimates the value used by Grackle).

• The only scenario where the hydro/mhd solvers use the same gamma value as Grackle is to compute pressure in the timestep method (although they still use the fixed adiabatic index in the formula for sound-speed). Note: in reality, we may actually want to use the fixed gamma to compute pressure in this scenario for self consistency within the hydro solver's courant condition - but maybe we should address this separately

There are 2 things to note:

1. The adjustments to the adiabatic index are generally small outside of cells where molecular Hydrogen makes up a large fraction of the gas density.
2. It may be better to move the detailed explanation of the inconsistencies to a different section of the documentation (maybe the developer guide?).

Other Related Changes

Before I refamiliared myself with exactly what the Ppm solver was doing (while writing website documentation), I originally thought the VL+CT solver should compute the pressure values used in reconstruction with the Grackle-suplied routines.

• This is the primary reason why EnzoMethodGrackle::calculate_pressure and the related functions accept stale_depth as an argument. For those who don't recall, this is a frequently used parameter in the VL+CT solver that is used to avoid calculations in the ghost zone that would involve invalid (or "stale") values.

• To maintain consistency with the Ppm solver, I eventually decided not to use the Grackle-supplied routines. However, I choose to keep the support for specifying stale_depth as it took a little work to implement and could be very useful in the future.

I also introduced a simple test called method_grackle_shock_tube.in that sets up a Sod shock tube and models radiative cooling. that is used for checking that enzo-e will not crash when used with the VL+CT solver. This test is and basically sets up a

• This test is not particularly well-motivated (it does not have a predictable result) and should be replaced with a better test in the future (ideally one taken from Enzo). This primarily exists to make sure Enzo-E doesn't crash when run with with Grackle and the VL+CT solver.

• Automating this test would require machinery from PR Grackle Tests #125. To avoid introducing a dependency on another PR and because of the point made in the previous bullet, I choose not to automate this test. If one of the reviewers feels strongly to the contrary, I can automate this test. Otherwise, my intention is open a new Issue documenting that this needs to be addressed after PR Grackle Tests #125 is merged.

• Introducing this test required a minor change to EnzoInitialShockTube::enforce_block in order to initialize the species fields to sensible values.

• While I was making this test, I was inspired to make a small tweak to EnzoMethodGrackle::update_grackle_density_fields to make the grackle_fields parameters an optional argument. In reality it might be better to do that in its own PR. If the reviewers prefer, I can spin that off into its own PR.

Future Direction

For a little while now I've been talking about refactoring how the equation of state stuff is handled in the VL+CT solver. This was first step towards achieving that (it took far more work than I expected). I intend to continue doing that.

As I've mentioned before, I think we can better integrate the machinery behind Field and EnzoEFltArrayMap. When we do that, it may simplify the implementation of EnzoFieldAdaptor (or even remove the need for it).

[mabruzzo]

@peeples Now that PR #151 is merged, this can definitely be reviewed

[gregbryan]

I think it could use a bit more discussion of the EnzoFieldAdaptor class (as you note, eventually it would be great to bring the EnzoEFltArray machinery to ppm and so maybe obviate this, but in the meantime, it is not entirely obvious why it exists and how it works). Otherwise, the code looks very nice to me!

[gregbryan]

Thanks for addressing my suggestions -- this looks good to go to me!