Compute and save the ISRF in every cell

docs/setup/setup_conf.rst

@@ -195,6 +195,58 @@ only after last iteration), or ``none`` (do not output). The default is to
 output only the last iteration of ``specific_energy``. To find out how to view
 these values, see :doc:`../postprocessing/postprocessing`
 
+Frequency-resolved specific energy
+----------------------------------
+
+In addition to the total specific energy absorbed in each cell
+(``output_specific_energy``), Hyperion can save the specific energy absorbed as
+a function of frequency. This is controlled in the same way as the other output
+quantities::
+
+    m.conf.output.output_specific_energy_spectrum = 'all'    # every iteration
+                                              # 'last'  -> only final iteration
+                                              # 'none'  -> not computed or saved (default)
+
+When it is not ``'none'``, two extra arrays are written out:
+
+* ``specific_energy_spectrum`` -- the specific energy absorbed in each cell as a
+  function of frequency, in erg/s/g.
+* ``specific_energy_spectrum_frequencies`` -- the frequencies (in Hz) corresponding to
+  the leading axis of ``specific_energy_spectrum``. These are bin centers, not edges,
+  so this array has exactly the same length as that axis (one entry per bin).
+
+This is the *absorbed* (deposited) energy spectrum, not the mean intensity: each
+contribution is weighted by the dust absorption opacity, so it is proportional
+to :math:`\kappa_\nu J_\nu`. To recover the mean intensity :math:`J_\nu` (the
+radiation field), divide by the dust absorption opacity at each frequency.
+Summed over frequency, ``specific_energy_spectrum`` recovers the total
+``specific_energy``.
+
+By default the binning uses the frequency grid of the first dust type. You can
+instead provide your own frequency grid (in Hz), independent of the dust
+properties::
+
+    import numpy as np
+    m.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 100))
+
+Photons are binned to the nearest of these frequencies in log space (so the
+supplied values act as bin centers). This works for all
+grid types, including AMR and Voronoi.
+
+.. warning:: The binning has no outer edges: the first and last bins collect
+             *all* photons below and above the outermost bin centers
+             respectively, however far away in frequency. This keeps the spectrum
+             energy-conserving (no photons are dropped), but it means a grid that
+             does not span the full frequency range will pile up out-of-range
+             flux in its end bins. When supplying a custom grid, make sure it
+             brackets the full range of frequencies present in the simulation
+             (the default dust-based grid already does this).
+
+``specific_energy_spectrum`` can be retrieved like other grid quantities, as an array
+with an extra leading frequency axis::
+
+    senu = m.get_quantities()['specific_energy_spectrum']
+
 Advanced parameters
 -------------------
 

hyperion/conf/conf_files.py

@@ -18,6 +18,7 @@ def __init__(self):
         self.output_density = 'none'
         self.output_density_diff = 'none'
         self.output_specific_energy = 'last'
+        self.output_specific_energy_spectrum = 'none'
         self.output_n_photons = 'none'
         self._freeze()
 
@@ -27,13 +28,18 @@ def read(cls, group):
         self.output_density = group.attrs['output_density'].decode('utf-8')
         self.output_density_diff = group.attrs['output_density_diff'].decode('utf-8')
         self.output_specific_energy = group.attrs['output_specific_energy'].decode('utf-8')
+        if 'output_specific_energy_spectrum' in group.attrs:
+            self.output_specific_energy_spectrum = group.attrs['output_specific_energy_spectrum'].decode('utf-8')
+        else:
+            self.output_specific_energy_spectrum = 'none'
         self.output_n_photons = group.attrs['output_n_photons'].decode('utf-8')
         return self
 
     def write(self, group):
         group.attrs['output_density'] = np.bytes_(self.output_density.encode('utf-8'))
         group.attrs['output_density_diff'] = np.bytes_(self.output_density_diff.encode('utf-8'))
         group.attrs['output_specific_energy'] = np.bytes_(self.output_specific_energy.encode('utf-8'))
+        group.attrs['output_specific_energy_spectrum'] = np.bytes_(self.output_specific_energy_spectrum.encode('utf-8'))
         group.attrs['output_n_photons'] = np.bytes_(self.output_n_photons.encode('utf-8'))
 
 
@@ -52,6 +58,8 @@ def __init__(self):
         self.set_max_reabsorptions(1000000)
         self.set_pda(False)
         self.set_mrw(False)
+        self.specific_energy_spectrum_frequencies = None
+
         self.set_convergence(False)
         self.set_kill_on_absorb(False)
         self.set_kill_on_scatter(False)
@@ -391,6 +399,43 @@ def _read_pda(self, group):
     def _write_pda(self, group):
         group.attrs['pda'] = bool2str(self.pda)
 
+    def _read_specific_energy_spectrum_frequencies(self, group):
+        if 'specific_energy_spectrum_frequencies' in group:
+            self.specific_energy_spectrum_frequencies = np.array(group['specific_energy_spectrum_frequencies']['nu'])
+        else:
+            self.specific_energy_spectrum_frequencies = None
+
+    def _write_specific_energy_spectrum_frequencies(self, group):
+        if self.specific_energy_spectrum_frequencies is not None:
+            group.create_dataset('specific_energy_spectrum_frequencies',
+                                 data=np.array(list(zip(self.specific_energy_spectrum_frequencies)),
+                                               dtype=[('nu', float)]))
+
+    def set_specific_energy_spectrum_frequencies(self, frequencies):
+
+        '''
+        Set the frequency grid onto which the frequency-resolved specific energy
+        (``specific_energy_spectrum``) is binned.
+
+        This is only relevant if ``conf.output.output_specific_energy_spectrum`` is set
+        to ``'all'`` or ``'last'``. If this method is not called, the frequency
+        grid of the first dust type is used. Photons are binned to the nearest
+        frequency in log space, so the supplied values are bin centers (not
+        edges) and ``specific_energy_spectrum`` has one entry per supplied frequency.
+
+        Parameters
+        ----------
+        frequencies : iterable of float
+            The frequencies (in Hz) onto which to bin ``specific_energy_spectrum``.
+        '''
+        frequencies = np.asarray(frequencies, dtype=float)
+        if frequencies.ndim != 1 or frequencies.size < 1:
+            raise ValueError("frequencies should be a 1-d array of at least one value")
+        if np.any(frequencies <= 0.):
+            raise ValueError("frequencies should be positive (in Hz)")
+        self.specific_energy_spectrum_frequencies = np.sort(frequencies)
+
+
     def set_mrw(self, mrw, gamma=1.0, inter_max=1000, warn=True):
         '''
         Set whether to use the Modified Random Walk (MRW) approximation
@@ -719,6 +764,7 @@ def read_run_conf(self, group):  # not a class method because inherited
         self._read_max_reabsorptions(group)
         self._read_pda(group)
         self._read_mrw(group)
+        self._read_specific_energy_spectrum_frequencies(group)
         self._read_convergence(group)
         self._read_kill_on_absorb(group)
         self._read_kill_on_scatter(group)
@@ -747,6 +793,7 @@ def write_run_conf(self, group):
         self._write_max_reabsorptions(group)
         self._write_pda(group)
         self._write_mrw(group)
+        self._write_specific_energy_spectrum_frequencies(group)
         self._write_convergence(group)
         self._write_kill_on_absorb(group)
         self._write_kill_on_scatter(group)

hyperion/conf/tests/test_conf_io.py

@@ -4,6 +4,7 @@
 
 from itertools import product
 
+import numpy as np
 from numpy.testing import assert_equal
 import pytest
 
@@ -13,7 +14,8 @@
 
 @pytest.mark.parametrize(('attribute', 'value'),
                          list(product(['output_density', 'output_density_diff',
-                                       'output_specific_energy', 'output_n_photons'],
+                                       'output_specific_energy', 'output_specific_energy_spectrum',
+                                       'output_n_photons'],
                                       ['none', 'last', 'all'])))
 def test_io_output_conf(attribute, value):
     o1 = OutputConf()
@@ -249,6 +251,16 @@ def test_io_run_conf_mrw(value):
     r2.read_run_conf(v)
     assert r2.mrw == r1.mrw
 
+def test_io_run_conf_specific_energy_spectrum_frequencies():
+    r1 = RunConf()
+    r1.set_specific_energy_spectrum_frequencies(np.logspace(11., 16., 10))
+    r1.set_n_photons(1, 2)
+    v = virtual_file()
+    r1.write_run_conf(v)
+    r2 = RunConf()
+    r2.read_run_conf(v)
+    np.testing.assert_allclose(r2.specific_energy_spectrum_frequencies,
+                               r1.specific_energy_spectrum_frequencies)
 
 @pytest.mark.parametrize(('value'), [False, True])
 def test_io_run_conf_convergence(value):

hyperion/grid/cartesian_grid.py

@@ -279,6 +279,8 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
+                if quantity == 'specific_energy_spectrum_frequencies':
+                    continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
                     if array.ndim == 4:  # if array is 4D, it is a list of 3D arrays

hyperion/grid/cylindrical_polar_grid.py

@@ -307,6 +307,8 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
+                if quantity == 'specific_energy_spectrum_frequencies':
+                    continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
                     if array.ndim == 4:  # if array is 4D, it is a list of 3D arrays

hyperion/grid/grid_helpers.py

@@ -22,7 +22,7 @@ def single_grid_dims(data, ndim=3):
     '''
 
     if type(data) in [list, tuple]:
-
+        
         n_pop = len(data)
         shape = None
         for item in data:
@@ -34,14 +34,21 @@ def single_grid_dims(data, ndim=3):
         if shape is not None and len(shape) != ndim:
             raise ValueError("Grids should be %i-dimensional" % ndim)
 
-    elif isinstance(data, np.ndarray):
 
+
+    elif isinstance(data, np.ndarray):
         if data.ndim == ndim:
             n_pop = None
             shape = data.shape
         elif data.ndim == ndim + 1:
             n_pop = data.shape[0]
             shape = data[0].shape
+
+        elif data.ndim == ndim + 2:
+            # Frequency-resolved quantity, stored as (n_freq, n_pop, *grid)
+            n_pop = data.shape[1]
+            shape = data.shape[-ndim:]
+
         else:
             raise Exception("Unexpected number of dimensions: %i" % data.ndim)
 
@@ -56,6 +63,9 @@ def single_grid_dims(data, ndim=3):
         elif len(shape) == ndim + 1:
             n_pop = shape[0]
             shape = shape[1:]
+        elif len(shape) == ndim + 2:
+            n_pop = shape[1]
+            shape = shape[-ndim:]
         else:
             raise Exception("Unexpected number of dimensions: %i" % len(shape))
     else:

hyperion/grid/octree_grid.py

@@ -372,6 +372,8 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
+                if quantity == 'specific_energy_spectrum_frequencies':
+                    continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
                     if array.ndim == 2:  # if array is 2D, it is a list of 1D arrays

hyperion/grid/spherical_polar_grid.py

@@ -317,6 +317,8 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
+                if quantity == 'specific_energy_spectrum_frequencies':
+                    continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
                     if array.ndim == 4:  # if array is 4D, it is a list of 3D arrays

hyperion/grid/voronoi_grid.py

@@ -396,6 +396,8 @@ def read_quantities(self, group, quantities='all'):
         # Read in physical quantities
         if quantities is not None:
             for quantity in group:
+                if quantity == 'specific_energy_spectrum_frequencies':
+                    continue  # per-frequency metadata, not a per-cell grid quantity
                 if quantities == 'all' or quantity in quantities:
                     array = np.array(group[quantity])
                     if array.ndim == 2:  # if array is 2D, it is a list of 1D arrays

hyperion/grid/yt3_wrappers.py

@@ -79,6 +79,9 @@ def amr_grid_to_yt_stream(levels, dust_id=0):
             grid_dict['level'] = ilevel
 
             for field in grid.quantities:
+                if not isinstance(grid.quantities[field], list):
+                    logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field))
+                    continue
                 grid_dict[('gas', field)] = grid.quantities[field][dust_id].transpose()
 
             grid_data.append(grid_dict)
@@ -158,6 +161,9 @@ def octree_grid_to_yt_stream(grid, dust_id=0):
 
     quantities = {}
     for field in grid.quantities:
+        if not isinstance(grid.quantities[field], list):
+            logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field))
+            continue
         quantities[('gas', field)] = np.atleast_2d(grid.quantities[field][dust_id][order][~refined]).transpose()
 
     bbox = np.array([[xmin, xmax], [ymin, ymax], [zmin, zmax]])
@@ -180,6 +186,9 @@ def cartesian_grid_to_yt_stream(grid, xmin, xmax, ymin, ymax, zmin, zmax, dust_i
     # Make data dict which should contain (array, unit) tuples
     data = {}
     for field in grid.quantities:
+        if not isinstance(grid.quantities[field], list):
+            logger.warning("Skipping frequency-resolved quantity '{0}' in yt export (select a single frequency first)".format(field))
+            continue
         data[field] = grid.quantities[field][dust_id].transpose(), ''
 
     # Load cartesian grid into yt

hyperion/model/model.py

@@ -223,7 +223,7 @@ def use_geometry(self, filename):
         # Close the file
         f.close()
 
-    def use_quantities(self, filename, quantities=['density', 'specific_energy'],
+    def use_quantities(self, filename, quantities=['density', 'specific_energy', 'specific_energy_spectrum'],
                        use_minimum_specific_energy=True, use_dust=True, copy=True,
                        only_initial=False):
         '''
@@ -297,6 +297,13 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'],
                 else:
                     quantities_path['specific_energy'] = last_iteration
 
+            if 'specific_energy_spectrum' in quantities:
+                if only_initial or last_iteration is None:
+                    if 'specific_energy_spectrum' in f['/Input/Grid/Quantities']:
+                        quantities_path['specific_energy_spectrum'] = '/Input/Grid/Quantities'
+                elif 'specific_energy_spectrum' in f[last_iteration]:
+                    quantities_path['specific_energy_spectrum'] = last_iteration
+
             # Minimum specific energy
             if use_minimum_specific_energy:
                 minimum_specific_energy_path = '/Input/Grid/Quantities'
@@ -313,6 +320,10 @@ def use_quantities(self, filename, quantities=['density', 'specific_energy'],
                 if 'specific_energy' in f['/Grid/Quantities']:
                     quantities_path['specific_energy'] = '/Grid/Quantities'
 
+            if 'specific_energy_spectrum' in quantities:
+                if 'specific_energy_spectrum' in f['/Grid/Quantities']:
+                    quantities_path['specific_energy_spectrum'] = '/Grid/Quantities'
+
             # Minimum specific energy
             if use_minimum_specific_energy:
                 minimum_specific_energy_path = '/Grid/Quantities'