respond to code review feedback

examples/parse_simulation.py

@@ -103,7 +103,7 @@
 # grid. Note that each profile is a `~astropy.units.Quantity` and can
 # be easily converted to any compatible unit.
 print(p.electron_temperature)
-print(p.ion_density)
+print(p.hydrogen_density)
 print(p.velocity)
 
 #################################################################

pydrad/__init__.py

@@ -10,12 +10,6 @@
 
 
 def _init_log():
-    """
-    Initializes the SunPy log.
-    In most circumstances this is called automatically when importing
-    SunPy. This code is based on that provided by Astropy see
-    "licenses/ASTROPY.rst".
-    """
     orig_logger_cls = logging.getLoggerClass()
     logging.setLoggerClass(AstropyLogger)
     try:

pydrad/configure/configure.py

@@ -425,8 +425,8 @@ def maximum_cells(self):
         refinement level and :math:`n_{min}` is the minimum allowed number of
         grid cells. Optionally, if the maximum number of cells is specified
         in ``config['grid']['maximum_cells']``, this value will take
-        precedence. In general, it is better to set this value explicitly to a sufficiently
-        large value.
+        precedence. Alternatively, using a value of 30000 will be sufficiently large for
+        any amount of refinement used in HYDRAD.
         """
         if 'maximum_cells' in self.config['grid']:
             return int(self.config['grid']['maximum_cells'])

pydrad/parse/parse.py

@@ -10,6 +10,7 @@
 import numpy as np
 from scipy.interpolate import splev, splrep
 
+import pydrad.util.constants
 from pydrad import log
 from pydrad.configure import Configure
 from pydrad.parse.util import (read_amr_file, read_hstate_file, read_ine_file,
@@ -343,52 +344,55 @@ def __repr__(self):
 
     def _read_amr(self):
         """
-        Parse the adaptive mesh grid file
+        Parse the adaptive mesh refinement (``.amr``) file
         """
         self._amr_data = read_amr_file(self._amr_filename)
 
     def _read_phy(self):
         """
-        Parse the hydrodynamics results file
+        Parse the physical variables (``.phy``) file and set attributes
+        for relevant quantities.
         """
         if not self._phy_filename.is_file():
-            log.warning(f'{self._phy_filename} not found. Skipping parsing of phy files. Set read_phy=False to suppress this warning.')
+            log.warning(f'{self._phy_filename} not found. Skipping parsing of .phy files. Set read_phy=False to suppress this warning.')
             return
         self._phy_data = read_phy_file(self._phy_filename)
         # NOTE: only adding three columns in this manner as the remaining columns are dealt with explicitly
         # below because of the overlap with the derived quantities from the .amr files.
-        for column in ['sound_speed', 'electron_heat_flux', 'ion_heat_flux']:
+        for column in ['sound_speed', 'electron_heat_flux', 'hydrogen_heat_flux']:
             _add_property(column, '_phy_data')
 
     def _read_trm(self):
         """
-        Parse the equation terms file
+        Parse the equation terms (``.trm``) file and set attributes
+        for relevant quantities.
         """
         if not self._trm_filename.is_file():
-            log.warning(f'{self._trm_filename} not found. Skipping parsing of trm files. Set read_trm=False to suppress this warning.')
+            log.warning(f'{self._trm_filename} not found. Skipping parsing of .trm files. Set read_trm=False to suppress this warning.')
             return
         self._trm_data = read_trm_file(self._trm_filename)
         for col in self._trm_data.colnames:
             _add_property(col, '_trm_data')
 
     def _read_ine(self):
         """
-        Parse non-equilibrium ionization population fraction files
+        Parse non-equilibrium ionization population fraction (``.ine``) file
         and set attributes for relevant quantities
         """
         if not self._ine_filename.is_file():
-            log.warning(f'{self._ine_filename} not found. Skipping parsing of ine files. Set read_ine=False to suppress this warning.')
+            log.warning(f'{self._ine_filename} not found. Skipping parsing of .ine files. Set read_ine=False to suppress this warning.')
             return
         self._ine_data = read_ine_file(self._ine_filename, self.coordinate.shape[0])
         for col in self._ine_data.colnames:
             _add_property(col, '_ine_data')
 
     def _read_hstate(self):
         """
-        Parse the hydrogen energy level populations file
+        Parse the hydrogen energy level populations (``.hstate``) file and set
+        the relevant attributes.
         """
         if not self._hstate_filename.is_file():
-            log.warning(f'{self._hstate_filename} not found. Skipping parsing of hstate files. Set read_hstate=False to suppress this warning.')
+            log.warning(f'{self._hstate_filename} not found. Skipping parsing of .hstate files. Set read_hstate=False to suppress this warning.')
             return
         self._hstate_data = read_hstate_file(self._hstate_filename)
         for col in self._hstate_data.colnames:
@@ -440,7 +444,7 @@ def grid_edges_right(self) -> u.cm:
     @u.quantity_input
     def coordinate(self) -> u.cm:
         """
-        Coordinate in the field-aligned direction
+        Spatial coordinate in the field-aligned direction
         """
         return self.grid_centers
 
@@ -453,12 +457,12 @@ def electron_mass_density(self) -> u.Unit('g cm-3'):
         # 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.ion_mass_density
+        return self.mass_density
 
     @property
     @u.quantity_input
-    def ion_mass_density(self) -> u.Unit('g cm-3'):
-        return self._amr_data['ion_mass_density']
+    def mass_density(self) -> u.Unit('g cm-3'):
+        return self._amr_data['mass_density']
 
     @property
     @u.quantity_input
@@ -472,51 +476,54 @@ def electron_energy_density(self) -> u.Unit('erg cm-3'):
 
     @property
     @u.quantity_input
-    def ion_energy_density(self) -> u.Unit('erg cm-3'):
-        return self._amr_data['ion_energy_density']
+    def hydrogen_energy_density(self) -> u.Unit('erg cm-3'):
+        return self._amr_data['hydrogen_energy_density']
 
     @property
     @u.quantity_input
     def velocity(self) -> u.cm/u.s:
         if hasattr(self, '_phy_data'):
             return self._phy_data['velocity']
-        return self.momentum_density / self.ion_mass_density
+        return self.momentum_density / self.mass_density
 
     @property
     @u.quantity_input
-    def ion_density(self) -> u.Unit('cm-3'):
+    def hydrogen_density(self) -> u.Unit('cm-3'):
         if hasattr(self, '_phy_data'):
-            return self._phy_data['ion_density']
-        # NOTE: This is pulled directly from the HYDRAD source code
-        # https://github.com/rice-solar-physics/HYDRAD/blob/master/Resources/source/constants.h
-        # and is the average ion mass for a H-He plasma computed by weighting the H and He ion
-        # masses with a particular set of abundances.
-        m_ion = 2.171e-24 * u.g
-        return self.ion_mass_density / m_ion
+            return self._phy_data['hydrogen_density']
+        return self.mass_density / pydrad.util.constants.m_avg_ion
 
     @property
     @u.quantity_input
     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
-        return self.ion_density
+        # choose that. Otherwise, the electron and ion densities are assumed
+        # to be the same.
+        return self.hydrogen_density
 
     @property
     @u.quantity_input
     def electron_pressure(self) -> u.Unit('dyne cm-2'):
         if hasattr(self, '_phy_data'):
             return self._phy_data['electron_pressure']
-        gamma = 5/3
-        return self.electron_energy_density / (gamma - 1)
+        return self.electron_energy_density / (pydrad.util.constants.gamma - 1)
 
     @property
     @u.quantity_input
-    def ion_pressure(self) -> u.Unit('dyne cm-2'):
+    def hydrogen_pressure(self) -> u.Unit('dyne cm-2'):
         if hasattr(self, '_phy_data'):
-            return self._phy_data['ion_pressure']
-        gamma = 5/3
-        return (self.ion_energy_density - self.momentum_density**2/(2*self.ion_mass_density))/(gamma - 1)
+            return self._phy_data['hydrogen_pressure']
+        return (
+            self.hydrogen_energy_density
+            - self.momentum_density**2/(2*self.mass_density)
+        )/(pydrad.util.constants.gamma - 1)
+
+    @property
+    @u.quantity_input
+    def total_pressure(self) -> u.Unit('dyne cm-2'):
+        return self.electron_pressure + self.hydrogen_pressure
 
     @property
     @u.quantity_input
@@ -527,10 +534,10 @@ def electron_temperature(self) -> u.Unit('K'):
 
     @property
     @u.quantity_input
-    def ion_temperature(self) -> u.Unit('K'):
+    def hydrogen_temperature(self) -> u.Unit('K'):
         if hasattr(self, '_phy_data'):
-            return self._phy_data['ion_temperature']
-        return self.ion_pressure / (const.k_B*self.ion_density)
+            return self._phy_data['hydrogen_temperature']
+        return self.hydrogen_pressure / (const.k_B*self.hydrogen_density)
 
     def spatial_average(self, quantity, bounds=None):
         """
@@ -577,7 +584,7 @@ def column_emission_measure(self, bins: u.K = None, bounds: u.cm = None):
             bins = 10.0**(np.arange(3.0, 8.0, 0.05)) * u.K
         if bounds is None:
             bounds = self.grid_edges[[0, -1]]
-        weights = self.electron_density * self.ion_density * self.grid_widths
+        weights = self.electron_density * self.hydrogen_density * self.grid_widths
         H, _, _ = np.histogram2d(self.grid_centers, self.electron_temperature,
                                  bins=(bounds, bins), weights=weights)
         return H.squeeze(), bins

pydrad/parse/tests/test_strand.py

@@ -14,15 +14,15 @@
     'grid_widths',
     'grid_centers',
     'electron_temperature',
-    'ion_temperature',
+    'hydrogen_temperature',
     'electron_density',
-    'ion_density',
+    'hydrogen_density',
     'electron_pressure',
-    'ion_pressure',
+    'hydrogen_pressure',
     'velocity',
     'sound_speed',
     'electron_heat_flux',
-    'ion_heat_flux',
+    'hydrogen_heat_flux',
     'mass_drhobydt',
     'mass_advection',
     'momentum_drho_vbydt',
@@ -77,11 +77,11 @@ def test_has_quantity(strand, quantity):
 
 @pytest.mark.parametrize('quantity', [
     'electron_temperature',
-    'ion_temperature',
+    'hydrogen_temperature',
     'electron_density',
-    'ion_density',
+    'hydrogen_density',
     'electron_pressure',
-    'ion_pressure',
+    'hydrogen_pressure',
     'velocity',
 ])
 def test_no_amr_run_has_quantity(strand_only_amr, quantity):

pydrad/parse/util.py

@@ -19,8 +19,11 @@
 
 def read_amr_file(filename):
     """
-    Parse ``.amr`` files containing grid parameters and conserved quantities as
-    a function of position
+    Parse the adaptive mesh refinement ``.amr`` files containing grid parameters
+    and conserved quantities as a function of position.
+
+    .. 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
@@ -29,18 +32,18 @@ def read_amr_file(filename):
     columns = [
         'grid_centers',
         'grid_widths',
-        'ion_mass_density',
+        'mass_density',
         'momentum_density',
         'electron_energy_density',
-        'ion_energy_density',
+        'hydrogen_energy_density',
     ]
     units = {
         'grid_centers': 'cm',
         'grid_widths': 'cm',
-        'ion_mass_density': 'g cm-3',
+        'mass_density': 'g cm-3',
         'momentum_density': 'g cm-2 s-1',
         'electron_energy_density': 'erg cm-3',
-        'ion_energy_density': 'erg cm-3',
+        'hydrogen_energy_density': 'erg cm-3',
     }
     table = astropy.table.QTable.read(
         filename,
@@ -61,33 +64,37 @@ def read_amr_file(filename):
 
 def read_phy_file(filename):
     """
-    Parse ``.phy`` files containing thermodynamic quantities as a function of position.
+    Parse physical variables ``.phy`` files containing thermodynamic
+    quantities as a function of position.
+
+    .. note:: This is not intended for direct use. Use the `pydrad.parse.Strand`
+              object instead.
     """
     columns = [
         'coordinate',
         'velocity',
         'sound_speed',
         'electron_density',
-        'ion_density',
+        'hydrogen_density',
         'electron_pressure',
-        'ion_pressure',
+        'hydrogen_pressure',
         'electron_temperature',
-        'ion_temperature',
+        'hydrogen_temperature',
         'electron_heat_flux',
-        'ion_heat_flux',
+        'hydrogen_heat_flux',
     ]
     units = {
         'coordinate': 'cm',
         'velocity': 'cm / s',
         'sound_speed': 'cm / s',
         'electron_density': 'cm-3',
-        'ion_density': 'cm-3',
+        'hydrogen_density': 'cm-3',
         'electron_pressure': 'dyne cm-2',
-        'ion_pressure': 'dyne cm-2',
+        'hydrogen_pressure': 'dyne cm-2',
         'electron_temperature': 'K',
-        'ion_temperature': 'K',
+        'hydrogen_temperature': 'K',
         'electron_heat_flux': 'erg s-1 cm-2',
-        'ion_heat_flux': 'erg s-1 cm-2',
+        'hydrogen_heat_flux': 'erg s-1 cm-2',
     }
     return astropy.table.QTable.read(
         filename,
@@ -99,7 +106,8 @@ def read_phy_file(filename):
 
 def read_ine_file(filename, n_s):
     """
-    Parse ``.ine`` files containing ionization fraction as a function of position
+    Parse the ionization non-equilibrium ``.ine`` files containing
+    ionization fraction as a function of position.
 
     .. note:: This is not intended for direct use. Use the `pydrad.parse.Strand`
               object instead.
@@ -220,7 +228,7 @@ def read_hstate_file(filename):
     """
     Parse the ``.hstate`` files containing the hydrogen level populations as a function of position
     """
-    columns = [f'level_population_hydrogen_{i}' for i in range(1,7)]
+    columns = [f'hydrogen_I_level_{i}' for i in range(1,6)] + ['hydrogen_II_fraction']
     return astropy.table.QTable.read(
         filename,
         format='ascii',

pydrad/visualize/animate.py

@@ -42,11 +42,11 @@ def animate_strand(strand, **kwargs):
     def update_plot(i):
         p = strand[i]
         l1a.set_data(p.coordinate.to(u.Mm), p.electron_temperature.to(u.MK))
-        l1b.set_data(p.coordinate.to(u.Mm), p.ion_temperature.to(u.MK))
+        l1b.set_data(p.coordinate.to(u.Mm), p.hydrogen_temperature.to(u.MK))
         l2a.set_data(p.coordinate.to(u.Mm), p.electron_density.to(u.cm**(-3)))
-        l2b.set_data(p.coordinate.to(u.Mm), p.ion_density.to(u.cm**(-3)))
+        l2b.set_data(p.coordinate.to(u.Mm), p.hydrogen_density.to(u.cm**(-3)))
         l3a.set_data(p.coordinate.to(u.Mm), p.electron_pressure.to(u.dyne/(u.cm**2)))
-        l3b.set_data(p.coordinate.to(u.Mm), p.ion_pressure.to(u.dyne/(u.cm**2)))
+        l3b.set_data(p.coordinate.to(u.Mm), p.hydrogen_pressure.to(u.dyne/(u.cm**2)))
         l4.set_data(p.coordinate.to(u.Mm), p.velocity.to(u.km/u.s))
         fig.suptitle(r'$t={:.0f}$ {}'.format(
             strand.time[i].value, strand.time[i].unit), y=0.905)