From cf7971104d2b06857093d3bf60aef317ec1bf196 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Sun, 3 Mar 2013 21:33:33 -0500 Subject: [PATCH 1/5] don't know what I added --- tardis/atomic.py | 35 ++++++++++++++++------------------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/tardis/atomic.py b/tardis/atomic.py index 6b983b2830b..faa3c4d5cb6 100644 --- a/tardis/atomic.py +++ b/tardis/atomic.py @@ -52,17 +52,17 @@ def read_hdf5_data(fname, dset_name): h5_file = h5py.File(fname) dataset = h5_file[dset_name] data = np.asarray(dataset) - data_units = dataset.attrs['units'] +# data_units = dataset.attrs['units'] data_table = table.Table(data) - for i, col_unit in enumerate(data_units): - if col_unit == 'n': - data_table.columns[i].units = None - elif col_unit == '1': - data_table.columns[i].units = units.Unit(1) - else: - data_table.columns[i].units = units.Unit(col_unit) +# for i, col_unit in enumerate(data_units): +# if col_unit == 'n': +# data_table.columns[i].units = None +# elif col_unit == '1': +# data_table.columns[i].units = units.Unit(1) +# else: +# data_table.columns[i].units = units.Unit(col_unit) h5_file.close() @@ -86,7 +86,7 @@ def read_basic_atom_data(fname=None): """ data_table = read_hdf5_data(fname, 'basic_atom_data') - data_table.columns['mass'].convert_units_to('g') +# data_table.columns['mass'] = units.Unit('u').to('g', data_table['mass']) return data_table @@ -110,7 +110,7 @@ def read_ionization_data(fname=None): """ data_table = read_hdf5_data(fname, 'ionization_data') - data_table.columns['ionization_energy'].convert_units_to('erg') + #data_table.columns['ionization_energy'] = units.Unit('eV').to('erg', data_table.columns['ionization_energy']) return data_table @@ -133,8 +133,9 @@ def read_levels_data(fname=None): """ data_table = read_hdf5_data(fname, 'levels_data') - data_table.columns['energy'].convert_units_to('erg') - #data_table.columns['ionization_energy'].convert_units_to('erg') + #data_table.columns['energy'].convert_units_to('erg') + #data_table.columns['energy'] = units.Unit('eV').to('erg', data_table.columns['energy']) + return data_table @@ -365,11 +366,14 @@ def __init__(self, atom_data, ionization_data, levels_data, lines_data, macro_at self.atom_data = DataFrame(atom_data.__array__()) self.atom_data.set_index('atomic_number', inplace=True) + self.atom_data.mass = units.Unit('u').to('g', self.atom_data.mass.values) self.ionization_data = DataFrame(ionization_data.__array__()) self.ionization_data.set_index(['atomic_number', 'ion_number'], inplace=True) + self.ionization_data.ionization_energy = units.Unit('eV').to('erg', self.ionization_data.ionization_energy.values) self.levels_data = DataFrame(levels_data.__array__()) + self.levels_data.energy = units.Unit('eV').to('erg', self.levels_data.energy.values) self.lines_data = DataFrame(lines_data.__array__()) self.lines_data.set_index('line_id', inplace=True) @@ -439,13 +443,6 @@ def prepare_atom_data(self, selected_atomic_numbers, line_interaction_type='scat self.atom_ion_index = None self.levels_index2atom_ion_index = None - einstein_coeff = (4 * np.pi ** 2 * constants.e.gauss.value ** 2) / ( - constants.m_e.cgs.value * constants.c.cgs.value) - - self.lines['B_lu'] = self.lines['f_lu'] * einstein_coeff / (constants.h.cgs.value * self.lines['nu']) - self.lines['B_ul'] = self.lines['f_ul'] * einstein_coeff / (constants.h.cgs.value * self.lines['nu']) - self.lines['A_ul'] = einstein_coeff * 2 * self.lines['nu'] ** 2 / constants.c.cgs.value ** 2 * self.lines[ - 'f_ul'] if self.has_macro_atom and not (line_interaction_type == 'scatter'): self.macro_atom_data = self.macro_atom_data_all[ From dcbe349d8fb29104dac719b1619ef9e2df27939a Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Tue, 5 Mar 2013 20:59:30 -0500 Subject: [PATCH 2/5] added a method simulate to radial1d_model. In addition, I switched to YAML for the configuration as this is much more flexible than our current config file. --- tardis/config_reader.py | 574 +++++++++++++++++++------------- tardis/model_radial_oned.py | 228 ++++++++----- tardis/montecarlo_multizone.pyx | 12 +- tardis/packet_source.py | 5 +- tardis/plasma.py | 13 +- tardis/simulation.py | 51 +-- 6 files changed, 514 insertions(+), 369 deletions(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index ab6bc6f9b28..4c74432d281 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -11,12 +11,185 @@ import re import pandas as pd from tardis import atomic +import yaml import tardis +import pdb + +import pprint logger = logging.getLogger(__name__) data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data')) +#File parsers for different file formats: + +density_structure_fileparser = {} + +class TardisConfigError(ValueError): + pass + + +def parse2quantity(quantity_string): + value_string, unit_string = quantity_string.split() + + value = float(value_string) + return units.Quantity(value, unit_string) + +def calculate_density_after_time(densities, time_0, time_explosion): + return densities * (time_explosion / time_0) ** -3 + + +def parse_density_file_section(density_file_dict, time_explosion): + + density_file_parser = {} + + def parse_artis_density(density_file_dict, time_explosion): + density_file = density_file_dict['name'] + for i, line in enumerate(file(density_file)): + if i == 0: + no_of_shells = int(line.strip()) + elif i == 1: + time_of_model = units.Quantity(float(line.strip()), 'day').to('s') + elif i == 2: + break + + velocities, mean_densities_0 = np.recfromtxt(density_file, skip_header=2, usecols=(1, 2), unpack=True) + #converting densities from log(g/cm^3) to g/cm^3 and stretching it to the current ti + velocities = units.Quantity(np.append([0], velocities), 'km/s').to('cm/s') + mean_densities_0 = units.Quantity(10**mean_densities_0, 'g/cm^3') + + mean_densities = calculate_density_after_time(time_of_model.value, mean_densities_0, + time_explosion) + + + #Verifying information + if len(mean_densities) == no_of_shells: + logger.debug('Verified ARTIS file %s (no_of_shells=length of dataset)', density_file) + else: + raise TardisConfigError('Error in ARTIS file %s - Number of shells not the same as dataset length' % density_file) + + min_shell = 0 + max_shell = no_of_shells + + v_inner = velocities[:-1] + v_outer = velocities[1:] + + volumes = (4 * np.pi /3) * (time_of_model**3) * ( v_outer**3 - v_inner**3) + masses = (volumes * mean_densities_0 / constants.M_sun).to(1) + + logger.info('Read ARTIS configuration file %s - found %d zones with total mass %g Msun', density_file, + no_of_shells, sum(masses.value)) + + + + if 'v_lowest' in density_file_dict: + v_lowest = parse2quantity(density_file_dict['v_lowest']).to('cm/s').value + min_shell = v_inner.value.searchsorted(v_lowest) + else: + min_shell = 0 + + if 'v_highest' in density_file_dict: + v_highest = parse2quantity(density_file_dict['v_highest']).to('cm/s').value + max_shell = v_outer.value.searchsorted(v_highest) + else: + max_shell = no_of_shells + + v_inner = v_inner[min_shell:max_shell] + v_outer = v_outer[min_shell:max_shell] + mean_densities = mean_densities[min_shell:max_shell] + + return v_inner.value, v_outer.value, mean_densities.value, min_shell, max_shell + + density_file_parser['artis'] = parse_artis_density + + try: + parser = density_file_parser[density_file_dict['type']] + except KeyError: + raise TardisConfigError('In abundance file section only types %s are allowed (supplied %s) ' % + (density_file_parser.keys(), density_file_dict['type'])) + + return parser(density_file_dict, time_explosion) + + + + + +def parse_velocity_section(velocity_dict, no_of_shells): + velocity_parser = {} + + def parse_linear_velocity(velocity_dict, no_of_shells): + v_inner = parse2quantity(velocity_dict['v_inner']).to('cm/s').value + v_outer = parse2quantity(velocity_dict['v_outer']).to('cm/s').value + velocities = np.linspace(v_inner, v_outer, no_of_shells + 1) + return velocities[:-1], velocities[1:] + + velocity_parser['linear'] = parse_linear_velocity + + try: + parser = velocity_parser[velocity_dict['type']] + except KeyError: + raise TardisConfigError('In velocity section only types %s are allowed (supplied %s) ' % + (velocity_parser.keys(), velocity_dict['type'])) + return parser(velocity_dict, no_of_shells) + +def parse_density_section(density_dict, no_of_shells, v_inner, v_outer, time_explosion): + density_parser = {} + + + #Parse density uniform + def parse_uniform(density_dict, no_of_shells, v_inner, v_outer, time_explosion): + return np.ones(no_of_shells) * parse2quantity(density_dict['value']).to('g cm^-3').value + density_parser['uniform'] = parse_uniform + + #Parse density branch85 w7 + def parse_branch85(density_dict, no_of_shells, v_inner, v_outer, time_explosion): + + time_0 = density_dict.pop('time_0', 19.9999584) + if isinstance(time_0, basestring): + time_0 = parse2quantity(time_0).to('s').value + else: + logger.debug('time_0 not supplied for density branch85 - using sensible default %g', time_0) + + density_coefficient = density_dict.pop('density_coefficient', None) + if density_coefficient is None: + density_coefficient = 3e29 + logger.debug('density_coefficient not supplied for density type branch85 - using sensible default %g', density_coefficient) + + velocities = 0.5 * (v_inner + v_outer) + densities = density_coefficient * (velocities * 1e-5) ** -7 + + densities = calculate_density_after_time(densities, time_0, time_explosion) + + return densities + density_parser['branch85_w7'] = parse_branch85 + + try: + parser = density_parser[density_dict['type']] + except KeyError: + raise TardisConfigError('In density section only types %s are allowed (supplied %s) ' % + (density_parser.keys(), density_dict['type'])) + return parser(density_dict, no_of_shells, v_inner, v_outer, time_explosion) + + +def parse_abundance_file_section(abundance_file_dict, abundances, min_shell, max_shell): + abundance_file_parser = {} + + def parse_artis(abundance_file_dict, abundances, min_shell, max_shell): + fname = abundance_file_dict['name'] + max_atom = 30 + logger.info("Parsing ARTIS Abundance section from shell %d to %d", min_shell, max_shell) + abundances.values[:,:max_atom] = np.loadtxt(fname)[min_shell:max_shell, 1:] + return abundances + abundance_file_parser['artis'] = parse_artis + + try: + parser = abundance_file_parser[abundance_file_dict['type']] + except KeyError: + raise TardisConfigError('In abundance file section only types %s are allowed (supplied %s) ' % + (abundance_file_parser.keys(), abundance_file_dict['type'])) + + return parser(abundance_file_dict, abundances, min_shell, max_shell) + def reformat_element_symbol(element_symbol): """ @@ -59,7 +232,9 @@ def calculate_w7_branch85_densities(velocities, time_explosion, time_0=19.999958 """ densities = density_coefficient * (velocities * 1e-5) ** -7 - densities *= (time_explosion / time_0) ** -3 + densities = calculate_density_after_time(densities, time_0, time_explosion) + + #densities *= (time_explosion / time_0) ** -3 return densities[1:] @@ -137,7 +312,8 @@ class TardisConfiguration(object): """ @classmethod - def from_file(cls, fname, args=None): + def from_ini(cls, fname, args=None): + print "This function is being deprecated and replaced by from_yaml classmethod" config_parse_object = ConfigParser() config_parse_object.read(fname) general_dict = dict(config_parse_object.items('general')) @@ -148,310 +324,246 @@ def from_file(cls, fname, args=None): config_object.parse_abundance_section(abundance_dict) return config_object + @classmethod + def from_yaml(cls, fname, args=None): + """ + Reading in from a YAML file and commandline args. Preferring commandline args when given - def __init__(self): - self.velocities = None - self.densities = None - self.abundances = None - self.initial_t_rad = 10000 - self.t_inner = 10000.0 # TODO calculate - self.ws = None - self.no_of_shells = None - - self._luminosity_outer = None - self.time_of_simulation = None + :param cls: + :param fname: + :param args: + :return: + """ + yaml_dict = yaml.load(file(fname)) + if yaml_dict['config_type'] not in ['simple1d']: + raise TardisConfigError('Only config_type=simple1d allowed at the moment.') - def parse_general_section(self, config_dict): - model_type = config_dict.pop('model_type') + config_dict = {} - log_level = config_dict.pop('log_level', "INFO") - tardis.logger.setLevel(getattr(logging, log_level.upper())) + #First let's see if we can find an atom_db anywhere: + if args is not None and args.atom_data is not None: + atom_data_fname = args.atom_data + if 'atom_data' in yaml_dict.keys(): + logger.warn('Ignoring atom_data given in config file (%s)', yaml_dict['atom_data']) + elif 'atom_data' in yaml_dict.keys(): + atom_data_fname = yaml_dict['atom_data'] + else: + raise TardisConfigError('No atom_data key found in config or command line') + + logger.info('Reading Atomic Data from %s', atom_data_fname) + atom_data = atomic.AtomData.from_hdf5(atom_data_fname) + config_dict['atom_data'] = atom_data + #Next finding the time of explosion + if 'time_explosion' not in yaml_dict.keys(): + raise TardisConfigError('No time_explosion found - essential for simulation') + else: + time_explosion = parse2quantity(yaml_dict['time_explosion']).to('s').value - if model_type != 'radial1d': - raise ValueError("Only supporting 'radial1d' at the moment") + config_dict['time_explosion'] = time_explosion - # reading time since explosion - time_explosion_value, time_explosion_unit = config_dict.pop('time_explosion').split() - self.time_explosion = units.Quantity(float(time_explosion_value), time_explosion_unit).to('s').value - # Reading luminosity, special unit log_l_sun is luminosity given in log10 - # of solar units - luminosity_value, luminosity_unit = config_dict.pop('luminosity').split() - if luminosity_unit == 'log_lsun': - self.luminosity_outer = 10 ** ( - float(luminosity_value) + np.log10(constants.L_sun.cgs.value)) + if 'log_lsun' in yaml_dict['luminosity']: + luminosity_value, luminosity_unit = yaml_dict['luminosity'].split() + config_dict['luminosity'] = 10 ** (float(luminosity_value) + np.log10(constants.L_sun.cgs.value)) else: - self.luminosity_outer = units.Quantity( - float(luminosity_value), luminosity_unit).to('erg/s').value + config_dict['luminosity'] = parse2quantity(yaml_dict['luminosity']) - # reading number of shells - no_of_shells = int(config_dict.pop('zones')) - self.no_of_shells = no_of_shells + #Trying to figure out the structure (number of shells) + structure_dict = yaml_dict['model'].pop('structure') - # reading velocities - # set of velocities currently supported are v_inner, v_outer and - # v_sampling linear + #first let's try to see if there's a file keyword + if 'file' in structure_dict.keys(): + density_file_section = structure_dict.pop('file') + v_inner, v_outer, mean_densities, min_shell, max_shell = parse_density_file_section(density_file_section, time_explosion) - v_inner_value, v_inner_unit = config_dict.pop('v_inner').split() - v_inner = units.Quantity( - float(v_inner_value), v_inner_unit).to('cm/s').value + no_of_shells = len(v_inner) + if structure_dict != {}: + logger.warn('Accepted file for structure (density/velocity) structure ignoring all other arguments: \n%s\n', + pprint.pformat(structure_dict, indent=4)) - v_outer_value, v_outer_unit = config_dict.pop('v_outer').split() - v_outer = units.Quantity( - float(v_outer_value), v_outer_unit).to('cm/s').value + else: + #requiring all keys: no_of_shells, velocity, density + if not all([item in structure_dict.keys() for item in ('no_of_shells', 'velocity', 'density')]): + raise TardisConfigError('If file-section is not given to structure-section, one needs to provide all: no_of_shells, velocity, density') - v_sampling = config_dict.pop('v_sampling') + no_of_shells = structure_dict['no_of_shells'] - self.set_velocities( - v_inner=v_inner, v_outer=v_outer, v_sampling=v_sampling) + v_inner, v_outer = parse_velocity_section(structure_dict['velocity'], no_of_shells) + mean_densities = parse_density_section(structure_dict['density'], no_of_shells, v_inner, v_outer, time_explosion) - density_set = config_dict.pop('density') + config_dict['v_inner'] = v_inner + config_dict['v_outer'] = v_outer + config_dict['mean_densities'] = mean_densities + config_dict['no_of_shells'] = no_of_shells - if density_set == 'w7_branch85': - self.densities = calculate_w7_branch85_densities( - self.velocities, - self.time_explosion) - elif density_set == 'exponential': - #TODO:Add here the function call which generates the exponential density profile. The easy way from tonight don't work as expected!! - if not (('exponential_n_factor' in config_dict) and ('exponential_rho0' in config_dict)): - raise ValueError( - 'If density=exponential is set the exponential_n_factor(float) and exponential_rho_0 have to be specified.') + #Now that the structure section is parsed we move on to the abundances - self.exponential_n_factor = float(config_dict.pop('exponential_n_factor')) - self.exponential_rho_0 = float(config_dict.pop('exponential_rho0')) + abundances_dict = yaml_dict['model']['abundances'].copy() + #TODO: columns are now until Z=120 + species_pattern = re.compile('\s*([a-zA-Z]*)(\d*)\s*') + abundances = pd.DataFrame(columns=np.arange(1, 120), index=pd.Index(np.arange(no_of_shells), name='shells')) - self.densities = calculate_exponential_densities(self.velocities, v_inner, - self.exponential_rho_0, self.exponential_n_factor) + if 'file' in abundances_dict.keys(): + abundance_file_dict = abundances_dict.pop('file') + parse_abundance_file_section(abundance_file_dict, abundances, min_shell, max_shell) + if 'abundance_set' in abundances_dict.keys(): + abundance_set_dict = abundances_dict.pop('abundance_set') + print "abundance set not implemented currently" +# abundance_set = abundance_dict.pop('abundance_set', None) +# if abundance_set == 'lucy99': +# abundances = read_lucy99_abundances() +# elif abundance_set is not None: +# raise ValueError("Currently only 'lucy99' abundance_set supported") - else: - try: - density = float(density_set) - except ValueError: - raise ValueError( - 'Currently only density = w7_branch85 or density = exponential,' - ' or specifying a uniform density (with a number) are supported') - self.densities = np.ones(self.no_of_shells) * density - - - # reading plasma type - self.plasma_type = config_dict.pop('plasma_type') - # - self.radiative_rates_type = config_dict.pop('radiative_rates_type') - - # reading initial t_rad - if 'initial_t_rad' in config_dict: - initial_t_rad_value, initial_t_rad_unit = config_dict.pop('initial_t_rad').split() - self.initial_t_rad = units.Quantity(float(initial_t_rad_value), initial_t_rad_unit).to('K').value - logger.info('Selected %g K as initial temperature', self.initial_t_rad) - else: - logger.warn('No initial shell temperature specified (initial_t_rad) - using default 10000 K') - # reading line interaction type - self.line_interaction_type = config_dict.pop('line_interaction_type') - atom_data_file = config_dict.pop('atom_data_file', None) - if atom_data_file is None: - raise ValueError("Please specify a filename with the keyword 'atom_data_file'") + nlte_species = [] + if 'nlte_species' in abundances_dict.keys(): + nlte_species_list = abundances_dict.pop('nlte_species') + for species_symbol in nlte_species_list: + species_match = species_pattern.match(species_symbol) + if species_match is None: + raise ValueError( + "'nlte_species' element %s could not be matched to a valid species notation (e.g. Ca2)") + species_element, species_ion = species_match.groups() + species_element = reformat_element_symbol(species_element) + if species_element not in atom_data.symbol2atomic_number: + raise ValueError("Element provided in NLTE species %s unknown" % species_element) + nlte_species.append((atom_data.symbol2atomic_number[species_element], int(species_ion) - 1)) - self.atom_data = atomic.AtomData.from_hdf5(atom_data_file) - logger.info("Loaded atom data with UUID=%s", self.atom_data.uuid1) - logger.info("Loaded atom data with MD5=%s", self.atom_data.md5) + for element in abundances_dict: + element_symbol = reformat_element_symbol(element) + if element_symbol not in atom_data.symbol2atomic_number: + raise ValueError('Element %s provided in config unknown' % element_symbol) - if self.plasma_type.lower == 'nebular' and not self.atom_data.has_zeta_data: - raise ValueError('The specified AtomData set does not contain zeta data needed for nebular approximation') - # reading number of packets and iterations - if 'iterations' in config_dict and 'no_of_packets' in config_dict: - self.iterations = int(float(config_dict.pop('iterations'))) - self.no_of_packets = int(float(config_dict.pop('no_of_packets'))) - else: - raise ValueError("'no_of_packets' and 'iterations' needs to be set in the configuration") + z = atom_data.symbol2atomic_number[element_symbol] - last_no_of_packets = config_dict.pop('last_no_of_packets', None) - if last_no_of_packets is not None: - self.last_no_of_packets = int(float(last_no_of_packets)) - logger.info('Last iteration will have %g packets', self.last_no_of_packets) - else: - self.last_no_of_packets = None + abundances[z] = float(abundances_dict[element]) - no_of_virtual_packets = config_dict.pop('no_of_virtual_packets', None) - if no_of_virtual_packets is not None: - self.no_of_virtual_packets = int(float(no_of_virtual_packets)) - logger.info('Activating Virtual packets for last iteration (%g)', self.no_of_virtual_packets) - else: - self.no_of_virtual_packets = None + config_dict['abundances'] = abundances + config_dict['nlte_species'] = nlte_species - spectrum_start_value, spectrum_end_unit = config_dict.pop( - 'spectrum_start').split() - spectrum_start = units.Quantity(float(spectrum_start_value), spectrum_end_unit).to('angstrom', - units.spectral()).value - spectrum_end_value, spectrum_end_unit = config_dict.pop('spectrum_end').split() - spectrum_end = units.Quantity(float(spectrum_end_value), spectrum_end_unit).to('angstrom', - units.spectral()).value + ########### DOING PLASMA SECTION ############### - self.spectrum_bins = int(float(config_dict.pop('spectrum_bins'))) + plasma_section = yaml_dict.pop('plasma') - if spectrum_end > spectrum_start: - logger.debug('Converted spectrum start/end to angstrom %.4g %.4g', spectrum_start, spectrum_end) - self.spectrum_start = spectrum_start - self.spectrum_end = spectrum_end - else: - logger.warn('Spectrum Start > Spectrum End in wavelength space - flipped them') + config_dict['initial_t_rad'] = parse2quantity(plasma_section['initial_t_rad']).to('K').value + config_dict['initial_t_inner'] = parse2quantity(plasma_section['initial_t_inner']).to('K').value - logger.debug('Converted spectrum start/end to angstrom %.4g %.4g', spectrum_end, spectrum_start) + if plasma_section['plasma_type'] not in ('nebular', 'lte'): + raise TardisConfigError('plasma_type only allowed to be "nebular" or "lte"') + config_dict['plasma_type'] = plasma_section['plasma_type'] - self.spectrum_start = spectrum_end - self.spectrum_end = spectrum_start + if plasma_section['radiative_rates_type'] not in ('nebular', 'lte', 'detailed'): + raise TardisConfigError('radiative_rates_types must be either "nebular", "lte", or "detailed"') + config_dict['radiative_rates_type'] = plasma_section['radiative_rates_type'] - self.spectrum_start_nu = units.Quantity(self.spectrum_end, 'angstrom').to('Hz', units.spectral()) - self.spectrum_end_nu = units.Quantity(self.spectrum_start, 'angstrom').to('Hz', units.spectral()) + if plasma_section['line_interaction_type'] not in ('scatter', 'downbranch', 'macroatom'): + raise TardisConfigError('radiative_rates_types must be either "scatter", "downbranch", or "macroatom"') + config_dict['line_interaction_type'] = plasma_section['line_interaction_type'] - sn_distance = config_dict.pop('sn_distance', None) - if sn_distance is not None: - if sn_distance.strip().lower() == 'lum_density': - logger.info('Luminosity density requested setting distance to sqrt(1/(4*pi))') - self.sn_distance = np.sqrt(1 / (4 * np.pi)) - else: - sn_distance_value, sn_distance_unit = sn_distance.split() - self.sn_distance = units.Quantity(sn_distance_value, sn_distance_unit).to('cm').value - else: - self.sn_distance = None - disable_electron_scattering = config_dict.pop('disable_electron_scattering', None) + config_dict.update(yaml_dict.pop('montecarlo', {})) + + disable_electron_scattering = plasma_section['disable_electron_scattering'] - if disable_electron_scattering is None or disable_electron_scattering.lower().startswith('f'): + if disable_electron_scattering is False: logger.info("Electron scattering switched on") - self.sigma_thomson = None + config_dict['sigma_thomson'] = None else: logger.warn('Disabling electron scattering - this is not physical') - self.sigma_thomson = 1e-200 + config_dict['sigma_thomson'] = 1e-200 + ##### spectrum section ###### + spectrum_section = yaml_dict.pop('spectrum') + spectrum_start = parse2quantity(spectrum_section['start']).to('angstrom', units.spectral()) + spectrum_end = parse2quantity(spectrum_section['end']).to('angstrom', units.spectral()) + spectrum_bins = int(spectrum_section['bins']) - if config_dict != {}: - logger.warn('Not all config options parsed - ignored %s' % config_dict) - def parse_abundance_section(self, abundance_dict): - abundances = {} - abundance_set = abundance_dict.pop('abundance_set', None) - if abundance_set == 'lucy99': - abundances = read_lucy99_abundances() - elif abundance_set is not None: - raise ValueError("Currently only 'lucy99' abundance_set supported") + if spectrum_end > spectrum_start: + logger.debug('Converted spectrum start/end to angstrom %.4g %.4g', spectrum_start, spectrum_end) + spectrum_start = spectrum_start + spectrum_end = spectrum_end - self.nlte_species = [] - species_pattern = re.compile('\s*([a-zA-Z]*)(\d*)\s*') - if 'nlte_species' in abundance_dict: - for species_symbol in abundance_dict.pop('nlte_species').split(','): - species_match = species_pattern.match(species_symbol) - if species_match is None: - raise ValueError( - "'nlte_species' element %s could not be matched to a valid species notation (e.g. Ca2)") - species_element, species_ion = species_match.groups() - species_element = reformat_element_symbol(species_element) - if species_element not in self.atom_data.symbol2atomic_number: - raise ValueError("Element provided in NLTE species %s unknown" % species_element) - self.nlte_species.append((self.atom_data.symbol2atomic_number[species_element], int(species_ion) - 1)) + else: + logger.warn('Spectrum Start > Spectrum End in wavelength space - flipped them') - for element in abundance_dict: - element_symbol = reformat_element_symbol(element) - if element_symbol not in self.atom_data.symbol2atomic_number: - raise ValueError('Element %s provided in config unknown' % element_symbol) - if element_symbol in abundances: - logger.debug('Element %s already in abundances - overwriting %g with %g', (element_symbol, - abundances[element_symbol], - abundance_dict[element])) - abundances[element_symbol] = float(abundance_dict[element]) + logger.debug('Converted spectrum start/end to angstrom %.4g %.4g', spectrum_end, spectrum_start) + tmp = spectrum_start + spectrum_start = spectrum_end + spectrum_end = tmp + config_dict['spectrum_start'] = spectrum_start + config_dict['spectrum_end'] = spectrum_end + config_dict['spectrum_bins'] = spectrum_bins - self.set_abundances(abundances) + config_dict['spectrum_start_nu'] = spectrum_end.to('Hz', units.spectral()) + config_dict['spectrum_end_nu'] = spectrum_start.to('Hz', units.spectral()) + sn_distance = spectrum_section.pop('sn_distance', None) - @property - def line_interaction_type(self): - return self._line_interaction_type + if sn_distance is not None: + if sn_distance.strip().lower() == 'lum_density': + logger.info('Luminosity density requested - setting distance to sqrt(1/(4*pi))') + config_dict['lum_density'] = True + config_dict['sn_distance'] = np.sqrt(1 / (4 * np.pi)) + else: + config_dict['sn_distance'] = parse2quantity(sn_distance).to('cm').value + config_dict['lum_density'] = False + else: + config_dict['sn_distance'] = None - @line_interaction_type.setter - def line_interaction_type(self, value): - if value not in ('scatter', 'downbranch', 'macroatom'): - raise ValueError('line_interaction_type can only be "scatter", "downbranch", or "macroatom"') - self._line_interaction_type = value + return cls(config_dict) - def set_velocities(self, velocities=None, v_inner=None, v_outer=None, v_sampling='linear'): - """ - Setting the velocities - :param velocities: - :param v_inner: - :param v_outer: - :param v_sampling: - """ - if self.no_of_shells is None: - raise ValueError('Can not set abundances before number of shells have been set') - if v_sampling == 'linear': - self.velocities = np.linspace( - v_inner, v_outer, self.no_of_shells + 1) - else: - raise ValueError('Currently only v_sampling = linear is possible') - def set_abundances(self, abundances): - """ - Setting the abundances - abundances: `dict` or `list` - if a dict is given the assumed mode is uniform, if a list is given it must contain only lists - """ - if self.no_of_shells is None: - raise ValueError('Can not set abundances before number of shells have been set') - if isinstance(abundances, dict): - self.abundances = [abundances] * self.no_of_shells + def __init__(self, config_dict): - @property - def luminosity_outer(self): - return self._luminosity_outer + for key in config_dict: + setattr(self, key, config_dict[key]) - @luminosity_outer.setter - def luminosity_outer(self, value): - self._luminosity_outer = value + self.number_densities = self.calculate_number_densities() - def set_densities(self, densities): - """ - :param densities: - :return: - """ - self.densities = densities - def final_preparation(self): - """ - Does the final preparation for the configuration object + def calculate_number_densities(self): + abundances = self.abundances + for atomic_number in abundances: + + if all(abundances[atomic_number].isnull()): + del abundances[atomic_number] + continue + else: + abundances[abundances[atomic_number].isnull()] == 0.0 + + #normalizing + abundances = abundances.divide(abundances.sum(axis=1), axis=0) + atom_mass = self.atom_data.atom_data.ix[abundances.columns].mass + number_densities = (abundances.mul(self.mean_densities, axis=0)).divide(atom_mass) + + return number_densities - Generates the atoms needed in the simulation - :return: - """ - self.selected_atoms = set() - for shell_abundance in self.abundances: - self.selected_atoms.update(shell_abundance.keys()) - # self.required_atomic_number = - # set([atom_data.symbol2atomic_number[item] for item in - # required_atoms]) diff --git a/tardis/model_radial_oned.py b/tardis/model_radial_oned.py index 24460a94667..c2dc555a0c1 100644 --- a/tardis/model_radial_oned.py +++ b/tardis/model_radial_oned.py @@ -1,12 +1,13 @@ # building of radial_oned_model import numpy as np -import plasma, atomic, packet_source +import plasma, packet_source import logging +import pylab import pandas as pd from astropy import constants, units -from copy import deepcopy +import montecarlo_multizone import os import yaml import pdb @@ -72,8 +73,6 @@ class Radial1DModel(object): def __init__(self, tardis_config): #final preparation for configuration object - tardis_config.final_preparation() - self.tardis_config = tardis_config self.atom_data = tardis_config.atom_data @@ -87,26 +86,23 @@ def __init__(self, tardis_config): self.time_explosion = tardis_config.time_explosion #initializing velocities and radii - self.v_inner = tardis_config.velocities[:-1] - self.v_outer = tardis_config.velocities[1:] + self.v_inner = tardis_config.v_inner + self.v_outer = tardis_config.v_outer self.r_inner = self.v_inner * self.time_explosion self.r_outer = self.v_outer * self.time_explosion self.r_middle = 0.5 * (self.r_inner + self.r_outer) self.volumes = (4. / 3) * np.pi * (self.r_outer ** 3 - self.r_inner ** 3) - #initializing densities - assert len(tardis_config.densities) == self.no_of_shells - - self.densities_middle = tardis_config.densities - self.time_of_simulation = tardis_config.time_of_simulation + self.mean_densities = tardis_config.mean_densities - self.t_inner = tardis_config.t_inner + self.t_inner = tardis_config.initial_t_inner - self.luminosity_outer = tardis_config.luminosity_outer + self.luminosity_outer = tardis_config.luminosity self.no_of_packets = tardis_config.no_of_packets + self.current_no_of_packets = tardis_config.no_of_packets self.iterations = tardis_config.iterations self.sigma_thomson = tardis_config.sigma_thomson @@ -139,26 +135,21 @@ def __init__(self, tardis_config): #initializing abundances self.abundances = tardis_config.abundances - self.number_densities = [] - for abundance, density in zip(self.abundances, self.densities_middle): - self.number_densities.append(calculate_atom_number_densities(self.atom_data, abundance, density)) - self.selected_atomic_numbers = self.number_densities[0].index.values.astype(int) + self.number_densities = tardis_config.number_densities + + self.selected_atomic_numbers = self.number_densities.columns self.line_interaction_type = tardis_config.line_interaction_type #setting dilution factors - if tardis_config.ws is None: - self.ws = 0.5 * (1 - np.sqrt(1 - self.r_inner[0] ** 2 / self.r_middle ** 2)) - else: - self.ws = np.array([(0.5 * (1 - np.sqrt(1 - self.r_inner[0] ** 2 / self.r_middle[i] ** 2))) if w < 0 \ - else w for i, w in enumerate(tardis_config.ws)]) + self.ws = 0.5 * (1 - np.sqrt(1 - self.r_inner[0] ** 2 / self.r_middle ** 2)) #initializing temperatures if np.isscalar(tardis_config.initial_t_rad): - self.t_rads = [tardis_config.initial_t_rad] * self.no_of_shells + self.t_rads = np.ones(self.no_of_shells) * tardis_config.initial_t_rad else: assert len(tardis_config.initial_t_rad) == self.no_of_shells self.t_rads = np.array(tardis_config.initial_t_rad, dtype=np.float64) @@ -189,6 +180,10 @@ def line_interaction_type(self, value): #final preparation for atom_data object - currently building data self.atom_data.prepare_atom_data(self.selected_atomic_numbers, line_interaction_type=self.line_interaction_type, max_ion_number=None) + if self.line_interaction_id in (1,2): + self.calculate_transition_probabilities() + + @property def t_inner(self): @@ -201,13 +196,13 @@ def t_inner(self, value): self.time_of_simulation = 1 / self.luminosity_inner - def create_packets(self, no_of_packets=None): + def create_packets(self): #Energy emitted from the inner boundary self.emitted_inner_energy = 4 * np.pi * constants.sigma_sb.cgs.value * self.r_inner[0] ** 2 * ( self.t_inner) ** 4 - if no_of_packets is None: - no_of_packets = self.no_of_packets + + no_of_packets = self.current_no_of_packets self.packet_src.create_packets(no_of_packets, self.t_inner) def initialize_plasmas(self): @@ -219,8 +214,8 @@ def initialize_plasmas(self): self.transition_probabilities = [] if self.plasma_type == 'lte': - for i, (current_abundances, current_t_rad) in \ - enumerate(zip(self.number_densities, self.t_rads)): + for i, ((tmp_index, current_abundances), current_t_rad) in \ + enumerate(zip(self.number_densities.iterrows(), self.t_rads)): current_plasma = self.plasma_class(current_abundances, self.atom_data, self.time_explosion, nlte_species=self.tardis_config.nlte_species, zone_id=i) logger.debug('Initializing Shell %d Plasma with T=%.3f' % (i, current_t_rad)) @@ -235,16 +230,13 @@ def initialize_plasmas(self): current_plasma.update_radiationfield(current_t_rad) self.tau_sobolevs[i] = current_plasma.tau_sobolevs - #TODO change this - if self.line_interaction_id in (1, 2): - self.transition_probabilities.append(current_plasma.update_macro_atom().values) self.plasmas.append(current_plasma) elif self.plasma_type == 'nebular': - for i, (current_abundances, current_t_rad, current_w) in \ - enumerate(zip(self.number_densities, self.t_rads, self.ws)): + for i, ((tmp_index, current_abundances), current_t_rad, current_w) in \ + enumerate(zip(self.number_densities.iterrows(), self.t_rads, self.ws)): current_plasma = self.plasma_class(current_abundances, self.atom_data, self.time_explosion, nlte_species=self.tardis_config.nlte_species, zone_id=i) logger.debug('Initializing Shell %d Plasma with T=%.3f W=%.4f' % (i, current_t_rad, current_w)) @@ -261,19 +253,24 @@ def initialize_plasmas(self): self.tau_sobolevs[i] = current_plasma.tau_sobolevs - if self.line_interaction_id in (1, 2): - self.transition_probabilities.append(current_plasma.update_macro_atom().values) - self.plasmas.append(current_plasma) self.tau_sobolevs = np.array(self.tau_sobolevs, dtype=float) self.j_blues = np.zeros_like(self.tau_sobolevs) if self.line_interaction_id in (1, 2): - self.transition_probabilities = np.array(self.transition_probabilities, dtype=np.float64) + self.calculate_transition_probabilities() # update plasmas + def calculate_transition_probabilities(self): + self.transition_probabilities = [] + + for current_plasma in self.plasmas: + self.transition_probabilities.append(current_plasma.update_macro_atom().values) + + self.transition_probabilities = np.array(self.transition_probabilities, dtype=np.float64) + def calculate_updated_radiationfield(self, nubar_estimator, j_estimator): """ Calculate an updated radiation field from the :math:`\\bar{nu}_\\textrm{estimator}` and :math:`\\J_\\textrm{estimator}` @@ -308,10 +305,9 @@ def normalize_j_blues(self): self.j_blues *= norm_factor - def update_plasmas(self, updated_t_rads, updated_ws=None): + def update_plasmas(self): if self.plasma_type == 'lte': - self.t_rads = updated_t_rads - for i, (current_plasma, new_trad) in enumerate(zip(self.plasmas, updated_t_rads)): + for i, (current_plasma, new_trad) in enumerate(zip(self.plasmas, self.t_rads)): logger.debug('Updating Shell %d Plasma with T=%.3f' % (i, new_trad)) if self.radiative_rates_type == 'lte': j_blues = plasma.intensity_black_body(self.atom_data.lines.nu.values, new_trad) @@ -327,11 +323,8 @@ def update_plasmas(self, updated_t_rads, updated_ws=None): self.tau_sobolevs[i] = current_plasma.tau_sobolevs elif self.plasma_type == 'nebular': - self.t_rads = updated_t_rads - self.ws = updated_ws - for i, (current_plasma, new_trad, new_ws) in enumerate(zip(self.plasmas, updated_t_rads, updated_ws)): + for i, (current_plasma, new_trad, new_ws) in enumerate(zip(self.plasmas, self.t_rads, self.ws)): logger.debug('Updating Shell %d Plasma with T=%.3f W=%.4f' % (i, new_trad, new_ws)) - if self.radiative_rates_type == 'lte': j_blues = plasma.intensity_black_body(self.atom_data.lines.nu.values, new_trad) elif self.radiative_rates_type == 'nebular': @@ -348,16 +341,19 @@ def update_plasmas(self, updated_t_rads, updated_ws=None): current_plasma.update_radiationfield(new_trad, new_ws) self.tau_sobolevs[i] = current_plasma.tau_sobolevs + if self.line_interaction_id in (1, 2): + self.calculate_transition_probabilities() + - def calculate_spectrum(self, out_nu, out_energy, distance=None): - self.out_nu = out_nu - self.out_energy = out_energy + def calculate_spectrum(self): - if distance is None: + if self.tardis_config.sn_distance is None: logger.info('Distance to supernova not selected assuming 10 pc for calculation of spectra') distance = units.Quantity(10, 'pc').to('cm').value - - self.spec_flux_nu = np.histogram(out_nu[out_nu > 0], weights=out_energy[out_nu > 0], bins=self.spec_nu_bins)[0] + else: + distance = self.tardis_config.sn_distance + self.spec_flux_nu = np.histogram(self.montecarlo_nu[self.montecarlo_nu > 0], + weights=self.montecarlo_energies[self.montecarlo_energies > 0], bins=self.spec_nu_bins)[0] flux_scale = self.time_of_simulation * (self.spec_nu[1] - self.spec_nu[0]) * (4 * np.pi * distance ** 2) @@ -366,7 +362,8 @@ def calculate_spectrum(self, out_nu, out_energy, distance=None): self.spec_virtual_flux_nu /= flux_scale self.spec_reabsorbed_nu = \ - np.histogram(out_nu[out_nu < 0], weights=out_energy[out_nu < 0], bins=self.spec_nu_bins)[0] + np.histogram(self.montecarlo_nu[self.montecarlo_nu < 0], + weights=self.montecarlo_energies[self.montecarlo_nu < 0], bins=self.spec_nu_bins)[0] self.spec_reabsorbed_nu /= flux_scale self.spec_angstrom = units.Unit('Hz').to('angstrom', self.spec_nu, units.spectral()) @@ -375,6 +372,81 @@ def calculate_spectrum(self, out_nu, out_energy, distance=None): self.spec_reabsorbed_angstrom = (self.spec_reabsorbed_nu * self.spec_nu ** 2 / constants.c.cgs / 1e8) self.spec_virtual_flux_angstrom = (self.spec_virtual_flux_nu * self.spec_nu ** 2 / constants.c.cgs / 1e8) + + def simulate(self, update_radiation_field=True, enable_virtual=False): + self.create_packets() + self.spec_virtual_flux_nu[:] = 0.0 + if enable_virtual: + self.montecarlo_nu, self.montecarlo_energies, self.j_estimators, self.nubar_estimators = \ + montecarlo_multizone.montecarlo_radial1d(self, + virtual_packet_flag=self.tardis_config.no_of_virtual_packets) + else: + self.montecarlo_nu, self.montecarlo_energies, self.j_estimators, self.nubar_estimators = \ + montecarlo_multizone.montecarlo_radial1d(self) + + self.normalize_j_blues() + + + self.calculate_spectrum() + + if update_radiation_field: + self.update_radiationfield() + self.update_plasmas() + + + + def update_radiationfield(self, update_mode='dampened', damping_constant=0.5, log_sampling=5): + """ + Updating radiatiantion field + + Parameters + ---------- + + nubar_estimators : ~np.ndarray + j_estimators : ~np.ndarray + update_mode : 'damped' or 'direct' + + damping_constant : + """ + + updated_t_rads, updated_ws = self.calculate_updated_radiationfield(self.nubar_estimators, self.j_estimators) + + + + if update_mode in ('dampened', 'direct'): + if update_mode == 'direct': + damping_constant = 1.0 + + old_t_rads = self.t_rads.copy() + old_ws = self.ws.copy() + self.t_rads += damping_constant * (updated_t_rads - self.t_rads) + self.ws += damping_constant * (updated_ws - self.ws) + + emitted_energy = self.emitted_inner_energy * \ + np.sum(self.montecarlo_energies[self.montecarlo_energies >= 0]) / 1. + absorbed_energy = self.emitted_inner_energy * \ + np.sum(self.montecarlo_energies[self.montecarlo_energies < 0]) / -1. + + + updated_t_inner = self.t_inner * (emitted_energy / self.luminosity_outer) ** -.25 + + self.t_inner += damping_constant * (updated_t_inner - self.t_inner) + + temperature_logging = pd.DataFrame({'t_rads': old_t_rads, 'updated_t_rads': updated_t_rads, 'new_trads': self.t_rads, + 'ws': old_ws, 'updated_ws': updated_ws, 'new_ws': self.ws}) + temperature_logging.index.name = 'Shell' + + temperature_logging = str(temperature_logging[::log_sampling]) + + temperature_logging = ''.join(['\t%s\n' % item for item in temperature_logging.split('\n')]) + + logger.info('Plasma stratification:\n%s\n', temperature_logging) + logger.info("Luminosity emitted = %.5e Luminosity absorbed = %.5e Luminosity requested = %.5e", emitted_energy, + absorbed_energy, + self.luminosity_outer) + logger.info('Calculating new t_inner = %.3f', updated_t_inner) + + def create_synpp_yaml(self, fname, lines_db=None): if not self.atom_data.has_synpp_refs: raise ValueError( @@ -418,40 +490,26 @@ def create_synpp_yaml(self, fname, lines_db=None): yaml.dump(yaml_reference, file(fname, 'w')) - -def calculate_atom_number_densities(atom_data, abundances, density): - """ - Calculates the atom number density, using the following formula, where Z is the atomic number - and X is the abundance fraction - - .. math:: - N_{Z} = \\frac{\\rho_\\textrm{total}\\times \\textrm{X}_\\textrm{Z}}{m_\\textrm{Z}} - - """ - - #Converting abundances - - - - abundance_fractions = pd.Series(abundances.values(), - index=pd.Index([atom_data.symbol2atomic_number[item] for item in abundances.keys()], - dtype=np.int, name='atomic_number'), name='abundance_fraction') - - - - #Normalizing Abundances - - abundance_sum = abundance_fractions.sum() - - if abs(abundance_sum - 1) > 1e-5: - logger.warn('Abundances do not add up to 1 (Sum = %.4f). Renormalizing', (abundance_sum)) - - abundance_fractions /= abundance_sum - - number_densities = (abundance_fractions * density) / \ - atom_data.atom_data.ix[abundance_fractions.index]['mass'] - - return pd.DataFrame({'abundance_fraction': abundance_fractions, 'number_density': number_densities}) + def plot_spectrum(self, ax=None, mode='wavelength', virtual=True): + if ax is None: + ax = pylab.gca() + if mode == 'wavelength': + x = self.spec_angstrom + if virtual: + y = self.spec_virtual_flux_angstrom + else: + y = self.spec_flux_angstrom + xlabel = 'Wavelength [\AA]' + if self.tardis_config.lum_density: + ylabel = 'Flux [erg s^-1 cm^-2 \AA^-1]' + else: + ylabel = 'Flux [erg s^-1 \AA^-1]' + + + + ax.plot(x, y) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) class ModelHistory(object): diff --git a/tardis/montecarlo_multizone.pyx b/tardis/montecarlo_multizone.pyx index b2235e0d4ec..9c06e4c6651 100644 --- a/tardis/montecarlo_multizone.pyx +++ b/tardis/montecarlo_multizone.pyx @@ -43,12 +43,18 @@ cdef extern from "randomkit/randomkit.h": void rk_seed(unsigned long seed, rk_state *state) float_type_t rk_double(rk_state *state) +cdef rk_state mt_state + cdef np.ndarray x cdef class StorageModel: """ Class for storing the arrays in a cythonized way (as pointers). This ensures fast access during the calculations. """ + + + + cdef np.ndarray packet_nus_a cdef float_type_t*packet_nus cdef np.ndarray packet_mus_a @@ -110,6 +116,9 @@ cdef class StorageModel: def __init__(self, model): + + rk_seed(250819801106, &mt_state) + cdef np.ndarray[float_type_t, ndim=1] packet_nus = model.packet_src.packet_nus self.packet_nus_a = packet_nus self.packet_nus = self.packet_nus_a.data @@ -237,8 +246,7 @@ logger = logging.getLogger(__name__) -cdef rk_state mt_state -rk_seed(250819801106, &mt_state) + diff --git a/tardis/packet_source.py b/tardis/packet_source.py index 92f785b82a1..4aea0c0db5c 100644 --- a/tardis/packet_source.py +++ b/tardis/packet_source.py @@ -45,7 +45,7 @@ def __init__(self, nu_start, nu_end, seed=250819801106, blackbody_sampling=int(1 np.random.seed(seed) - def create_packets(self, number_of_packets, t_rad): + def create_packets(self, number_of_packets, t_rad, seed=None): """ Creating a new random number of packets, with a certain temperature @@ -59,6 +59,9 @@ def create_packets(self, number_of_packets, t_rad): radiation temperature """ + if seed is not None: + np.random.seed(seed) + number_of_packets = int(number_of_packets) self.packet_nus = self.random_blackbody_nu(t_rad, number_of_packets) self.packet_mus = np.sqrt(np.random.random(size=number_of_packets)) diff --git a/tardis/plasma.py b/tardis/plasma.py index 2b75b503a1b..0188b8608f2 100644 --- a/tardis/plasma.py +++ b/tardis/plasma.py @@ -160,7 +160,7 @@ def update_radiationfield(self, t_rad, t_electron=None, n_e_convergence_threshol phis = self.calculate_saha() #initialize electron density with the sum of number densities - electron_density = self.abundances['number_density'].sum() + electron_density = self.abundances.sum() n_e_iterations = 0 @@ -292,7 +292,7 @@ def calculate_ionization_balance(self, phis, electron_density): current_phis = groups.values / electron_density phis_product = np.cumproduct(current_phis) - neutral_atom_density = self.abundances.ix[atomic_number]['number_density'] / (1 + np.sum(phis_product)) + neutral_atom_density = self.abundances.ix[atomic_number] / (1 + np.sum(phis_product)) ion_densities = [neutral_atom_density] + list(neutral_atom_density * phis_product) self.ion_number_density.ix[atomic_number] = ion_densities @@ -523,7 +523,7 @@ def update_radiationfield(self, t_rad, w, t_electron=None, n_e_convergence_thres phis = self.calculate_saha() #initialize electron density with the sum of number densities - electron_density = self.abundances['number_density'].sum() + electron_density = self.abundances.sum() n_e_iterations = 0 @@ -635,7 +635,12 @@ def calculate_saha(self): zeta = Series(index=phis.index) for idx in zeta.index: - zeta.ix[idx] = self.atom_data.zeta_data[idx](self.t_rad) + try: + current_zeta = self.atom_data.zeta_data[idx](self.t_rad) + except KeyError: + current_zeta = 1.0 + + zeta.ix[idx] = current_zeta phis *= self.w * (delta.ix[phis.index] * zeta + self.w * (1 - zeta)) * \ (self.t_electron / self.t_rad) ** .5 diff --git a/tardis/simulation.py b/tardis/simulation.py index 768937983ac..2749e3fb931 100644 --- a/tardis/simulation.py +++ b/tardis/simulation.py @@ -13,62 +13,21 @@ logger = logging.getLogger(__name__) -def run_single_radial1d(radial1d_model): - """ - Run a single 1D radial simulation - - Parameters - ---------- - - radial1d_model : `~tardis.model_radial_oned.Radial1DModel` - - - """ - - out_nu, out_energy, j_estimators, nubar_estimators = montecarlo_multizone.montecarlo_radial1d(radial1d_model, 0) - def run_radial1d(radial1d_model, save_history=None): for i in xrange(radial1d_model.iterations - 1): logger.info('At run %d of %d', i + 1, radial1d_model.iterations) - radial1d_model.create_packets() - out_nu, out_energy, j_estimators, nubar_estimators = montecarlo_multizone.montecarlo_radial1d(radial1d_model) - radial1d_model.normalize_j_blues() + radial1d_model.simulate() #spec_nu_flux = np.histogram(out_nu, weights=out_energy, bins=radial1d_model.spec_virt_nu) if save_history is not None: save_history.store_all(radial1d_model, i) - updated_t_rads, updated_ws = radial1d_model.calculate_updated_radiationfield(nubar_estimators, j_estimators) - - new_trads = 0.5 * (radial1d_model.t_rads + updated_t_rads) - new_ws = 0.5 * (radial1d_model.ws + updated_ws) - - print pd.DataFrame({'t_rads': radial1d_model.t_rads, 'updated_t_rads': updated_t_rads, 'new_trads': new_trads, - 'ws': radial1d_model.ws, 'updated_ws': updated_ws, 'new_ws': new_ws})[::5] - - emitted_energy = radial1d_model.emitted_inner_energy * np.sum(out_energy[out_energy >= 0]) / 1. - absorbed_energy = radial1d_model.emitted_inner_energy * np.sum(out_energy[out_energy < 0]) / 1. - - print "Luminosity emitted = %.5e Luminosity absorbed = %.5e Luminosity requested = %.5e" % (emitted_energy, - absorbed_energy, - radial1d_model.luminosity_outer) - new_t_inner = radial1d_model.t_inner * (emitted_energy / radial1d_model.luminosity_outer) ** -.25 - print "new t_inner = %.2f" % (new_t_inner,) - - radial1d_model.t_inner = 0.5 * (new_t_inner + radial1d_model.t_inner) - - radial1d_model.update_plasmas(new_trads, new_ws) - #spec_nu_flux = np.histogram(out_nu, weights=out_energy, bins=radial1d_model.spec_virt_nu) - #Finished second to last loop running one more time logger.info('Doing last run') if radial1d_model.tardis_config.last_no_of_packets is not None: - radial1d_model.create_packets(radial1d_model.tardis_config.last_no_of_packets) - else: - radial1d_model.create_packets() - out_nu, out_energy, j_estimators, nubar_estimators = montecarlo_multizone.montecarlo_radial1d(radial1d_model, - virtual_packet_flag=radial1d_model.tardis_config.no_of_virtual_packets) - radial1d_model.normalize_j_blues() - radial1d_model.calculate_spectrum(out_nu, out_energy, distance=radial1d_model.tardis_config.sn_distance) + radial1d_model.current_no_of_packets = radial1d_model.tardis_config.last_no_of_packets + + radial1d_model.simulate(enable_virtual=True) + if save_history is not None: save_history.store_all(radial1d_model, radial1d_model.iterations - 1) From bf7e1db74b62d7f8be35c135e346e61f1da21f8d Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Tue, 5 Mar 2013 21:02:46 -0500 Subject: [PATCH 3/5] added the example configuration in the new YAML format --- tardis/data/tardis_example_config.yml | 99 +++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 tardis/data/tardis_example_config.yml diff --git a/tardis/data/tardis_example_config.yml b/tardis/data/tardis_example_config.yml new file mode 100644 index 00000000000..96ec0a03618 --- /dev/null +++ b/tardis/data/tardis_example_config.yml @@ -0,0 +1,99 @@ +#New configuration for TARDIS based on YAML +#IMPORTANT any pure floats need to have a +/- after the e e.g. 2e+5 +#Hopefully pyyaml will fix this soon. +--- +#Currently only simple1d is allowed +config_type: simple1d + +#luminosity any astropy.unit convertible to erg/s +#special unit log_lsun(log(luminosity) - log(L_sun) +luminosity: 9.44 log_lsun + +#time since explosion +time_explosion: 18 day + + +atom_data: kurucz_atom_tmp.h5 + +plasma: + initial_t_inner: 10000 K + initial_t_rad: 10000 K + disable_electron_scattering: no + plasma_type: nebular + radiative_rates_type: nebular + #line interaction type - currently supported are scatter, downbranch and macroatom + line_interaction_type : macroatom + + + +model: + structure: + no_of_shells : 20 + file: + type : artis + name : artis_model.dat + v_lowest: 10000.0 km/s + v_highest: 20000.0 km/s + + + velocity: + type : linear + v_inner : 1e4 km/s + v_outer : 2e4 km/s + + + density: + #showing different configuration options separated by comments + #simple uniform: + #--------------- + type : uniform + value : 1e-12 g/cm^3 + #--------------- + + #branch85_w7 - fit of seven order polynomial to W7 (like Branch 85): + #--------------- + #type : branch85_w7 + #value : 1e-12 + #time_0 : 19.9999584 s# default, no need to change! + #density_coefficient : 3e29 # default, no need to change! + #--------------- + + + abundances: + file: + type : artis + name : artis_abundances.dat + +# abundance_set: +# type: uniform +# name: lucy99 +# nlte_species : [] +# sI : 1.5 +# ca: 1.5 + + + +montecarlo: + seed: 23111963171620 + no_of_packets : 2.e+4 + iterations: 1 + last_no_of_packets: 2.0e+4 + no_of_virtual_packets: 10 + + +spectrum: + start : 500 angstrom + end : 20000 angstrom + bins : 1000 + sn_distance : lum_density + + + + + + + + + + + From e2aec9382ca5fc4e2a2ff50add4e248d15edea76 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 6 Mar 2013 11:02:32 -0500 Subject: [PATCH 4/5] minor change to config_reader --- tardis/config_reader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index 4c74432d281..e3f33828e88 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -12,7 +12,7 @@ import pandas as pd from tardis import atomic import yaml -import tardis + import pdb import pprint From 0216c966f0cfa98449cce423e2b3ad21f60d2d12 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 6 Mar 2013 13:55:41 -0500 Subject: [PATCH 5/5] fixed the calculate density after time routine in artis parser. --- tardis/config_reader.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tardis/config_reader.py b/tardis/config_reader.py index e3f33828e88..d59b7bdc22b 100644 --- a/tardis/config_reader.py +++ b/tardis/config_reader.py @@ -58,8 +58,7 @@ def parse_artis_density(density_file_dict, time_explosion): velocities = units.Quantity(np.append([0], velocities), 'km/s').to('cm/s') mean_densities_0 = units.Quantity(10**mean_densities_0, 'g/cm^3') - mean_densities = calculate_density_after_time(time_of_model.value, mean_densities_0, - time_explosion) + mean_densities = calculate_density_after_time(mean_densities_0, time_of_model.value, time_explosion) #Verifying information