diff --git a/tardis/io/config_reader.py b/tardis/io/config_reader.py index 6fc603d5e9f..77a72d6c32d 100644 --- a/tardis/io/config_reader.py +++ b/tardis/io/config_reader.py @@ -460,11 +460,6 @@ def parse_supernova_section(supernova_dict): config_dict['time_explosion'] = parse_quantity(supernova_dict['time_explosion']).to('s') - if 'distance' in supernova_dict: - config_dict['distance'] = parse_quantity(supernova_dict['distance']) - else: - config_dict['distance'] = None - if 'luminosity_wavelength_start' in supernova_dict: config_dict['luminosity_nu_end'] = parse_quantity(supernova_dict['luminosity_wavelength_start']). \ to('Hz', u.spectral()) @@ -480,35 +475,6 @@ def parse_supernova_section(supernova_dict): return config_dict -def parse_spectrum_list2dict(spectrum_list): - """ - Parse the spectrum list [start, stop, num] to a list - """ - if 'start' in spectrum_list and 'stop' in spectrum_list \ - and 'num' in spectrum_list: - spectrum_list = [spectrum_list['start'], spectrum_list['stop'], - spectrum_list['num']] - if spectrum_list[0].unit.physical_type != 'length' and \ - spectrum_list[1].unit.physical_type != 'length': - raise ValueError('start and end of spectrum need to be a length') - - - spectrum_config_dict = {} - spectrum_config_dict['start'] = spectrum_list[0] - spectrum_config_dict['end'] = spectrum_list[1] - spectrum_config_dict['bins'] = spectrum_list[2] - - spectrum_frequency = quantity_linspace( - spectrum_config_dict['end'].to('Hz', u.spectral()), - spectrum_config_dict['start'].to('Hz', u.spectral()), - num=spectrum_config_dict['bins'] + 1) - - spectrum_config_dict['frequency'] = spectrum_frequency - - return spectrum_config_dict - - - def parse_convergence_section(convergence_section_dict): """ Parse the convergence section dictionary @@ -943,22 +909,6 @@ def from_config_dict(cls, config_dict, atom_data=None, test_parser=False, else: plasma_section['t_rads'] = None - if plasma_section['disable_electron_scattering'] is False: - logger.debug("Electron scattering switched on") - validated_config_dict['montecarlo']['sigma_thomson'] = 6.652486e-25 / (u.cm ** 2) - else: - logger.warn('Disabling electron scattering - this is not physical') - validated_config_dict['montecarlo']['sigma_thomson'] = 1e-200 / (u.cm ** 2) - - if plasma_section['helium_treatment'] == 'recomb-nlte': - validated_config_dict['plasma']['helium_treatment'] == 'recomb-nlte' - else: - validated_config_dict['plasma']['helium_treatment'] == 'dilute-lte' - - - - - ##### NLTE subsection of Plasma start nlte_validated_config_dict = {} nlte_species = [] @@ -1015,21 +965,9 @@ def from_config_dict(cls, config_dict, atom_data=None, test_parser=False, black_body_section['stop'] montecarlo_section['black_body_sampling']['samples'] = \ black_body_section['num'] - virtual_spectrum_section = montecarlo_section['virtual_spectrum_range'] - montecarlo_section['virtual_spectrum_range'] = {} - montecarlo_section['virtual_spectrum_range']['start'] = \ - virtual_spectrum_section['start'] - montecarlo_section['virtual_spectrum_range']['end'] = \ - virtual_spectrum_section['stop'] - montecarlo_section['virtual_spectrum_range']['samples'] = \ - virtual_spectrum_section['num'] ###### END of convergence section reading - - validated_config_dict['spectrum'] = parse_spectrum_list2dict( - validated_config_dict['spectrum']) - return cls(validated_config_dict, atom_data) diff --git a/tardis/io/tests/test_config_reader.py b/tardis/io/tests/test_config_reader.py index 073f8b39a68..56aafe9ddbf 100644 --- a/tardis/io/tests/test_config_reader.py +++ b/tardis/io/tests/test_config_reader.py @@ -31,14 +31,6 @@ def test_quantity_linspace(): assert len(quantity_linspace) == 1000 -def test_spectrum_list2_dict(): - spectrum_dict = config_reader.parse_spectrum_list2dict( - [200*u.angstrom, 10000 * u.angstrom, 100]) - assert_almost_equal(spectrum_dict['start'].to(u.angstrom).value, 200) - assert_almost_equal(spectrum_dict['end'].to(u.angstrom).value, 10000) - assert_almost_equal(spectrum_dict['bins'], 100) - - def test_convergence_section_parser(): test_convergence_section = {'type': 'damped', 'lock_t_inner_cyles': 1, diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index 9ce42a283e6..c870834672c 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -7,13 +7,14 @@ from scipy.special import zeta from spectrum import TARDISSpectrum +from tardis.util import quantity_linspace from tardis.montecarlo import montecarlo, packet_source from tardis.io.util import to_hdf import numpy as np import pandas as pd -logger = logging.getLevelName(__name__) +logger = logging.getLogger(__name__) class MontecarloRunner(object): """ @@ -28,9 +29,18 @@ class MontecarloRunner(object): t_rad_estimator_constant = ((np.pi**4 / (15 * 24 * zeta(5, 1))) * (const.h / const.k_B)).cgs.value - def __init__(self, seed, spectrum_frequency, distance=None): + def __init__(self, seed, spectrum_frequency, virtual_spectrum_range, + sigma_thomson, enable_reflective_inner_boundary, + inner_boundary_albedo, line_interaction_type, distance=None): + + self.seed = seed self.packet_source = packet_source.BlackBodySimpleSource(seed) self.spectrum_frequency = spectrum_frequency + self.virtual_spectrum_range = virtual_spectrum_range # TODO: linspace handling + self.sigma_thomson = sigma_thomson + self.enable_reflective_inner_boundary = enable_reflective_inner_boundary + self.inner_boundary_albedo = inner_boundary_albedo + self.line_interaction_type = line_interaction_type self.spectrum = TARDISSpectrum(spectrum_frequency, distance) self.spectrum_virtual = TARDISSpectrum(spectrum_frequency, distance) self.spectrum_reabsorbed = TARDISSpectrum(spectrum_frequency, distance) @@ -53,18 +63,18 @@ def _initialize_estimator_arrays(self, no_of_shells, tau_sobolev_shape): self.Edotlu_estimator = np.zeros(tau_sobolev_shape) - def _initialize_geometry_arrays(self, structure): + def _initialize_geometry_arrays(self, model): """ Generate the cgs like geometry arrays for the montecarlo part Parameters ---------- - structure: ~ConfigurationNameSpace + model : model.Radial1DModel """ - self.r_inner_cgs = structure.r_inner.to('cm').value - self.r_outer_cgs = structure.r_outer.to('cm').value - self.v_inner_cgs = structure.v_inner.to('cm/s').value + self.r_inner_cgs = model.r_inner.to('cm').value + self.r_outer_cgs = model.r_outer.to('cm').value + self.v_inner_cgs = model.v_inner.to('cm/s').value def _initialize_packets(self, T, no_of_packets, no_of_virtual_packets=None): nus, mus, energies = self.packet_source.create_packets(T, no_of_packets) @@ -110,7 +120,7 @@ def legacy_update_spectrum(self, no_of_virtual_packets): self.spectrum_virtual.update_luminosity( self.montecarlo_virtual_luminosity) - def run(self, model, no_of_packets, no_of_virtual_packets=0, nthreads=1,last_run=False): + def run(self, model, plasma, no_of_packets, no_of_virtual_packets=0, nthreads=1,last_run=False): """ Running the TARDIS simulation @@ -118,21 +128,22 @@ def run(self, model, no_of_packets, no_of_virtual_packets=0, nthreads=1,last_run ---------- :param model: + :param plasma: :param no_of_virtual_packets: :param nthreads: :return: """ - self.time_of_simulation = model.time_of_simulation - self.volume = model.tardis_config.structure.volumes + self.time_of_simulation = self.calculate_time_of_simulation(model) + self.volume = model.volume self._initialize_estimator_arrays(self.volume.shape[0], - model.plasma.tau_sobolevs.shape) - self._initialize_geometry_arrays(model.tardis_config.structure) + plasma.tau_sobolevs.shape) + self._initialize_geometry_arrays(model) self._initialize_packets(model.t_inner.value, no_of_packets) montecarlo.montecarlo_radial1d( - model, self, virtual_packet_flag=no_of_virtual_packets, + model, plasma, self, virtual_packet_flag=no_of_virtual_packets, nthreads=nthreads,last_run=last_run) # Workaround so that j_blue_estimator is in the right ordering # They are written as an array of dimension (no_of_shells, no_of_lines) @@ -277,6 +288,13 @@ def calculate_radiationfield_properties(self): return t_rad * u.K, w + def calculate_luminosity_inner(self, model): + return (4 * np.pi * const.sigma_sb.cgs * + model.r_inner[0] ** 2 * model.t_inner ** 4).to('erg/s') + + def calculate_time_of_simulation(self, model): + return (1.0 * u.erg / self.calculate_luminosity_inner(model)) + def calculate_f_nu(self, frequency): pass @@ -308,3 +326,27 @@ def to_hdf(self, path_or_buf, path=''): 'spectrum_virtual') self.spectrum_reabsorbed.to_hdf(path_or_buf, runner_path, 'spectrum_reabsorbed') + + @classmethod + def from_config(cls, config): + if config.plasma.disable_electron_scattering: + logger.warn('Disabling electron scattering - this is not physical') + sigma_thomson = 1e-200 / (u.cm ** 2) + else: + logger.debug("Electron scattering switched on") + sigma_thomson = 6.652486e-25 / (u.cm ** 2) + + spectrum_frequency = quantity_linspace( + config.spectrum.stop.to('Hz', u.spectral()), + config.spectrum.start.to('Hz', u.spectral()), + num=config.spectrum.num + 1) + + return cls(seed=config.montecarlo.seed, + spectrum_frequency=spectrum_frequency, + # TODO: Linspace handling for virtual_spectrum_range + virtual_spectrum_range=config.montecarlo.virtual_spectrum_range, + sigma_thomson=sigma_thomson, + enable_reflective_inner_boundary=config.montecarlo.enable_reflective_inner_boundary, + inner_boundary_albedo=config.montecarlo.inner_boundary_albedo, + line_interaction_type=config.plasma.line_interaction_type, + distance=config.supernova.get('distance', None)) \ No newline at end of file diff --git a/tardis/montecarlo/montecarlo.pyx b/tardis/montecarlo/montecarlo.pyx index 8e412fde4b3..82f060b6af5 100644 --- a/tardis/montecarlo/montecarlo.pyx +++ b/tardis/montecarlo/montecarlo.pyx @@ -93,7 +93,7 @@ cdef extern from "src/cmontecarlo.h": -cdef initialize_storage_model(model, runner, storage_model_t *storage): +cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage): """ Initializing the storage struct. @@ -105,8 +105,7 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage): storage.packet_energies = PyArray_DATA(runner.input_energy) # Setup of structure - structure = model.tardis_config.structure - storage.no_of_shells = structure.no_of_shells + storage.no_of_shells = model.no_of_shells storage.r_inner = PyArray_DATA(runner.r_inner_cgs) @@ -115,15 +114,14 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage): # Setup the rest # times - storage.time_explosion = model.tardis_config.supernova.time_explosion.to( - 's').value + storage.time_explosion = model.time_explosion.to('s').value storage.inverse_time_explosion = 1.0 / storage.time_explosion #electron density storage.electron_densities = PyArray_DATA( - model.plasma.electron_densities.values) + plasma.electron_densities.values) runner.inverse_electron_densities = ( - 1.0 / model.plasma.electron_densities.values) + 1.0 / plasma.electron_densities.values) storage.inverse_electron_densities = PyArray_DATA( runner.inverse_electron_densities) # Switch for continuum processes @@ -146,10 +144,10 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage): storage.l_pop_r = l_pop_r.data # Line lists - storage.no_of_lines = model.atom_data.lines.nu.values.size - storage.line_list_nu = PyArray_DATA(model.atom_data.lines.nu.values) + storage.no_of_lines = plasma.atomic_data.lines.nu.values.size + storage.line_list_nu = PyArray_DATA(plasma.atomic_data.lines.nu.values) runner.line_lists_tau_sobolevs = ( - model.plasma.tau_sobolevs.values.flatten(order='F') + plasma.tau_sobolevs.values.flatten(order='F') ) storage.line_lists_tau_sobolevs = PyArray_DATA( runner.line_lists_tau_sobolevs @@ -161,30 +159,30 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage): runner.Edotlu_estimator) storage.line_interaction_id = runner.get_line_interaction_id( - model.tardis_config.plasma.line_interaction_type) + runner.line_interaction_type) # macro atom & downbranch if storage.line_interaction_id >= 1: runner.transition_probabilities = ( - model.plasma.transition_probabilities.values.flatten(order='F') + plasma.transition_probabilities.values.flatten(order='F') ) storage.transition_probabilities = PyArray_DATA( runner.transition_probabilities ) storage.transition_probabilities_nd = ( - model.plasma.transition_probabilities.values.shape[0]) + plasma.transition_probabilities.values.shape[0]) storage.line2macro_level_upper = PyArray_DATA( - model.atom_data.lines_upper2macro_reference_idx) + plasma.atomic_data.lines_upper2macro_reference_idx) storage.macro_block_references = PyArray_DATA( - model.atom_data.macro_atom_references['block_references'].values) + plasma.atomic_data.macro_atom_references['block_references'].values) storage.transition_type = PyArray_DATA( - model.atom_data.macro_atom_data['transition_type'].values) + plasma.atomic_data.macro_atom_data['transition_type'].values) # Destination level is not needed and/or generated for downbranch storage.destination_level_id = PyArray_DATA( - model.atom_data.macro_atom_data['destination_level_idx'].values) + plasma.atomic_data.macro_atom_data['destination_level_idx'].values) storage.transition_line_id = PyArray_DATA( - model.atom_data.macro_atom_data['lines_idx'].values) + plasma.atomic_data.macro_atom_data['lines_idx'].values) storage.output_nus = PyArray_DATA(runner._output_nu) storage.output_energies = PyArray_DATA(runner._output_energy) @@ -203,25 +201,25 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage): storage.js = PyArray_DATA(runner.j_estimator) storage.nubars = PyArray_DATA(runner.nu_bar_estimator) - storage.spectrum_start_nu = model.tardis_config.spectrum.frequency.value.min() - storage.spectrum_end_nu = model.tardis_config.spectrum.frequency.value.max() - storage.spectrum_virt_start_nu = model.tardis_config.montecarlo.virtual_spectrum_range.end.to('Hz', units.spectral()).value - storage.spectrum_virt_end_nu = model.tardis_config.montecarlo.virtual_spectrum_range.start.to('Hz', units.spectral()).value - storage.spectrum_delta_nu = model.tardis_config.spectrum.frequency.value[1] - model.tardis_config.spectrum.frequency.value[0] + storage.spectrum_start_nu = runner.spectrum_frequency.value.min() + storage.spectrum_end_nu = runner.spectrum_frequency.value.max() + # TODO: Linspace handling for virtual_spectrum_range + storage.spectrum_virt_start_nu = runner.virtual_spectrum_range.stop.to('Hz', units.spectral()).value + storage.spectrum_virt_end_nu = runner.virtual_spectrum_range.start.to('Hz', units.spectral()).value + storage.spectrum_delta_nu = runner.spectrum_frequency.value[1] - runner.spectrum_frequency.value[0] storage.spectrum_virt_nu = PyArray_DATA( runner.legacy_montecarlo_virtual_luminosity) - storage.sigma_thomson = ( - model.tardis_config.montecarlo.sigma_thomson.cgs.value) + storage.sigma_thomson = runner.sigma_thomson.cgs.value storage.inverse_sigma_thomson = 1.0 / storage.sigma_thomson - storage.reflective_inner_boundary = model.tardis_config.montecarlo.enable_reflective_inner_boundary - storage.inner_boundary_albedo = model.tardis_config.montecarlo.inner_boundary_albedo + storage.reflective_inner_boundary = runner.enable_reflective_inner_boundary + storage.inner_boundary_albedo = runner.inner_boundary_albedo # Data for continuum implementation - cdef np.ndarray[double, ndim=1] t_electrons = model.plasma.t_electrons + cdef np.ndarray[double, ndim=1] t_electrons = plasma.t_electrons storage.t_electrons = t_electrons.data -def montecarlo_radial1d(model, runner, int_type_t virtual_packet_flag=0, +def montecarlo_radial1d(model, plasma, runner, int_type_t virtual_packet_flag=0, int nthreads=4,last_run=False): """ Parameters @@ -253,10 +251,9 @@ def montecarlo_radial1d(model, runner, int_type_t virtual_packet_flag=0, cdef storage_model_t storage - initialize_storage_model(model, runner, &storage) + initialize_storage_model(model, plasma, runner, &storage) - montecarlo_main_loop(&storage, virtual_packet_flag, nthreads, - model.tardis_config.montecarlo.seed) + montecarlo_main_loop(&storage, virtual_packet_flag, nthreads, runner.seed) cdef np.ndarray[double, ndim=1] virt_packet_nus = np.zeros(storage.virt_packet_count, dtype=np.float64) cdef np.ndarray[double, ndim=1] virt_packet_energies = np.zeros(storage.virt_packet_count, dtype=np.float64) cdef np.ndarray[double, ndim=1] virt_packet_last_interaction_in_nu = np.zeros(storage.virt_packet_count, dtype=np.float64) @@ -298,8 +295,8 @@ def montecarlo_radial1d(model, runner, int_type_t virtual_packet_flag=0, def postprocess(model, runner): Edotlu_norm_factor = (1 / - (model.time_of_simulation * model.tardis_config.structure.volumes)) + (runner.time_of_simulation * model.volume)) exptau = 1 - np.exp(- runner.line_lists_tau_sobolevs.reshape(-1, runner.j_estimator.shape[0]) ) - model.Edotlu = Edotlu_norm_factor*exptau*runner.Edotlu_estimator + runner.Edotlu = Edotlu_norm_factor*exptau*runner.Edotlu_estimator