Merge pull request #772 from vg3095/cmf_density

[TEP007] CMFGEN density parser

tardis/io/model_reader.py

@@ -42,12 +42,21 @@ def read_density_file(filename, filetype):
 
     """
     file_parsers = {'artis': read_artis_density,
-                    'simple_ascii': read_simple_ascii_density}
+                    'simple_ascii': read_simple_ascii_density,
+                    'tardis_model': read_cmfgen_density}
+
+    electron_densities = None
+    temperature = None
+    if filetype == 'tardis_model':
+        (time_of_model, velocity,
+         unscaled_mean_densities, electron_densities, temperature) = read_cmfgen_density(filename)
+    else:
+        (time_of_model, velocity,
+         unscaled_mean_densities) = file_parsers[filetype](filename)
 
-    (time_of_model, velocity,
-     unscaled_mean_densities) = file_parsers[filetype](filename)
     v_inner = velocity[:-1]
     v_outer = velocity[1:]
+
     invalid_volume_mask = (v_outer - v_inner) <= 0
     if invalid_volume_mask.sum() > 0:
         message = "\n".join(["cell {0:d}: v_inner {1:s}, v_outer "
@@ -59,7 +68,7 @@ def read_density_file(filename, filetype):
         raise ConfigurationError("Invalid volume of following cell(s):\n"
                                  "{:s}".format(message))
 
-    return time_of_model, velocity, unscaled_mean_densities
+    return time_of_model, velocity, unscaled_mean_densities, electron_densities, temperature
 
 def read_abundances_file(abundance_filename, abundance_filetype,
                          inner_boundary_index=None, outer_boundary_index=None):
@@ -235,6 +244,54 @@ def read_artis_density(fname):
     return time_of_model, velocity, mean_density
 
 
+def read_cmfgen_density(fname):
+    """
+    Reading a density file of the following structure (example; lines starting with a hash will be ignored):
+    The first density describes the mean density in the center of the model and is not used.
+    The file consists of a header row and next row contains unit of the respective attributes
+    velocity densities electron_densities temperature
+    km/s g/cm^3 /cm^3 K
+    871.66905 4.2537191e-09 2.5953807e+14 7.6395577
+    877.44269 4.2537191e-09 2.5953807e+14 7.6395577
+
+    Rest columns contain abundances of elements and isotopes
+
+    Parameters
+    ----------
+
+    fname: str
+        filename or path with filename
+
+
+    Returns
+    -------
+
+    time_of_model: ~astropy.units.Quantity
+        time at which the model is valid
+
+    velocity: ~np.ndarray
+    mean_density: ~np.ndarray
+    electron_densities: ~np.ndarray
+    temperature: ~np.ndarray
+
+    """
+    df = pd.read_csv(fname, comment='#', delimiter='\s+', skiprows=[0, 2])
+
+    with open(fname) as fh:
+        for row_index, line in enumerate(fh):
+            if row_index == 0:
+                time_of_model_string = line.strip().replace('t0:', '')
+                time_of_model = parse_quantity(time_of_model_string)
+            elif row_index == 2:
+                quantities = line.split()
+
+    velocity = u.Quantity(df['velocity'].values, quantities[0]).to('cm/s')
+    temperature = u.Quantity(df['temperature'].values, quantities[1])[1:]
+    mean_density = u.Quantity(df['densities'].values, quantities[2])[1:]
+    electron_densities = u.Quantity(
+        df['electron_densities'].values, quantities[3])[1:]
+
+    return time_of_model, velocity, mean_density, electron_densities, temperature
 
 def read_simple_ascii_abundances(fname):
     """
@@ -267,28 +324,55 @@ def read_simple_ascii_abundances(fname):
 
 
 def read_simple_isotope_abundances(fname, delimiter='\s+'):
-    df = pd.read_csv(fname, comment='#', delimiter=delimiter)
+    """
+    Reading an abundance file of the following structure (example; lines starting with hash will be ignored):
+    The first line of abundances describe the abundances in the center of the model and are not used.
+    First 4 columns contain values related to velocity, density, electron_density and temperature.
+    From 5th column onwards, abundances of elements and isotopes begin.
+    The file consists of a header row and next row contains unit of the respective attributes
+    Since abundance fractions are unitless , its unit row is filled with ones
+    Example 
+    velocity...temperature C O Ni56
+    km/s.........K         1 1 1
+    ...................... 0.4 0.3 0.2
+
+    Parameters
+    ----------
+
+    fname: str
+        filename or path with filename
+
+    Returns
+    -------
+
+    index: ~np.ndarray
+    abundances: ~pandas.DataFrame
+    isotope_abundance: ~pandas.MultiIndex    
+    """
+    df = pd.read_csv(fname, comment='#',
+                     delimiter=delimiter, skiprows=[0, 2])
     df = df.transpose()
 
-    abundance = pd.DataFrame(columns=np.arange(df.shape[1]),
+    abundance = pd.DataFrame(columns=np.arange(df.shape[1] - 1),
                              index=pd.Index([],
                                             name='atomic_number'),
                              dtype=np.float64)
 
     isotope_index = pd.MultiIndex(
         [[]] * 2, [[]] * 2, names=['atomic_number', 'mass_number'])
-    isotope_abundance = pd.DataFrame(columns=np.arange(df.shape[1]),
+    isotope_abundance = pd.DataFrame(columns=np.arange(df.shape[1] - 1),
                                      index=isotope_index,
                                      dtype=np.float64)
-
-    for element_symbol_string in df.index:
+
+    #First 4 columns related to density parser (e.g. velocity)
+    for element_symbol_string in df.index[4:]:
         if element_symbol_string in nucname.name_zz:
             z = nucname.name_zz[element_symbol_string]
-            abundance.loc[z, :] = df.loc[element_symbol_string].tolist()
+            abundance.loc[z, :] = df.loc[element_symbol_string].tolist()[1:]
         else:
             z = nucname.znum(element_symbol_string)
             mass_no = nucname.anum(element_symbol_string)
             isotope_abundance.loc[(
-                z, mass_no), :] = df.loc[element_symbol_string].tolist()
+                z, mass_no), :] = df.loc[element_symbol_string].tolist()[1:]
 
     return abundance.index, abundance, isotope_abundance

tardis/io/schemas/model.yml

@@ -104,6 +104,7 @@ definitions:
           enum:
           - simple_ascii
           - artis
+          - tardis_model
           description: file type
         v_inner_boundary:
           type: quantity

tardis/io/tests/data/tardis_configv1_tardis_model_format.yml

@@ -0,0 +1,36 @@
+tardis_config_version: v1.0
+
+supernova:
+    luminosity_requested: 2.8e9 solLum
+    time_explosion: 13 day
+
+atom_data: kurucz_cd23_chianti_H_He.h5
+model:
+
+    structure:
+        type: file
+        filename: tardis_model_format.csv
+        filetype: tardis_model
+
+    abundances:
+        type: file
+        filename: tardis_model_format.csv
+        filetype: tardis_model
+
+plasma:
+    ionization: lte
+    excitation: lte
+    radiative_rates_type: dilute-blackbody
+    line_interaction_type: macroatom
+
+montecarlo:
+    seed: 23111963
+    no_of_packets : 2.0e+5
+    iterations: 5
+    last_no_of_packets: 5.0e+5
+    no_of_virtual_packets: 5
+
+spectrum:
+    start: 500 angstrom
+    stop: 20000 angstrom
+    num: 10000

tardis/io/tests/data/tardis_model_format.csv

@@ -0,0 +1,13 @@
+t0: 0.976 day       
+velocity temperature densities electron_densities C O Mg Si Ni56 Ni58
+km/s K g/cm^3 /cm^3 1 1 1 1 1 1
+871.66905 76395.577 4.2537191E-09 2.60E+14 0 0 0 0.6 0.4 0
+877.44269 76395.577 4.2537191E-09 2.60E+14 0 0 0 0.1 0.5 0.4
+894.99407 76395.631 4.2537191E-09 2.60E+14 0 0 0 0.3 0 0.7
+931.1571 76396.057 4.2537191E-09 2.60E+14 0 0.2 0.8 0 0 0
+990.30752 76399.042 4.2537271E-09 2.60E+14 0 0.3 0.7 0 0 0
+1050.8676 76411.983 4.2539537E-09 2.60E+14 0 0.2 0.8 0 0 0
+1115.1545 76459.592 4.2563604E-09 2.60E+14 0 0.2 0.8 0 0 0
+1183.3741 76633.367 4.2683079E-09 2.61E+14 0 0.2 0.8 0 0 0
+1255.767 77312.12 4.290997E-09 2.64E+14 0.5 0.5 0 0 0 0
+1332.5886 79602.375 4.3396835E-09 2.72E+14 0.5 0.5 0 0 0 0

tardis/io/tests/test_model_reader.py

@@ -7,7 +7,7 @@
 import tardis
 from tardis.io.config_reader import Configuration
 from tardis.io.model_reader import (
-    read_artis_density, read_simple_ascii_abundances, read_simple_isotope_abundances, read_uniform_abundances)
+    read_artis_density, read_simple_ascii_abundances, read_simple_isotope_abundances, read_uniform_abundances, read_cmfgen_density)
 
 data_path = os.path.join(tardis.__path__[0], 'io', 'tests', 'data')
 
@@ -21,8 +21,8 @@ def artis_abundances_fname():
 
 
 @pytest.fixture
-def tardis_model_abundance_fname():
-    return os.path.join(data_path, 'tardis_model_abund.csv')
+def tardis_model_fname():
+    return os.path.join(data_path, 'tardis_model_format.csv')
 
 
 @pytest.fixture
@@ -48,14 +48,14 @@ def test_read_simple_ascii_abundances(artis_abundances_fname):
     assert np.isclose(abundances[23].ix[2], 2.672351e-08 , atol=1.e-12)
 
 
-def test_read_simple_isotope_abundances(tardis_model_abundance_fname):
+def test_read_simple_isotope_abundances(tardis_model_fname):
     index, abundances, isotope_abundance = read_simple_isotope_abundances(
-        tardis_model_abundance_fname)
-    assert np.isclose(abundances.loc[6, 9], 0.5, atol=1.e-12)
-    assert np.isclose(abundances.loc[12, 6], 0.8, atol=1.e-12)
-    assert np.isclose(abundances.loc[14, 2], 0.3, atol=1.e-12)
-    assert np.isclose(isotope_abundance.loc[(28, 56), 1], 0.5, atol=1.e-12)
-    assert np.isclose(isotope_abundance.loc[(28, 58), 2], 0.7, atol=1.e-12)
+        tardis_model_fname)
+    assert np.isclose(abundances.loc[6, 8], 0.5, atol=1.e-12)
+    assert np.isclose(abundances.loc[12, 5], 0.8, atol=1.e-12)
+    assert np.isclose(abundances.loc[14, 1], 0.3, atol=1.e-12)
+    assert np.isclose(isotope_abundance.loc[(28, 56), 0], 0.5, atol=1.e-12)
+    assert np.isclose(isotope_abundance.loc[(28, 58), 1], 0.7, atol=1.e-12)
 
 
 def test_read_uniform_abundances(isotope_uniform_abundance):
@@ -65,3 +65,16 @@ def test_read_uniform_abundances(isotope_uniform_abundance):
     assert np.isclose(abundances.loc[20, 5], 0.03, atol=1.e-12)
     assert np.isclose(isotope_abundance.loc[(28, 56), 15], 0.05, atol=1.e-12)
     assert np.isclose(isotope_abundance.loc[(28, 58), 2], 0.05, atol=1.e-12)
+
+
+def test_simple_read_cmfgen_density(tardis_model_fname):
+    time_of_model, velocity, mean_density, electron_densities, temperature = read_cmfgen_density(
+        tardis_model_fname)
+
+    assert np.isclose(0.976 * u.day, time_of_model, atol=1e-7 * u.day)
+    assert np.isclose(mean_density[4], 4.2539537e-09 * u.g / u.cm**3, atol=1.e-6
+                      * u.g / u.cm**3)
+    assert np.isclose(electron_densities[5], 2.6e+14 * u.cm**-3, atol=1.e-6
+                      * u.cm**-3)
+    assert len(mean_density) == 9
+    assert len(velocity) == len(mean_density) + 1

tardis/model/base.py

@@ -63,7 +63,7 @@ class Radial1DModel(HDFWriterMixin):
 
     def __init__(self, velocity, homologous_density, abundance, isotope_abundance,
                  time_explosion, t_inner, luminosity_requested=None, t_radiative=None,
-                 dilution_factor=None, v_boundary_inner=None, v_boundary_outer=None):
+                 dilution_factor=None, v_boundary_inner=None, v_boundary_outer=None, electron_densities=None):
         self._v_boundary_inner = None
         self._v_boundary_outer = None
         self._velocity = None
@@ -73,6 +73,7 @@ def __init__(self, velocity, homologous_density, abundance, isotope_abundance,
         self.homologous_density = homologous_density
         self._abundance = abundance
         self.time_explosion = time_explosion
+        self._electron_densities = electron_densities
 
         self.raw_abundance = self._abundance
         self.raw_isotope_abundance = isotope_abundance 
@@ -289,6 +290,8 @@ def from_config(cls, config):
         time_explosion = config.supernova.time_explosion.cgs
 
         structure = config.model.structure
+        electron_densities = None
+        temperature = None
         if structure.type == 'specific':
             velocity = quantity_linspace(structure.velocity.start,
                                          structure.velocity.stop,
@@ -301,7 +304,7 @@ def from_config(cls, config):
                 structure_fname = os.path.join(config.config_dirname,
                                                structure.filename)
 
-            time_0, velocity, density_0 = read_density_file(
+            time_0, velocity, density_0, electron_densities, temperature = read_density_file(
                 structure_fname, structure.filetype)
             density_0 = density_0.insert(0, 0)
             homologous_density = HomologousDensity(density_0, time_0)
@@ -311,7 +314,9 @@ def from_config(cls, config):
         #       v boundaries.
         no_of_shells = len(velocity) - 1
 
-        if config.plasma.initial_t_rad > 0 * u.K:
+        if temperature:
+            t_radiative = temperature
+        elif config.plasma.initial_t_rad > 0 * u.K:
             t_radiative = np.ones(no_of_shells) * config.plasma.initial_t_rad
         else:
             t_radiative = None
@@ -363,4 +368,5 @@ def from_config(cls, config):
                    luminosity_requested=luminosity_requested,
                    dilution_factor=None,
                    v_boundary_inner=structure.get('v_inner_boundary', None),
-                   v_boundary_outer=structure.get('v_outer_boundary', None))
+                   v_boundary_outer=structure.get('v_outer_boundary', None),
+                   electron_densities=electron_densities)