Skip to content

Commit

Permalink
Fully parse .trm files (#169)
Browse files Browse the repository at this point in the history 
* fully parse .trm files

* add new variables to tests

* swap if for elif

* swap else for elif -- i % 0 is not tested

* revamp to remove all indices

* only call len() one time for each

* change - to +

* add some unit tests

* strand -> profile

* import numpy in parse/tests/test_strand'

* add tolerance to allclose

* switch gravity to collisions -- gravity is 0 at time zero

* small typo

* remove absolute tolerance

---------

Co-authored-by: Jeffrey Reep <reep@atrcw15.IfA.Hawaii.Edu>
Co-authored-by: Jeffrey Reep <reep@kahiku.local>
Co-authored-by: Jeffrey Reep <reep@atrcl164.IfA.Hawaii.edu>
  • Loading branch information
@jwreep
4 people committed Jun 17, 2024
1 parent 3ed98ff commit dd41f6a
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 42 deletions.
105 changes: 68 additions & 37 deletions pydrad/parse/parse.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
import matplotlib.pyplot as plt
import numpy as np
import plasmapy.particles
from pandas import read_csv
from scipy.interpolate import splev, splrep

from pydrad import log
Expand Down Expand Up @@ -358,42 +359,72 @@ def _read_amr(self):
def _read_trm(self):
"""
Parse the equation terms file
"""
try:
with self._trm_filename.open() as f:
lines = f.readlines()
except FileNotFoundError:
log.debug(f'{self._trm_filename} not found')
return
n_elements = int(len(lines)/5)
self._trm_data = np.zeros([n_elements, 3])

# The files come in sets of 5 rows
# -- Loop coordinate, and at each one:
# -- Terms of mass equation
# -- Terms of momentum equation
# -- Terms of electron energy equation
# -- Terms of hydrogen energy equation
# Right now, we only read 3 values from this:
# the electron heating, hydrogen heating,
# and bolometric radiative losses
for i in range(len(lines)):
j = int(i/5)
line = lines[i].strip().split()
# Electron heating and radiative loss terms from the
# electron energy equation
if i % 5 == 3:
self._trm_data[j, 0] = float(line[5])
self._trm_data[j, 1] = float(line[6])
# Hydrogen heating term from the hydrogen energy
# equation
if i % 5 == 4:
self._trm_data[j, 2] = float(line[5])

properties = [('electron_heating_term', '_trm_data', 0, 'erg cm-3 s-1'),
('hydrogen_heating_term', '_trm_data', 2, 'erg cm-3 s-1'),
('radiative_loss_term', '_trm_data', 1, 'erg cm-3 s-1')]
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)
Expand Down Expand Up @@ -609,8 +640,8 @@ def property_template(self):
('ion_pressure', '_phy_data', 6, 'dyne cm-2'),
('electron_temperature', '_phy_data', 7, 'K'),
('ion_temperature', '_phy_data', 8, 'K'),
('electron_conduction', '_phy_data', 9, 'erg s-1 cm-2'),
('ion_conduction', '_phy_data', 10, 'erg s-1 cm-2'),
('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:
Expand Down
51 changes: 46 additions & 5 deletions pydrad/parse/tests/test_strand.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""
import astropy.units as u
import h5py
import numpy as np
import pytest

from pydrad.parse import Strand
Expand All @@ -19,20 +20,41 @@
'electron_pressure',
'ion_pressure',
'velocity',
'electron_heating_term',
'hydrogen_heating_term',
'radiative_loss_term',
'sound_speed',
'electron_heat_flux',
'ion_heat_flux',
'mass_drhobydt',
'mass_advection',
'momentum_drho_vbydt',
'momentum_advection',
'momentum_pressure_gradient',
'momentum_gravity',
'momentum_viscous_stress',
'momentum_numerical_viscosity',
'electron_dTEKEbydt',
'electron_enthalpy',
'electron_conduction',
'ion_conduction',
'electron_collisions',
'electron_heating',
'electron_radiative_loss',
'electron_electric_field',
'electron_ionization',
'hydrogen_dTEKEbydt',
'hydrogen_enthalpy',
'hydrogen_conduction',
'hydrogen_gravity',
'hydrogen_collisions',
'hydrogen_heating',
'hydrogen_electric_field',
'hydrogen_viscous_stress',
'hydrogen_numerical_viscosity',
]


@pytest.fixture
def strand(hydrad):
return Strand(hydrad)


def test_parse_initial_conditions(strand):
assert hasattr(strand, 'initial_conditions')

Expand Down Expand Up @@ -70,3 +92,22 @@ def test_emission_measure(strand):
assert isinstance(em, u.Quantity)
assert isinstance(bins, u.Quantity)
assert len(bins) == len(em) + 1

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,
)
# 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,
)

def test_term_file_units(strand):
assert strand[0].mass_advection.unit == u.Unit('g s-1 cm-3')
assert strand[0].momentum_gravity.unit == u.Unit('dyne s-1 cm-3')
assert strand[0].electron_viscous_stress.unit == u.Unit('erg s-1 cm-3')
assert strand[0].hydrogen_collisions.unit == u.Unit('erg s-1 cm-3')

0 comments on commit dd41f6a

Please sign in to comment.