Add tests comparing with HYDRAD results (#85)

examples/plot.py

@@ -1,13 +1,16 @@
 """
 Visualize ebtel++ example results
 """
+import astropy.units as u
+from astropy.visualization import quantity_support
 import matplotlib.pyplot as plt
 from matplotlib import gridspec
 from matplotlib.ticker import MaxNLocator
 import seaborn
 
 
 def make_figure(results, filename):
+    quantity_support()
     seaborn.set_context('notebook', font_scale=1.1)
     seaborn.set(font='serif')
     seaborn.set_style("white", {
@@ -59,8 +62,8 @@ def make_figure(results, filename):
     ax2.legend(loc='best')
     ax4.set_xlabel(r'$T$ (K)')
     ax4.set_ylabel(r'$\mathrm{DEM}$ (cm$^{-5}$ K$^{-1}$)')
-    ax4.set_xlim([10**(4.5), 10**(7.5)])
-    ax4.set_ylim([10**(20.0), 10**(23.5)])
+    ax4.set_xlim([10**(4.5), 10**(7.5)]*u.K)
+    ax4.set_ylim([10**(20.0), 10**(23.5)]*u.Unit('cm-5 K-1'))
     ax4.set_xscale('log')
     ax4.set_yscale('log')
     ax4.legend(loc='best')

examples/util.py

@@ -9,6 +9,7 @@
 import xml.etree.ElementTree as ET
 import xml.dom.minidom as xdm
 
+import astropy.units as u
 import numpy as np
 
 __all__ = ['run_ebtel', 'read_xml', 'write_xml']
@@ -48,14 +49,14 @@ def run_ebtel(config, ebtel_dir):
         data = np.loadtxt(results_filename)
 
         results = {
-            'time': data[:, 0],
-            'electron_temperature': data[:, 1],
-            'ion_temperature': data[:, 2],
-            'density': data[:, 3],
-            'electron_pressure': data[:, 4],
-            'ion_pressure': data[:, 5],
-            'velocity': data[:, 6],
-            'heat': data[:, 7],
+            'time': data[:, 0]*u.s,
+            'electron_temperature': data[:, 1]*u.K,
+            'ion_temperature': data[:, 2]*u.K,
+            'density': data[:, 3]*u.cm**(-3),
+            'electron_pressure': data[:, 4]*u.dyne/(u.cm**2),
+            'ion_pressure': data[:, 5]*u.dyne/(u.cm**2),
+            'velocity': data[:, 6]*u.cm/u.s,
+            'heat': data[:, 7]*u.erg/(u.cm**3*u.s),
         }
 
         results_dem = {}
@@ -65,9 +66,9 @@ def run_ebtel(config, ebtel_dir):
             results_dem['dem_corona'] = np.loadtxt(
                 config['output_filename'] + '.dem_corona')
             # The first row of both is the temperature bins
-            results_dem['dem_temperature'] = results_dem['dem_tr'][0, :]
-            results_dem['dem_tr'] = results_dem['dem_tr'][1:, :]
-            results_dem['dem_corona'] = results_dem['dem_corona'][1:, :]
+            results_dem['dem_temperature'] = results_dem['dem_tr'][0, :]*u.K
+            results_dem['dem_tr'] = results_dem['dem_tr'][1:, :]*u.Unit('cm-5 K-1')
+            results_dem['dem_corona'] = results_dem['dem_corona'][1:, :]*u.Unit('cm-5 K-1')
 
     return {**results, **results_dem}
 

requirements/requirements-test.txt

@@ -1,3 +1,5 @@
 -r requirements.txt
 pytest
-hissw
+hissw
+h5py
+astropy

tests/helpers.py

@@ -2,17 +2,23 @@
 Helper functions for tests
 """
 import os
+import pathlib
 import sys
 
+import astropy.units as u
+import h5py
 import numpy as np
 
-TOPDIR = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
-sys.path.append(os.path.join(TOPDIR, 'examples'))
+TEST_DIR = pathlib.Path(__file__).parent.resolve()
+DATA_DIR = TEST_DIR / 'data'
+REPO_DIR = TEST_DIR.parent
+sys.path.append(str(REPO_DIR / 'examples'))
 from util import run_ebtel
 
 
 def run_ebtelplusplus(config):
-    return run_ebtel(config, TOPDIR)
+    res = run_ebtel(config, REPO_DIR)
+    return res
 
 
 def generate_idl_test_data(ebtel_idl_path, config):
@@ -55,6 +61,7 @@ def generate_idl_test_data(ebtel_idl_path, config):
 
 def read_idl_test_data(data_filename, ebtel_idl_path, config):
     varnames = ['time', 'temperature', 'density', 'pressure', 'velocity']
+    varunits = ['s', 'K', 'cm-3', 'dyne cm-2', 'cm s-1']
     # Generate and save if it does not exist
     if not os.path.isfile(data_filename):
         data = generate_idl_test_data(ebtel_idl_path, config)
@@ -64,4 +71,18 @@ def read_idl_test_data(data_filename, ebtel_idl_path, config):
         np.savetxt(data_filename, data_array)
     # Load data into a dictionary
     data = np.loadtxt(data_filename)
-    return {v: data[:, i] for i, v in enumerate(varnames)}
+    return {v: u.Quantity(data[:, i], vu) for i, (v,vu) in enumerate(zip(varnames, varunits))}
+
+
+def read_hydrad_test_data(data_filename, tau, heating):
+    data = {}
+    with h5py.File(data_filename, 'r') as hf:
+        grp = hf[f'/{heating}/tau{tau:.0f}']
+        data['time'] = u.Quantity(hf['time'], hf['time'].attrs['unit'])
+        data['electron_temperature'] = u.Quantity(grp['electron_temperature'],
+                                                  grp[f'electron_temperature'].attrs['unit'])
+        data['ion_temperature'] = u.Quantity(grp['ion_temperature'],
+                                             grp['ion_temperature'].attrs['unit'])
+        data['density'] = u.Quantity(grp['density'],
+                                     grp['density'].attrs['unit'])
+    return data

tests/test_compare_hydrad.py

@@ -0,0 +1,76 @@
+"""
+Compare output of EBTEL IDL and ebtel++
+"""
+import pytest
+from collections import OrderedDict
+
+import astropy.units as u
+from scipy.interpolate import interp1d
+
+from .helpers import DATA_DIR, read_hydrad_test_data, run_ebtelplusplus
+
+
+@pytest.fixture
+def base_config():
+    base_config = {
+        'total_time': 5e3,
+        'tau': 0.1,
+        'tau_max': 10,
+        'loop_length': 4e9,
+        'saturation_limit': 1,
+        'force_single_fluid': False,
+        'use_c1_loss_correction': True,
+        'use_c1_grav_correction': True,
+        'use_flux_limiting': True,
+        'use_power_law_radiative_losses': True,
+        'calculate_dem': False,
+        'save_terms': False,
+        'use_adaptive_solver': True,
+        'adaptive_solver_error': 1e-8,
+        'adaptive_solver_safety': 0.5,
+        'c1_cond0': 6.0,
+        'c1_rad0': 0.6,
+        'helium_to_hydrogen_ratio': 0.075,
+        'surface_gravity': 1.0,
+        'heating': OrderedDict({
+            'partition': 1.,
+            'background': 3.5e-5,}),
+    }
+    return base_config
+
+
+@pytest.mark.parametrize('tau', [200.0, 500.0])
+@pytest.mark.parametrize('heating_type', ['single', 'electron', 'ion'])
+def test_compare_hydrad_single_event_peak_values(base_config, tau, heating_type):
+    config = base_config.copy()
+    if heating_type == 'single':
+        config['force_single_fluid'] = True
+        config['heating']['partition'] = 0.5
+    elif heating_type == 'electron':
+        config['force_single_fluid'] = False
+        config['heating']['partition'] = 1.0
+    elif heating_type == 'ion':
+        config['force_single_fluid'] = False
+        config['heating']['partition'] = 0.0
+    else:
+        raise ValueError('Unrecognized heating type')
+    total_energy = 10.0
+    config['heating']['events'] = [
+        {'event': {
+            'magnitude': 2 * total_energy/tau,
+            'rise_start': 0.,
+            'rise_end': tau/2.,
+            'decay_start': tau/2.,
+            'decay_end': tau,
+        }}
+    ]
+    r_ebtel = run_ebtelplusplus(config)
+    r_hydrad = read_hydrad_test_data(DATA_DIR / 'hydrad_results.h5', tau, heating_type)
+    # Require all quantities at the peak to be <20% in accordance with the comparisons
+    # done in Barnes et al. (2016a)
+    for name in ['electron_temperature', 'ion_temperature', 'density']:
+        f_interp = interp1d(r_ebtel['time'].to_value('s'), r_ebtel[name].to_value(r_hydrad[name].unit), 
+                            fill_value='extrapolate')
+        x_ebtel_interp = u.Quantity(f_interp(r_hydrad['time'].to_value('s')), r_hydrad[name].unit)
+        i_peak = r_hydrad[name].argmax()
+        assert u.allclose(x_ebtel_interp[i_peak], r_hydrad[name][i_peak], rtol=0.20)

tests/test_compare_idl.py

@@ -1,13 +1,12 @@
 """
 Compare output of EBTEL IDL and ebtel++
 """
-import os
 import pytest
 from collections import OrderedDict
 
-import numpy as np
+import astropy.units as u
 
-from .helpers import read_idl_test_data, run_ebtelplusplus
+from .helpers import DATA_DIR, read_idl_test_data, run_ebtelplusplus
 
 
 @pytest.fixture
@@ -47,18 +46,11 @@ def test_compare_idl_single_event(base_config, tmpdir, ebtel_idl_path):
             'rise_start': 0.0, 'rise_end': 100.0, 'decay_start': 100.0,
             'decay_end': 200.0, 'magnitude': 0.1}}]
     r_cpp = run_ebtelplusplus(config)
-    r_idl = read_idl_test_data('tests/data/idl_single_event.json',
+    r_idl = read_idl_test_data(DATA_DIR / 'data' / 'idl_single_event.json',
                                ebtel_idl_path, config)
-    # Temperature
-    assert np.allclose(r_cpp['electron_temperature'], r_idl['temperature'],
-                       atol=0., rtol=1e-2)
-    assert np.allclose(r_cpp['ion_temperature'], r_idl['temperature'],
-                       atol=0., rtol=1e-2)
-    # Density
-    assert np.allclose(r_cpp['density'], r_idl['density'], atol=0., rtol=1e-2)
-    # Pressure
-    assert np.allclose(r_cpp['electron_pressure']+r_cpp['ion_pressure'],
-                       r_idl['pressure'], atol=0., rtol=1e-2)
-    # Velocity
-    assert np.allclose(r_cpp['velocity'], r_idl['velocity'], atol=0.,
-                       rtol=1e-2)
+    assert u.allclose(r_cpp['electron_temperature'], r_idl['temperature'], rtol=1e-2)
+    assert u.allclose(r_cpp['ion_temperature'], r_idl['temperature'], rtol=1e-2)
+    assert u.allclose(r_cpp['density'], r_idl['density'], rtol=1e-2)
+    assert u.allclose(r_cpp['electron_pressure']+r_cpp['ion_pressure'], r_idl['pressure'], 
+                      rtol=1e-2)
+    assert u.allclose(r_cpp['velocity'], r_idl['velocity'], rtol=1e-2)

tests/test_single_fluid.py

@@ -4,8 +4,8 @@
 import os
 from collections import OrderedDict
 
+import astropy.units as u
 import pytest
-import numpy as np
 
 from .helpers import run_ebtelplusplus
 
@@ -39,27 +39,20 @@ def base_config():
     return base_config
 
 
-def test_single_fluid_gentle(base_config, tmpdir):
+@pytest.mark.parametrize(['rise_start', 'rise_end', 'decay_start', 'decay_end', 'magnitude'], [
+    (0, 250, 1000, 2000, 0.005),  # gentle case
+    (0, 100, 100, 200, 0.1),  # more impulsive case
+])
+def test_single_fluid_gentle(base_config, rise_start, rise_end, decay_start, decay_end, magnitude):
     base_config['heating']['events'] = [
-        {'event': {'rise_start': 0.0, 'rise_end': 250.0, 'decay_start': 1000.0,
-                   'decay_end': 2000.0, 'magnitude': 0.005}}]
+        {'event': {
+            'rise_start': rise_start,
+            'rise_end': rise_end,
+            'decay_start': decay_start,
+            'decay_end': decay_end,
+            'magnitude': magnitude,
+        }}
+    ]
     results = run_ebtelplusplus(base_config)
-    # Temperature
-    assert np.allclose(results['electron_temperature'],
-                       results['ion_temperature'], atol=0., rtol=1e-10)
-    # Pressure
-    assert np.allclose(results['electron_pressure'], results['ion_pressure'],
-                       atol=0., rtol=1e-10)
-
-
-def test_single_fluid_impulsive(base_config, tmpdir):
-    base_config['heating']['events'] = [
-        {'event': {'rise_start': 0.0, 'rise_end': 100.0, 'decay_start': 100.0,
-                   'decay_end': 200.0, 'magnitude': 0.1}}]
-    results = run_ebtelplusplus(base_config)
-    # Temperature
-    assert np.allclose(results['electron_temperature'],
-                       results['ion_temperature'], atol=0., rtol=1e-10)
-    # Pressure
-    assert np.allclose(results['electron_pressure'], results['ion_pressure'],
-                       atol=0., rtol=1e-10)
+    assert u.allclose(results['electron_temperature'], results['ion_temperature'], rtol=1e-10)
+    assert u.allclose(results['electron_pressure'], results['ion_pressure'], rtol=1e-10)

tests/test_solver.py

@@ -9,7 +9,7 @@
 from .helpers import run_ebtelplusplus
 
 
-@pytest.fixture
+@pytest.fixture(scope='module')
 def base_config():
     base_config = {
         'total_time': 5e3,
@@ -45,51 +45,32 @@ def base_config():
     return base_config
 
 
-@pytest.fixture
-def adaptive_results(base_config, tmpdir):
+@pytest.fixture(scope='module')
+def adaptive_results(base_config):
     config = base_config.copy()
     config['use_adaptive_solver'] = True
     return run_ebtelplusplus(config)    
 
 
-@pytest.fixture
-def static_results(base_config, tmpdir):
+@pytest.fixture(scope='module')
+def static_results(base_config):
     config = base_config.copy()
     config['use_adaptive_solver'] = False
     return run_ebtelplusplus(config)
 
-
-def test_temperature_equal_adaptive_static(adaptive_results, static_results):
-    t_static = static_results['time']
-    t_adaptive = adaptive_results['time']
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['electron_temperature'])
+@pytest.mark.parametrize(['name', 'atol'], [
+    ('electron_temperature', 1e4),
+    ('ion_temperature', 1e4),
+    ('density', 1e7),
+    ('electron_pressure', 1e-4),
+    ('ion_pressure', 1e-4),
+    ('velocity', 1e4),
+])
+def test_quantities_equal_adaptive_static(adaptive_results, static_results, name, atol):
+    t_static = static_results['time'].to_value('s')
+    t_adaptive = adaptive_results['time'].to_value('s')
+    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results[name].value)
     # NOTE: Skip the first 5 steps b/c there is always one anomalous point that gives
     # an error > 10%; due to static case not rising fast enough
-    assert np.allclose(adapt_interp[5:],
-                       static_results['electron_temperature'][5:], rtol=1e-2, atol=1e4)
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['ion_temperature'])
-    assert np.allclose(adapt_interp[5:], static_results['ion_temperature'][5:], rtol=1e-2, atol=1e4)
-
-
-def test_density_equal_adaptive_static(adaptive_results, static_results):
-    t_static = static_results['time']
-    t_adaptive = adaptive_results['time']
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['density'])
-    assert np.allclose(adapt_interp[5:], static_results['density'][5:], rtol=1e-2, atol=1e7)
-
-
-def test_pressure_equal_adaptive_static(adaptive_results, static_results):
-    t_static = static_results['time']
-    t_adaptive = adaptive_results['time']
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['electron_pressure'])
-    assert np.allclose(adapt_interp[5:], static_results['electron_pressure'][5:],
-                       rtol=1e-2, atol=1e-4)
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['ion_pressure'])
-    assert np.allclose(adapt_interp[5:], static_results['ion_pressure'][5:], rtol=1e-2, atol=1e-4)
-
-
-def test_velocity_equal_adaptive_static(adaptive_results, static_results):
-    t_static = static_results['time']
-    t_adaptive = adaptive_results['time']
-    adapt_interp = np.interp(t_static, t_adaptive, adaptive_results['velocity'])
-    assert np.allclose(adapt_interp[5:], static_results['velocity'][5:], rtol=1e-2, atol=1e4)
+    assert np.allclose(adapt_interp[5:], static_results[name][5:].to_value(adaptive_results[name].unit),
+                       rtol=1e-2, atol=atol)