Account for possibility of electron mass density in .amr files (#189)

* parse electron mass density; fix reading of .Hstate

* update return logic for electron density and mass density

* add a few basic tests

* add Hstate files to tests, and test them

* remove Hstate tests; fix bug in rho_e

* add clarifying comments

* pre-commit fixes

Author: jwreep

examples/configure_beam_simulation.py

@@ -52,7 +52,7 @@
 # For consistency with VAL C, we also need to set a few boundary conditions:
 config['general']['footpoint_height'] = 2.26 * u.Mm
 config['initial_conditions']['footpoint_temperature'] = 24000 * u.K
-config['initial_conditions']['footpoint_density'] = 4.2489e9 * u.cm**(-3)
+config['initial_conditions']['footpoint_density'] = 4.2486e9 * u.cm**(-3)
 config['solver']['minimum_radiation_temperature'] = 24000 * u.K
 config['solver']['minimum_temperature'] = 4170 * u.K
 
@@ -62,8 +62,10 @@
 # solve an approximation to radiative transfer for hydrogen with the caveat that
 # it will slow the code by well over an order of magnitude.  It is most
 # useful for users who wish to synthesize chromospheric line profiles.
-# Change the value to True if you would like to try using it.
+# Change the value to True if you would like to try using it, and add
+# the minimum_density_limit parameter.
 config['radiation']['nlte_chromosphere'] = False
+config['radiation']['minimum_density_limit'] = 4.2486e9*u.cm**(-3)
 
 #################################################################
 # This takes care of the chromosphere, but we should also set up the radiation

pydrad/parse/parse.py

@@ -470,13 +470,11 @@ def coordinate(self) -> u.cm:
     @property
     @u.quantity_input
     def electron_mass_density(self) -> u.Unit('g cm-3'):
-        # TODO: Account for possible presence of both electron
-        # and ion mass densities. If the electron mass density
-        # is present in this file, it will shift all of the
-        # indices in the .amr file by one.
         if 'electron_mass_density' in self._amr_data.colnames:
             return self._amr_data['electron_mass_density']
-        return self.mass_density
+        if hasattr(self, '_phy_data'):
+            return self._phy_data['electron_density'] * const.m_e
+        return (self.mass_density / pydrad.util.constants.m_avg_ion) * const.m_e
 
     @property
     @u.quantity_input
@@ -548,9 +546,8 @@ def electron_density(self) -> u.Unit('cm-3'):
         """
         if hasattr(self, '_phy_data'):
             return self._phy_data['electron_density']
-        # FIXME: If this exists as a separate column in the .amr file then
-        # choose that. Otherwise, the electron and ion densities are assumed
-        # to be the same.
+        if 'electron_mass_density' in self._amr_data.colnames:
+            return self._amr_data['electron_mass_density'] / const.m_e
         return self.hydrogen_density
 
     @property

pydrad/parse/tests/test_strand.py

@@ -13,6 +13,8 @@
     'grid_edges',
     'grid_widths',
     'grid_centers',
+    'electron_mass_density',
+    'mass_density',
     'electron_temperature',
     'hydrogen_temperature',
     'electron_density',
@@ -91,6 +93,7 @@ def test_has_quantity(strand, quantity):
     'electron_pressure',
     'hydrogen_pressure',
     'velocity',
+    'electron_mass_density',
 ])
 def test_no_amr_run_has_quantity(strand_only_amr, quantity):
     # check that Strand falls back to deriving thermodynamic quantities from
@@ -159,3 +162,9 @@ def test_scale_file_units(strand):
     assert strand[0].advective_timescale.unit == u.Unit('s')
     assert strand[0].electron_conductive_timescale.unit == u.Unit('s')
     assert strand[0].collisional_timescale.unit == u.Unit('s')
+
+def test_amr_file_units(strand, strand_only_amr):
+    assert strand[0].mass_density.unit == u.Unit('g cm-3')
+    assert strand_only_amr[0].mass_density.unit == u.Unit('g cm-3')
+    assert strand[0].electron_mass_density.unit == u.Unit('g cm-3')
+    assert strand_only_amr[0].electron_mass_density.unit == u.Unit('g cm-3')

pydrad/parse/util.py

@@ -26,13 +26,10 @@ def read_amr_file(filename):
     .. note:: This is not intended for direct use. Use the `pydrad.parse.Strand`
               object instead.
     """
-    # TODO: Account for possible presence of both electron
-    # and ion mass densities. If the electron mass density
-    # is present in this file, it will shift all of the
-    # columns from the index=2 position to the right.
     columns = [
         'grid_centers',
         'grid_widths',
+        'electron_mass_density',
         'mass_density',
         'momentum_density',
         'electron_energy_density',
@@ -41,6 +38,7 @@ def read_amr_file(filename):
     units = {
         'grid_centers': 'cm',
         'grid_widths': 'cm',
+        'electron_mass_density': 'g cm-3',
         'mass_density': 'g cm-3',
         'momentum_density': 'g cm-2 s-1',
         'electron_energy_density': 'erg cm-3',
@@ -54,6 +52,15 @@ def read_amr_file(filename):
     # NOTE: This is done after creating the table because the
     # remaining number of columns can be variable and thus we
     # cannot assign all of the column names at once.
+
+    # The columns we care about are doubles in HYDRAD, while the
+    # other columns are integers with information about the
+    # refinement level of the grid cell.  As a result, if electron
+    # mass density is not present in the .amr file, then the
+    # seventh column is an integer.
+    if table.dtype[len(columns)-1] == np.int64:
+        columns.remove('electron_mass_density')
+        del units['electron_mass_density']
     table.rename_columns(
         table.colnames[:len(columns)],
         columns,
@@ -62,7 +69,6 @@ def read_amr_file(filename):
         table[column].unit = units[column]
     return table
 
-
 def read_phy_file(filename):
     """
     Parse physical variables ``.phy`` files containing thermodynamic
@@ -229,12 +235,17 @@ def read_hstate_file(filename):
     """
     Parse the ``.hstate`` files containing the hydrogen level populations as a function of position
     """
-    columns = [f'hydrogen_I_level_{i}' for i in range(1,6)] + ['hydrogen_II_fraction']
-    return astropy.table.QTable.read(
+    columns = ['coordinate'] + [f'hydrogen_I_level_{i}' for i in range(1,6)] + ['hydrogen_II_fraction']
+    table = astropy.table.QTable.read(
         filename,
-        format='ascii',
+        format='ascii.no_header',
+        delimiter='\t',
         names=columns,
     )
+    # The coordinate is already stored from the .amr file, so we don't need to save it.
+    # However, we need to parse the correct number of columns.
+    table.remove_column('coordinate')
+    return table
 
 def read_scl_file(filename):
     """