diff --git a/.gitignore b/.gitignore index 2413634..f85386d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ site/ docs/api/* docs/config_tables.rst docs/generated +docs/sg_execution_times.rst # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/examples/configure_simulation.py b/examples/configure_simulation.py index b88b3f5..3d298d1 100644 --- a/examples/configure_simulation.py +++ b/examples/configure_simulation.py @@ -27,7 +27,7 @@ # Mm loop lasting 5000 s heated by a single 200 s nanoflare # solved on an adaptive grid. # A complete list of configuration parameters can be found in -# the `configuration-tables`_ page. +# the :ref:`configuration-tables` page. config_dict = { 'general': { 'loop_length': 80*u.Mm, diff --git a/examples/parse_simulation.py b/examples/parse_simulation.py index 45ae81a..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) ################################################################# @@ -125,7 +125,7 @@ # single time step as a function of field-aligned coordinate. # The `~pydrad.parse.Profile` object provides a simple quicklook # method for this, -s[0].peek() +p.peek() ################################################################# # Similarly, we can call this method on a `~pydrad.parse.Strand` to 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 186a765..6846fd2 100644 --- a/pydrad/configure/configure.py +++ b/pydrad/configure/configure.py @@ -352,10 +352,13 @@ def poly_fit_magnetic_field(self): Sixth-order polynomial fit coefficients for computing flux tube expansion """ - fit = self._fit_poly_domains('poly_fit_magnetic_field', 'G') + fit_results = self._fit_poly_domains( + **self.config['general']['poly_fit_magnetic_field'], + target_unit='G', + ) return self.env.get_template('coefficients.cfg').render( date=self.date, - **fit, + **fit_results, ) @property @@ -364,30 +367,37 @@ def poly_fit_gravity(self): Sixth-order polynomial fit coefficients for computing gravitational acceleration """ - fit = self._fit_poly_domains('poly_fit_gravity', 'cm s-2') + fit_results = self._fit_poly_domains( + **self.config['general']['poly_fit_gravity'], + target_unit='cm s-2' + ) return self.env.get_template('coefficients.cfg').render( date=self.date, - **fit, + **fit_results, ) - def _fit_poly_domains(self, name, unit): + def _fit_poly_domains(self, x, y, domains, order, target_unit): """ Perform polynomial fit to quantity as a function of field aligned coordinate over multiple domains and return fitting coefficients. """ - # TODO: refactor to be independent of dictionary - fit = copy.deepcopy(self.config['general'][name]) - x = (fit['x'] / self.config['general']['loop_length']).decompose().to(u.dimensionless_unscaled).value - y = fit['y'].to(unit).value + x /= self.config['general']['loop_length'] + x = x.decompose().to_value(u.dimensionless_unscaled) + y = y.to_value(target_unit) coefficients = [] minmax = [] - for i in range(len(fit['domains'])-1): - i_d = np.where(np.logical_and(x>=fit['domains'][i], x<=fit['domains'][i+1])) - coefficients.append(np.polyfit(x[i_d], y[i_d], fit['order'])[::-1]) + for i in range(len(domains)-1): + i_d = np.where(np.logical_and(x>=domains[i], x<=domains[i+1])) + # NOTE: The order is reversed because HYDRAD expects the opposite order of + # what NumPy outputs + coefficients.append(np.polyfit(x[i_d], y[i_d], order)[::-1]) minmax.append([y[i_d].min(), y[i_d].max()]) - fit['minmax'] = minmax - fit['coefficients'] = coefficients - return fit + return { + 'domains': domains, + 'order': order, + 'minmax': minmax, + 'coefficients': coefficients, + } @property def minimum_cells(self): @@ -415,7 +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. + 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 3ee8ec1..5458553 100644 --- a/pydrad/parse/parse.py +++ b/pydrad/parse/parse.py @@ -1,20 +1,20 @@ """ Interface for easily parsing HYDRAD results """ -import glob -import os import pathlib +import astropy.constants as const import astropy.units as u import h5py import matplotlib.pyplot as plt import numpy as np -import plasmapy.particles -from pandas import read_csv 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, + read_phy_file, read_trm_file) from pydrad.visualize import (animate_strand, plot_histogram, plot_profile, plot_strand, plot_time_distance, plot_time_mesh) @@ -172,10 +172,15 @@ def initial_conditions(self): `~pydrad.parse.Profile` for the solutions to the hydrostatic equations used as initial conditions. """ + profile_kwargs = self._profile_kwargs.copy() + # NOTE: These files do not exist for the initial conditions + profile_kwargs['read_hstate'] = False + profile_kwargs['read_ine'] = False + profile_kwargs['read_trm'] = False return InitialProfile(self.hydrad_root, 0*u.s, master_time=[0]*u.s, - **self._profile_kwargs) + **profile_kwargs) def peek(self, **kwargs): """ @@ -216,11 +221,6 @@ def get_uniform_grid(self, delta_s: u.cm) -> u.cm: return np.arange( 0, self.loop_length.to(u.cm).value, delta_s.to(u.cm).value)*u.cm - @u.quantity_input - def get_unique_grid(self) -> u.cm: - all_coordinates = [p.coordinate.to(u.cm).value for p in self] - return np.unique(np.concatenate(all_coordinates).ravel()) * u.cm - @u.quantity_input def to_constant_grid(self, name, grid: u.cm, order=1): """ @@ -280,7 +280,7 @@ class Profile(object): Parameters ---------- - hydrad_root : `str` + hydrad_root : path-like Path to HYDRAD directory time : `~astropy.units.Quantity` Timestep corresponding to the profile of interest @@ -297,12 +297,13 @@ def __init__(self, hydrad_root, time: u.s, **kwargs): self._master_time = get_master_time(self.hydrad_root, read_from_cfg=kwargs.get('read_from_cfg', False)) # Read results files - self._read_phy() self._read_amr() - if kwargs.get('read_hstate', True): - self._read_hstate() + if kwargs.get('read_phy', True): + self._read_phy() if kwargs.get('read_ine', True): self._read_ine() + if kwargs.get('read_hstate', True): + self._read_hstate() if kwargs.get('read_trm', True): self._read_trm() @@ -336,152 +337,61 @@ def __repr__(self): Filename: {self._phy_filename} Timestep #: {self._index}""" - def _read_phy(self): + def _read_amr(self): """ - Parse the hydrodynamics results file + Parse the adaptive mesh refinement (``.amr``) file """ - self._phy_data = np.loadtxt(self._phy_filename) + self._amr_data = read_amr_file(self._amr_filename) - def _read_amr(self): + def _read_phy(self): """ - Parse the adaptive mesh grid file + Parse the physical variables (``.phy``) file and set attributes + for relevant quantities. """ - # TODO: Read other quantities from .amr file - with self._amr_filename.open() as f: - lines = f.readlines() - self._grid_centers = np.zeros((int(lines[3]),)) - self._grid_widths = np.zeros((int(lines[3]),)) - for i, l in enumerate(lines[4:]): - tmp = np.array(l.split(), dtype=float) - self._grid_centers[i] = tmp[0] - self._grid_widths[i] = tmp[1] + 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.') + 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', 'hydrogen_heat_flux']: + _add_property(column, '_phy_data') def _read_trm(self): """ - Parse the equation terms file - - The files come in sets of 5 rows with variable number of columns: - -- Loop coordinate (1 column), and at each position: - -- Terms of mass equation (2 columns) - -- Terms of momentum equation (6 columns) - -- Terms of electron energy equation (11 columns) - -- Terms of hydrogen energy equation (11 columns) - """ - mass_columns = ['mass_drhobydt','mass_advection'] - momentum_columns = ['momentum_drho_vbydt','momentum_advection','momentum_pressure_gradient', - 'momentum_gravity','momentum_viscous_stress','momentum_numerical_viscosity'] - electron_columns = ['electron_dTEKEbydt','electron_enthalpy','electron_conduction', - 'electron_gravity','electron_collisions','electron_heating','electron_radiative_loss', - 'electron_electric_field','electron_viscous_stress','electron_numerical_viscosity', - 'electron_ionization'] - hydrogen_columns = ['hydrogen_dTEKEbydt','hydrogen_enthalpy','hydrogen_conduction', - 'hydrogen_gravity','hydrogen_collisions','hydrogen_heating','hydrogen_radiative_loss', - 'hydrogen_electric_field','hydrogen_viscous_stress','hydrogen_numerical_viscosity', - 'hydrogen_ionization'] - - n_mass = len(mass_columns) - n_momentum = len(momentum_columns) - n_electron = len(electron_columns) - n_hydrogen = len(hydrogen_columns) - - # Terms from the mass equation: - mass_terms = read_csv(self._trm_filename, sep='\t', header=None, names=mass_columns, - skiprows=lambda x: x % 5 != 1 ) - # Terms from the momentum equation: - momentum_terms = read_csv(self._trm_filename, sep='\t', header=None, names=momentum_columns, - skiprows=lambda x: x % 5 != 2 ) - # Terms from the electron energy equation: - electron_terms = read_csv(self._trm_filename, sep='\t', header=None, names=electron_columns, - skiprows=lambda x: x % 5 != 3 ) - - # Terms from the hydrogen energy equation: - hydrogen_terms = read_csv(self._trm_filename, sep='\t', header=None, names=hydrogen_columns, - skiprows=lambda x: x % 5 != 4 ) - - offsets = [n_mass, - n_mass+n_momentum, - n_mass+n_momentum+n_electron - ] - - n_elements = len(mass_terms) - self._trm_data = np.zeros([n_elements, offsets[2]+n_hydrogen]) - - for i in range(n_elements): - for j in range(n_mass): - self._trm_data[i,j] = mass_terms[mass_columns[j]][i] - for j in range(n_momentum): - self._trm_data[i,j+offsets[0]] = momentum_terms[momentum_columns[j]][i] - for j in range(n_electron): - self._trm_data[i,j+offsets[1]] = electron_terms[electron_columns[j]][i] - for j in range(n_hydrogen): - self._trm_data[i,j+offsets[2]] = hydrogen_terms[hydrogen_columns[j]][i] - - properties = [] - for i in range(n_mass): - properties += [(mass_columns[i], '_trm_data', i, 'g cm-3 s-1')] - for i in range(n_momentum): - properties += [(momentum_columns[i], '_trm_data', i+offsets[0], 'dyne cm-3 s-1')] - for i in range(n_electron): - properties += [(electron_columns[i], '_trm_data', i+offsets[1], 'erg cm-3 s-1')] - for i in range(n_hydrogen): - properties += [(hydrogen_columns[i], '_trm_data', i+offsets[2], 'erg cm-3 s-1')] - - for p in properties: - add_property(*p) + 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.') + 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 """ - # TODO: clean this up somehow? I've purposefully included - # a lot of comments because the format of this file makes - # the parsing code quite opaque - try: - with self._ine_filename.open() as f: - lines = f.readlines() - except FileNotFoundError: - log.debug(f'{self._ine_filename} not found') + 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.') return - # First parse all of the population fraction arrays - n_s = self.coordinate.shape[0] - # NOTE: Have to calculate the number of elements we have - # computed population fractions for as we do not necessarily - # know this ahead of time - n_e = int(len(lines)/n_s - 1) - # The file is arranged in n_s groups of n_e+1 lines each where the first - # line is the coordinate and the subsequent lines are the population fraction - # for each element, with each column corresponding to an ion of that element - # First, separate by coordinate - pop_lists = [lines[i*(n_e+1)+1:(i+1)*(n_e+1)] for i in range(n_s)] - # Convert each row of each group into a floating point array - pop_lists = [[np.array(l.split(), dtype=float) for l in p] for p in pop_lists] - # NOTE: each row has Z+2 entries as the first entry is the atomic number Z - # Get these from just the first group as the number of elements is the same - # for each - Z = np.array([p[0] for p in pop_lists[0]], dtype=int) - pop_arrays = [np.zeros((n_s, z+1))for z in Z] - for i, p in enumerate(pop_lists): - for j, l in enumerate(p): - pop_arrays[j][i, :] = l[1:] # Skip first entry, it is the atomic number - - # Then set attributes for each ion of each element - for z, p in zip(Z, pop_arrays): - name = plasmapy.particles.element_name(z) - attr = f'_population_fraction_{name}' - setattr(self, attr, p) - for p in [(f'{attr[1:]}_{i+1}', attr, i, '') for i in range(z+1)]: - add_property(*p) + 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. """ - try: - self._hstate_data = np.loadtxt(self._hstate_filename) - except OSError: - log.debug(f'{self._hstate_filename} not found') - self._hstate_data = None + 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.') + return + self._hstate_data = read_hstate_file(self._hstate_filename) + for col in self._hstate_data.colnames: + _add_property(col, '_hstate_data') @property @u.quantity_input @@ -489,7 +399,7 @@ def grid_centers(self) -> u.cm: """ Spatial location of the grid centers """ - return u.Quantity(self._grid_centers, 'cm') + return self._amr_data['grid_centers'] @property @u.quantity_input @@ -497,7 +407,7 @@ def grid_widths(self) -> u.cm: """ Spatial width of each grid cell """ - return u.Quantity(self._grid_widths, 'cm') + return self._amr_data['grid_widths'] @property @u.quantity_input @@ -525,6 +435,153 @@ def grid_edges_right(self) -> u.cm: """ return self.grid_edges[1:] + @property + @u.quantity_input + def coordinate(self) -> u.cm: + """ + Spatial coordinate, :math:`s`, in the field-aligned direction. + + An alias for `~pydrad.parse.Profile.grid_centers`. + """ + return self.grid_centers + + @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 + + @property + @u.quantity_input + def mass_density(self) -> u.Unit('g cm-3'): + r""" + Mass density, :math:`\rho`, as a function of :math:`s`. + + .. note:: This is a conserved quantity in HYDRAD. + """ + return self._amr_data['mass_density'] + + @property + @u.quantity_input + def momentum_density(self) -> u.Unit('g cm-2 s-1'): + r""" + Momentum density, :math:`\rho v`, as a function of :math:`s`. + + .. note:: This is a conserved quantity in HYDRAD. + """ + return self._amr_data['momentum_density'] + + @property + @u.quantity_input + def electron_energy_density(self) -> u.Unit('erg cm-3'): + r""" + Electron energy density, :math:`E_e`, as a function of :math:`s`. + + .. note:: This is a conserved quantity in HYDRAD. + """ + return self._amr_data['electron_energy_density'] + + @property + @u.quantity_input + def hydrogen_energy_density(self) -> u.Unit('erg cm-3'): + r""" + Hydrogen energy density, :math:`E_H`, as a function of :math:`s`. + + .. note:: This is a conserved quantity in HYDRAD. + """ + return self._amr_data['hydrogen_energy_density'] + + @property + @u.quantity_input + def velocity(self) -> u.cm/u.s: + r""" + Bulk velocity, :math:`v`, as a function of :math:`s`. + """ + if hasattr(self, '_phy_data'): + return self._phy_data['velocity'] + return self.momentum_density / self.mass_density + + @property + @u.quantity_input + def hydrogen_density(self) -> u.Unit('cm-3'): + r""" + Hydrogen density, :math:`n_H=\rho/\bar{m_i}`, as a function of :math:`s`, + where :math:`\bar{m_i} is the average ion mass of a H-He plasma. + """ + if hasattr(self, '_phy_data'): + 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'): + r""" + Electron density, :math:`n_e`, as a function of :math:`s`. + In nearly all cases, this is equal to `~pydrad.parse.Profile.hydrogen_density`. + """ + 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. + return self.hydrogen_density + + @property + @u.quantity_input + def electron_pressure(self) -> u.Unit('dyne cm-2'): + r""" + Electron pressure, :math:`P_e=(\gamma - 1)E_e`, as a function of :math:`s`. + """ + if hasattr(self, '_phy_data'): + return self._phy_data['electron_pressure'] + return self.electron_energy_density / (pydrad.util.constants.gamma - 1) + + @property + @u.quantity_input + def hydrogen_pressure(self) -> u.Unit('dyne cm-2'): + r""" + Hydrogen pressure, :math:`P_H = (\gamma - 1)(E_H - \rho v^2/2)`, as a function of :math:`s`. + """ + if hasattr(self, '_phy_data'): + 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'): + r""" + Total pressure, :math:`P = P_e + P_H`, as a function of :math:`s`. + """ + return self.electron_pressure + self.hydrogen_pressure + + @property + @u.quantity_input + def electron_temperature(self) -> u.Unit('K'): + r""" + Electron temperature, :math:`T_e = P_e / (k_B n_e)`, as a function of :math:`s`. + """ + if hasattr(self, '_phy_data'): + return self._phy_data['electron_temperature'] + return self.electron_pressure / (const.k_B*self.electron_density) + + @property + @u.quantity_input + def hydrogen_temperature(self) -> u.Unit('K'): + r""" + Hydrogen temperature, :math:`T_H = P_H / (k_B n_H)`, as a function of :math:`s`. + """ + if hasattr(self, '_phy_data'): + return self._phy_data['hydrogen_temperature'] + return self.hydrogen_pressure / (const.k_B*self.hydrogen_density) + def spatial_average(self, quantity, bounds=None): """ Compute a spatial average of a specific quantity @@ -570,7 +627,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 @@ -616,7 +673,7 @@ def _phy_filename(self): return self.hydrad_root / 'Initial_Conditions' / 'profiles' / 'initial.amr.phy' -def add_property(name, attr, index, unit): +def _add_property(name, attr,): """ Auto-generate properties for various pieces of data """ @@ -624,25 +681,7 @@ def property_template(self): data = getattr(self, attr) if data is None: raise ValueError(f'No data available for {name}') - return u.Quantity(data[:, index], unit) + return u.Quantity(data[name]) property_template.__doc__ = f'{" ".join(name.split("_")).capitalize()} as a function of :math:`s`' property_template.__name__ = name setattr(Profile, property_template.__name__, property(property_template)) - - -properties = [ - ('coordinate', '_phy_data', 0, 'cm'), - ('velocity', '_phy_data', 1, 'cm / s'), - ('sound_speed', '_phy_data', 2, 'cm / s'), - ('electron_density', '_phy_data', 3, 'cm-3'), - ('ion_density', '_phy_data', 4, 'cm-3'), - ('electron_pressure', '_phy_data', 5, 'dyne cm-2'), - ('ion_pressure', '_phy_data', 6, 'dyne cm-2'), - ('electron_temperature', '_phy_data', 7, 'K'), - ('ion_temperature', '_phy_data', 8, 'K'), - ('electron_heat_flux', '_phy_data', 9, 'erg s-1 cm-2'), - ('ion_heat_flux', '_phy_data', 10, 'erg s-1 cm-2'), -] -properties += [(f'level_population_hydrogen_{i}', '_hstate_data', i, '') for i in range(1, 7)] -for p in properties: - add_property(*p) diff --git a/pydrad/parse/tests/test_strand.py b/pydrad/parse/tests/test_strand.py index a6de280..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', @@ -55,6 +55,15 @@ def strand(hydrad): return Strand(hydrad) +@pytest.fixture +def strand_only_amr(hydrad): + return Strand(hydrad, + read_phy=False, + read_ine=False, + read_trm=False, + read_hstate=False) + + def test_parse_initial_conditions(strand): assert hasattr(strand, 'initial_conditions') @@ -66,6 +75,23 @@ def test_has_quantity(strand, quantity): assert isinstance(getattr(p, quantity), u.Quantity) +@pytest.mark.parametrize('quantity', [ + 'electron_temperature', + 'hydrogen_temperature', + 'electron_density', + 'hydrogen_density', + 'electron_pressure', + 'hydrogen_pressure', + 'velocity', +]) +def test_no_amr_run_has_quantity(strand_only_amr, quantity): + # check that Strand falls back to deriving thermodynamic quantities from + # conserved quantities when we do not read the phy files + for p in strand_only_amr: + assert hasattr(p, quantity) + assert isinstance(getattr(p, quantity), u.Quantity) + + def test_time_arrays_same(hydrad, strand): """Check that reading time arrays different ways yields same answer""" strand2 = Strand(hydrad, read_from_cfg=True) @@ -96,15 +122,17 @@ def test_emission_measure(strand): def test_term_file_output(strand): for p in strand: # The electron energy equation's numerical viscosity term is always 0: - assert u.allclose(p.electron_numerical_viscosity, - np.zeros_like(p.electron_numerical_viscosity), - rtol=0.0, - ) + assert u.allclose( + p.electron_numerical_viscosity, + np.zeros_like(p.electron_numerical_viscosity), + rtol=0.0, + ) # The hydrogen energy equation's collision rate is never 0 at all positions: - assert not u.allclose(p.hydrogen_collisions, - np.zeros_like(p.hydrogen_collisions), - rtol=0.0, - ) + assert not u.allclose( + p.hydrogen_collisions, + np.zeros_like(p.hydrogen_collisions), + rtol=0.0, + ) def test_term_file_units(strand): assert strand[0].mass_advection.unit == u.Unit('g s-1 cm-3') diff --git a/pydrad/parse/util.py b/pydrad/parse/util.py new file mode 100644 index 0000000..23d9780 --- /dev/null +++ b/pydrad/parse/util.py @@ -0,0 +1,236 @@ +""" +Utilities related to parsing HYDRAD results +""" +import pathlib + +import astropy.table +import numpy as np +import plasmapy.particles +from pandas import read_csv + +__all__ = [ + 'read_amr_file', + 'read_phy_file', + 'read_ine_file', + 'read_trm_file', + 'read_hstate_file', +] + + +def read_amr_file(filename): + """ + 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 + # 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', + 'mass_density', + 'momentum_density', + 'electron_energy_density', + 'hydrogen_energy_density', + ] + units = { + 'grid_centers': 'cm', + 'grid_widths': 'cm', + 'mass_density': 'g cm-3', + 'momentum_density': 'g cm-2 s-1', + 'electron_energy_density': 'erg cm-3', + 'hydrogen_energy_density': 'erg cm-3', + } + table = astropy.table.QTable.read( + filename, + format='ascii', + data_start=4, + ) + # 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. + table.rename_columns( + table.colnames[:len(columns)], + columns, + ) + for column in columns: + table[column].unit = units[column] + return table + + +def read_phy_file(filename): + """ + 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', + 'hydrogen_density', + 'electron_pressure', + 'hydrogen_pressure', + 'electron_temperature', + 'hydrogen_temperature', + 'electron_heat_flux', + 'hydrogen_heat_flux', + ] + units = { + 'coordinate': 'cm', + 'velocity': 'cm / s', + 'sound_speed': 'cm / s', + 'electron_density': 'cm-3', + 'hydrogen_density': 'cm-3', + 'electron_pressure': 'dyne cm-2', + 'hydrogen_pressure': 'dyne cm-2', + 'electron_temperature': 'K', + 'hydrogen_temperature': 'K', + 'electron_heat_flux': 'erg s-1 cm-2', + 'hydrogen_heat_flux': 'erg s-1 cm-2', + } + return astropy.table.QTable.read( + filename, + format='ascii', + names=columns, + units=units, + ) + + +def read_ine_file(filename, n_s): + """ + 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. + + Parameters + ---------- + filename: path-like + n_s: `int` + """ + # TODO: clean this up somehow? I've purposefully included + # a lot of comments because the format of this file makes + # the parsing code quite opaque + with pathlib.Path(filename).open() as f: + lines = f.readlines() + # First parse all of the population fraction arrays + # NOTE: Have to calculate the number of elements we have + # computed population fractions for as we do not necessarily + # know this ahead of time + n_e = int(len(lines)/n_s - 1) + # The file is arranged in n_s groups of n_e+1 lines each where the first + # line is the coordinate and the subsequent lines are the population fraction + # for each element, with each column corresponding to an ion of that element + # First, separate by coordinate + pop_lists = [lines[i*(n_e+1)+1:(i+1)*(n_e+1)] for i in range(n_s)] + # Convert each row of each group into a floating point array + pop_lists = [[np.array(l.split(), dtype=float) for l in p] for p in pop_lists] + # NOTE: each row has Z+2 entries as the first entry is the atomic number Z + # Get these from just the first group as the number of elements is the same + # for each + Z = np.array([p[0] for p in pop_lists[0]], dtype=int) + pop_arrays = [np.zeros((n_s, z+1)) for z in Z] + for i, p in enumerate(pop_lists): + for j, line in enumerate(p): + pop_arrays[j][i, :] = line[1:] # Skip first entry, it is the atomic number + columns = [] + for z in Z: + el_name = plasmapy.particles.element_name(z) + columns += [f'{el_name}_{i+1}' for i in range(z+1)] + data = np.hstack([p for p in pop_arrays]) + return astropy.table.QTable(data=data, names=columns) + + +def read_trm_file(filename): + """ + Parse ``.trm`` files with hydrodynamic equation terms as a function of position. + + The files come in sets of 5 rows with variable number of columns: + -- Loop coordinate (1 column), and at each position: + -- Terms of mass equation (2 columns) + -- Terms of momentum equation (6 columns) + -- Terms of electron energy equation (11 columns) + -- Terms of hydrogen energy equation (11 columns) + """ + units = { + 'mass': 'g cm-3 s-1', + 'momentum': 'dyne cm-3 s-1', + 'electron': 'erg cm-3 s-1', + 'hydrogen': 'erg cm-3 s-1', + } + columns = { + 'mass': [ + 'mass_drhobydt', + 'mass_advection', + ], + 'momentum': [ + 'momentum_drho_vbydt', + 'momentum_advection', + 'momentum_pressure_gradient', + 'momentum_gravity', + 'momentum_viscous_stress', + 'momentum_numerical_viscosity', + ], + 'electron': [ + 'electron_dTEKEbydt', + 'electron_enthalpy', + 'electron_conduction', + 'electron_gravity', + 'electron_collisions', + 'electron_heating', + 'electron_radiative_loss', + 'electron_electric_field', + 'electron_viscous_stress', + 'electron_numerical_viscosity', + 'electron_ionization' + ], + 'hydrogen': [ + 'hydrogen_dTEKEbydt', + 'hydrogen_enthalpy', + 'hydrogen_conduction', + 'hydrogen_gravity', + 'hydrogen_collisions', + 'hydrogen_heating', + 'hydrogen_radiative_loss', + 'hydrogen_electric_field', + 'hydrogen_viscous_stress', + 'hydrogen_numerical_viscosity', + 'hydrogen_ionization', + ] + } + tables = [] + for i,(k,v) in enumerate(columns.items()): + tables.append( + astropy.table.QTable.from_pandas( + read_csv( + filename, + sep='\t', + header=None, + names=v, + skiprows=lambda x: x % 5 != (i+1) + ), + units={c: units[k] for c in v} + ) + ) + return astropy.table.hstack(tables) + + +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( + filename, + format='ascii', + names=columns, + ) diff --git a/pydrad/util/__init__.py b/pydrad/util/__init__.py new file mode 100644 index 0000000..7691e8d --- /dev/null +++ b/pydrad/util/__init__.py @@ -0,0 +1,3 @@ +""" +Utility functions useful to the whole `pydrad` package +""" diff --git a/pydrad/util/constants.py b/pydrad/util/constants.py new file mode 100644 index 0000000..eaff23f --- /dev/null +++ b/pydrad/util/constants.py @@ -0,0 +1,36 @@ +""" +Constants used throughout `pydrad` +""" +from astropy.constants import Constant + +__all__ = [ + 'gamma', + 'm_avg_ion', +] + + +gamma = Constant( + abbrev='gamma', + name='heat capacity ratio of a monatomic ideal gas', + value=5/3, + unit='', + uncertainty=0, + reference='https://en.wikipedia.org/wiki/Heat_capacity_ratio', + system=None, +) + + +# 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_avg_ion = m_H + \frac{n_{He}}{n_{H}} m_{He} +m_avg_ion = Constant( + abbrev='m_avg_ion', + name='average ion mass for a H-He plasma', + value=2.171e-24, + unit='g', + uncertainty=0, + reference='https://github.com/rice-solar-physics/HYDRAD/blob/master/Resources/source/constants.h', + system='cgs', +) 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='--' ) diff --git a/setup.cfg b/setup.cfg index 1aec17c..f94b50b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -66,6 +66,7 @@ filterwarnings = ignore:numpy.ndarray size changed:RuntimeWarning ignore:The unit 'erg' has been deprecated in the VOUnit standard.:astropy.units.core.UnitsWarning ignore:datetime.datetime.utcnow\(\) is deprecated and scheduled for removal in a future version.:DeprecationWarning + ignore:pkg_resources is deprecated as an API.:DeprecationWarning [coverage:run] omit =