diff --git a/examples/parse_simulation.py b/examples/parse_simulation.py index dcf72ef..e10249a 100644 --- a/examples/parse_simulation.py +++ b/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) ################################################################# diff --git a/pydrad/__init__.py b/pydrad/__init__.py index 2255a3a..002016c 100644 --- a/pydrad/__init__.py +++ b/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: diff --git a/pydrad/configure/configure.py b/pydrad/configure/configure.py index 6f32a44..6846fd2 100644 --- a/pydrad/configure/configure.py +++ b/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']) diff --git a/pydrad/parse/parse.py b/pydrad/parse/parse.py index 09f8f97..2af2eb8 100644 --- a/pydrad/parse/parse.py +++ b/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,29 +344,31 @@ 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: @@ -373,11 +376,11 @@ def _read_trm(self): 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: @@ -385,10 +388,11 @@ def _read_ine(self): 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,27 +476,22 @@ 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 @@ -500,23 +499,31 @@ 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 diff --git a/pydrad/parse/tests/test_strand.py b/pydrad/parse/tests/test_strand.py index 7abfa12..97321ef 100644 --- a/pydrad/parse/tests/test_strand.py +++ b/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): diff --git a/pydrad/parse/util.py b/pydrad/parse/util.py index 978d1aa..23d9780 100644 --- a/pydrad/parse/util.py +++ b/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', diff --git a/pydrad/visualize/animate.py b/pydrad/visualize/animate.py index dc2f5c8..cbd9e2b 100644 --- a/pydrad/visualize/animate.py +++ b/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) diff --git a/pydrad/visualize/plot.py b/pydrad/visualize/plot.py index 51b2a44..117b637 100644 --- a/pydrad/visualize/plot.py +++ b/pydrad/visualize/plot.py @@ -252,7 +252,7 @@ def _plot_profile(profile, axes, **kwargs): ) line1b, = axes[0, 0].plot( profile.coordinate.to(u.Mm), - profile.ion_temperature.to(u.MK), + profile.hydrogen_temperature.to(u.MK), **kwargs, ls='--' ) @@ -264,7 +264,7 @@ def _plot_profile(profile, axes, **kwargs): ) line2b, = axes[0, 1].plot( profile.coordinate.to(u.Mm), - profile.ion_density.to(u.cm**(-3)), + profile.hydrogen_density.to(u.cm**(-3)), **kwargs, ls='--' ) @@ -276,7 +276,7 @@ def _plot_profile(profile, axes, **kwargs): ) line3b, = axes[1, 0].plot( profile.coordinate.to(u.Mm), - profile.ion_pressure.to(u.dyne / (u.cm**2)), + profile.hydrogen_pressure.to(u.dyne / (u.cm**2)), **kwargs, ls='--' )