Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Restructure the MontecarloRunner #636

Closed
wants to merge 15 commits into from
Closed

Restructure the MontecarloRunner #636

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
b5854bc
Change MontecarloRunner so it doesn't depend on a Configuration object.
ftsamis Jul 26, 2016
5f83ef5
Add a plasma parameter to MontecarloRunner.run
ftsamis Jul 26, 2016
9f4385e
Move and adapt calculate_j_blues from the old Radial1DModel to Montec…
ftsamis Aug 3, 2016
fc8da95
Add the time_of_simulation property to the MontecarloRunner.
ftsamis Aug 4, 2016
f777f8f
Add the plasma as an argument to montecarlo.montecarlo_radial1d.
ftsamis Aug 5, 2016
097454c
Remove every Configuration object usage from the montecarlo.pyx
ftsamis Aug 5, 2016
894510c
Add a from_config classmethod to the MontecarloRunner.
ftsamis Aug 5, 2016
f69c995
montecarlo.pyx: Access atomic_data from the plasma, instead of model.
ftsamis Aug 8, 2016
0e7cf7d
Revert "Move and adapt calculate_j_blues from the old Radial1DModel t…
ftsamis Aug 9, 2016
f410a60
Add a calculate_time_of_simulation and a calculate_luminosity_inner t…
ftsamis Aug 9, 2016
3895e14
Move the sigma_thomson calculation from config_reader to montecarlo.
ftsamis Aug 13, 2016
52071c8
Move the spectrum frequency calculation from config_reader to monteca…
ftsamis Aug 13, 2016
790737c
Remove supernova distance checking from config_dict
ftsamis Aug 13, 2016
931751e
Don't rename virtual_spectrum_range's stop and num to end and samples.
ftsamis Aug 13, 2016
dde78a8
Fix a bug when initializing the logger for montecarlo/base.py
ftsamis Aug 16, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 0 additions & 62 deletions tardis/io/config_reader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand All @@ -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
Expand Down Expand Up @@ -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 = []
Expand Down Expand Up @@ -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)


Expand Down
8 changes: 0 additions & 8 deletions tardis/io/tests/test_config_reader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
68 changes: 55 additions & 13 deletions tardis/montecarlo/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -110,29 +120,30 @@ 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

Parameters
----------

: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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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))
65 changes: 31 additions & 34 deletions tardis/montecarlo/montecarlo.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -105,8 +105,7 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage):
storage.packet_energies = <double*> 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 = <double*> PyArray_DATA(runner.r_inner_cgs)
Expand All @@ -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 = <double*> 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 = <double*> PyArray_DATA(
runner.inverse_electron_densities)
# Switch for continuum processes
Expand All @@ -146,10 +144,10 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage):
storage.l_pop_r = <double*> l_pop_r.data

# Line lists
storage.no_of_lines = model.atom_data.lines.nu.values.size
storage.line_list_nu = <double*> PyArray_DATA(model.atom_data.lines.nu.values)
storage.no_of_lines = plasma.atomic_data.lines.nu.values.size
storage.line_list_nu = <double*> 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 = <double*> PyArray_DATA(
runner.line_lists_tau_sobolevs
Expand All @@ -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 = <double*> 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 = <int_type_t*> PyArray_DATA(
model.atom_data.lines_upper2macro_reference_idx)
plasma.atomic_data.lines_upper2macro_reference_idx)
storage.macro_block_references = <int_type_t*> PyArray_DATA(
model.atom_data.macro_atom_references['block_references'].values)
plasma.atomic_data.macro_atom_references['block_references'].values)
storage.transition_type = <int_type_t*> 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 = <int_type_t*> 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 = <int_type_t*> PyArray_DATA(
model.atom_data.macro_atom_data['lines_idx'].values)
plasma.atomic_data.macro_atom_data['lines_idx'].values)

storage.output_nus = <double*> PyArray_DATA(runner._output_nu)
storage.output_energies = <double*> PyArray_DATA(runner._output_energy)
Expand All @@ -203,25 +201,25 @@ cdef initialize_storage_model(model, runner, storage_model_t *storage):
storage.js = <double*> PyArray_DATA(runner.j_estimator)
storage.nubars = <double*> 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 = <double*> 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 = <double*> 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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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