diff --git a/tardis/io/model_reader.py b/tardis/io/model_reader.py index 540848220c3..4c95ae842b6 100644 --- a/tardis/io/model_reader.py +++ b/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 diff --git a/tardis/io/schemas/model.yml b/tardis/io/schemas/model.yml index e317d21b68b..c6e76ed0b65 100644 --- a/tardis/io/schemas/model.yml +++ b/tardis/io/schemas/model.yml @@ -104,6 +104,7 @@ definitions: enum: - simple_ascii - artis + - tardis_model description: file type v_inner_boundary: type: quantity diff --git a/tardis/io/tests/data/tardis_configv1_tardis_model_format.yml b/tardis/io/tests/data/tardis_configv1_tardis_model_format.yml new file mode 100755 index 00000000000..b3b6a54f6c3 --- /dev/null +++ b/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 diff --git a/tardis/io/tests/data/tardis_model_abund.csv b/tardis/io/tests/data/tardis_model_abund.csv deleted file mode 100644 index 46200ee5b77..00000000000 --- a/tardis/io/tests/data/tardis_model_abund.csv +++ /dev/null @@ -1,11 +0,0 @@ -C O Mg Si Ni56 Ni58 -0 0 0 0.6 0.4 0 -0 0 0 0.1 0.5 0.4 -0 0 0 0.3 0 0.7 -0 0.2 0.8 0 0 0 -0 0.3 0.7 0 0 0 -0 0.2 0.8 0 0 0 -0 0.2 0.8 0 0 0 -0 0.2 0.8 0 0 0 -0.5 0.5 0 0 0 0 -0.5 0.5 0 0 0 0 diff --git a/tardis/io/tests/data/tardis_model_format.csv b/tardis/io/tests/data/tardis_model_format.csv new file mode 100644 index 00000000000..4e7dffee474 --- /dev/null +++ b/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 diff --git a/tardis/io/tests/test_model_reader.py b/tardis/io/tests/test_model_reader.py index b44aa60494c..8bb7f08517e 100644 --- a/tardis/io/tests/test_model_reader.py +++ b/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 diff --git a/tardis/model/base.py b/tardis/model/base.py index d57338c192f..a75e17cc9e7 100644 --- a/tardis/model/base.py +++ b/tardis/model/base.py @@ -62,7 +62,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 @@ -72,6 +72,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 @@ -285,6 +286,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, @@ -297,7 +300,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) @@ -307,7 +310,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 @@ -357,4 +362,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) diff --git a/tardis/plasma/properties/ion_population.py b/tardis/plasma/properties/ion_population.py index e7fb8f530d7..a50831ffcb7 100644 --- a/tardis/plasma/properties/ion_population.py +++ b/tardis/plasma/properties/ion_population.py @@ -184,10 +184,11 @@ class IonNumberDensity(ProcessingPlasmaProperty): outputs = ('ion_number_density', 'electron_densities') latex_name = ('N_{i,j}','n_{e}',) - def __init__(self, plasma_parent, ion_zero_threshold=1e-20): + def __init__(self, plasma_parent, ion_zero_threshold=1e-20, electron_densities=None): super(IonNumberDensity, self).__init__(plasma_parent) self.ion_zero_threshold = ion_zero_threshold self.block_ids = None + self._electron_densities = electron_densities @staticmethod def calculate_with_n_electron(phi, partition_function, @@ -223,31 +224,39 @@ def _calculate_block_ids(phi): return calculate_block_ids_from_dataframe(phi) def calculate(self, phi, partition_function, number_density): - n_e_convergence_threshold = 0.05 - n_electron = number_density.sum(axis=0) - n_electron_iterations = 0 - - while True: + if self._electron_densities is None: + n_e_convergence_threshold = 0.05 + n_electron = number_density.sum(axis=0) + n_electron_iterations = 0 + + while True: + ion_number_density, self.block_ids = \ + self.calculate_with_n_electron( + phi, partition_function, number_density, n_electron, + self.block_ids, self.ion_zero_threshold) + ion_numbers = ion_number_density.index.get_level_values(1).values + ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) + new_n_electron = (ion_number_density.values * ion_numbers).sum( + axis=0) + if np.any(np.isnan(new_n_electron)): + raise PlasmaIonizationError('n_electron just turned "nan" -' + ' aborting') + n_electron_iterations += 1 + if n_electron_iterations > 100: + logger.warn('n_electron iterations above 100 ({0}) -' + ' something is probably wrong'.format( + n_electron_iterations)) + if np.all(np.abs(new_n_electron - n_electron) + / n_electron < n_e_convergence_threshold): + break + n_electron = 0.5 * (new_n_electron + n_electron) + else: + n_electron = self._electron_densities ion_number_density, self.block_ids = \ self.calculate_with_n_electron( - phi, partition_function, number_density, n_electron, + phi, partition_function, number_density, n_electron, self.block_ids, self.ion_zero_threshold) - ion_numbers = ion_number_density.index.get_level_values(1).values - ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) - new_n_electron = (ion_number_density.values * ion_numbers).sum( - axis=0) - if np.any(np.isnan(new_n_electron)): - raise PlasmaIonizationError('n_electron just turned "nan" -' - ' aborting') - n_electron_iterations += 1 - if n_electron_iterations > 100: - logger.warn('n_electron iterations above 100 ({0}) -' - ' something is probably wrong'.format( - n_electron_iterations)) - if np.all(np.abs(new_n_electron - n_electron) - / n_electron < n_e_convergence_threshold): - break - n_electron = 0.5 * (new_n_electron + n_electron) + return ion_number_density, n_electron class IonNumberDensityHeNLTE(ProcessingPlasmaProperty): @@ -270,10 +279,11 @@ class IonNumberDensityHeNLTE(ProcessingPlasmaProperty): 'helium_population_updated') latex_name = ('N_{i,j}','n_{e}',) - def __init__(self, plasma_parent, ion_zero_threshold=1e-20): + def __init__(self, plasma_parent, ion_zero_threshold=1e-20, electron_densities=None): super(IonNumberDensityHeNLTE, self).__init__(plasma_parent) self.ion_zero_threshold = ion_zero_threshold self.block_ids = None + self._electron_densities = electron_densities def update_he_population(self, helium_population, n_electron, number_density): @@ -291,36 +301,44 @@ def update_he_population(self, helium_population, n_electron, def calculate(self, phi, partition_function, number_density, helium_population): - n_e_convergence_threshold = 0.05 - n_electron = number_density.sum(axis=0) - n_electron_iterations = 0 - while True: + if self._electron_densities is None: + n_e_convergence_threshold = 0.05 + n_electron = number_density.sum(axis=0) + n_electron_iterations = 0 + while True: + ion_number_density, self.block_ids = \ + IonNumberDensity.calculate_with_n_electron( + phi, partition_function, number_density, n_electron, + self.block_ids, self.ion_zero_threshold) + helium_population_updated = self.update_he_population( + helium_population, n_electron, number_density) + ion_number_density.ix[2].ix[0].update(helium_population_updated.ix[ + 0].sum(axis=0)) + ion_number_density.ix[2].ix[1].update(helium_population_updated.ix[ + 1].sum(axis=0)) + ion_number_density.ix[2].ix[2].update(helium_population_updated.ix[ + 2].ix[0]) + ion_numbers = ion_number_density.index.get_level_values(1).values + ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) + new_n_electron = (ion_number_density.values * ion_numbers).sum( + axis=0) + if np.any(np.isnan(new_n_electron)): + raise PlasmaIonizationError('n_electron just turned "nan" -' + ' aborting') + n_electron_iterations += 1 + if n_electron_iterations > 100: + logger.warn('n_electron iterations above 100 ({0}) -' + ' something is probably wrong'.format( + n_electron_iterations)) + if np.all(np.abs(new_n_electron - n_electron) + / n_electron < n_e_convergence_threshold): + break + n_electron = 0.5 * (new_n_electron + n_electron) + else: + n_electron = self._electron_densities ion_number_density, self.block_ids = \ IonNumberDensity.calculate_with_n_electron( - phi, partition_function, number_density, n_electron, + phi, partition_function, number_density, n_electron, self.block_ids, self.ion_zero_threshold) - helium_population_updated = self.update_he_population( - helium_population, n_electron, number_density) - ion_number_density.ix[2].ix[0].update(helium_population_updated.ix[ - 0].sum(axis=0)) - ion_number_density.ix[2].ix[1].update(helium_population_updated.ix[ - 1].sum(axis=0)) - ion_number_density.ix[2].ix[2].update(helium_population_updated.ix[ - 2].ix[0]) - ion_numbers = ion_number_density.index.get_level_values(1).values - ion_numbers = ion_numbers.reshape((ion_numbers.shape[0], 1)) - new_n_electron = (ion_number_density.values * ion_numbers).sum( - axis=0) - if np.any(np.isnan(new_n_electron)): - raise PlasmaIonizationError('n_electron just turned "nan" -' - ' aborting') - n_electron_iterations += 1 - if n_electron_iterations > 100: - logger.warn('n_electron iterations above 100 ({0}) -' - ' something is probably wrong'.format( - n_electron_iterations)) - if np.all(np.abs(new_n_electron - n_electron) - / n_electron < n_e_convergence_threshold): - break - n_electron = 0.5 * (new_n_electron + n_electron) - return ion_number_density, n_electron, helium_population_updated \ No newline at end of file + + return ion_number_density, n_electron, helium_population_updated diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index d837878dfd1..bc2f6cad7ce 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -18,7 +18,7 @@ JBluesDiluteBlackBody, JBluesDetailed, RadiationFieldCorrection, StimulatedEmissionFactor, - HeliumNumericalNLTE) + HeliumNumericalNLTE, IonNumberDensity, IonNumberDensityHeNLTE) logger = logging.getLogger(__name__) @@ -147,6 +147,16 @@ def assemble_plasma(config, model, atom_data=None): heating_rate_data_file=config.plasma.heating_rate_data_file) else: plasma_modules += helium_lte_properties + + if model._electron_densities: + electron_densities = pd.Series(model._electron_densities.cgs.value) + if config.plasma.helium_treatment == 'numerical-nlte': + property_kwargs[IonNumberDensityHeNLTE] = dict( + electron_densities=electron_densities) + else: + property_kwargs[IonNumberDensity] = dict( + electron_densities=electron_densities) + kwargs['helium_treatment'] = config.plasma.helium_treatment plasma = BasePlasma(plasma_properties=plasma_modules, diff --git a/tardis/plasma/tests/test_tardis_model_density_config.py b/tardis/plasma/tests/test_tardis_model_density_config.py new file mode 100644 index 00000000000..22d0a1130b5 --- /dev/null +++ b/tardis/plasma/tests/test_tardis_model_density_config.py @@ -0,0 +1,34 @@ +import pytest +import os +from tardis.io.config_reader import Configuration +from tardis.model import Radial1DModel +from tardis.plasma.standard_plasmas import assemble_plasma +from numpy.testing import assert_almost_equal + +data_path = os.path.join('tardis', 'io', 'tests', 'data') + + +@pytest.fixture +def tardis_model_density_config(): + filename = 'tardis_configv1_tardis_model_format.yml' + return Configuration.from_yaml(os.path.join(data_path, filename)) + + +@pytest.fixture() +def raw_model(tardis_model_density_config): + return Radial1DModel.from_config(tardis_model_density_config) + + +@pytest.fixture() +def raw_plasma(tardis_model_density_config, raw_model, kurucz_atomic_data): + return assemble_plasma(tardis_model_density_config, raw_model, kurucz_atomic_data) + + +def test_electron_densities(raw_plasma): + assert_almost_equal(raw_plasma.electron_densities[8], 2.72e+14) + assert_almost_equal(raw_plasma.electron_densities[3], 2.6e+14) + + +def test_t_rad(raw_plasma): + assert_almost_equal(raw_plasma.t_rad[5], 76459.592) + assert_almost_equal(raw_plasma.t_rad[3], 76399.042)