From 4d54577a4ba0fc5523977302fe0f79dde08a038f Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Tue, 7 Jul 2015 21:02:55 +0100 Subject: [PATCH 01/13] New NLTE system. --- tardis/model.py | 6 +- tardis/plasma/base.py | 3 +- tardis/plasma/properties/atomic.py | 11 +- tardis/plasma/properties/level_population.py | 71 +------- .../plasma/properties/partition_function.py | 152 +++++++++++++++++- tardis/plasma/properties/plasma_input.py | 9 +- .../plasma/properties/property_collections.py | 13 +- .../plasma/properties/radiative_properties.py | 19 ++- tardis/plasma/standard_plasmas.py | 26 ++- 9 files changed, 222 insertions(+), 88 deletions(-) diff --git a/tardis/model.py b/tardis/model.py index 1ad6f75db11..c9567ebf1e6 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -135,6 +135,8 @@ def __init__(self, tardis_config): line_interaction_type=tardis_config.plasma.line_interaction_type, link_t_rad_t_electron=0.9) + self.previous_iteration_beta_sobolevs = np.ones((len(self.plasma_array.lines), len(self.plasma_array.t_electron))) + self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_virtual = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_reabsorbed = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) @@ -236,8 +238,10 @@ def calculate_updated_radiationfield(self, nubar_estimator, j_estimator): def update_plasmas(self, initialize_nlte=False): self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, j_blues=self.j_blues, - initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05) + initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05, + previous_iteration_beta_sobolevs = self.previous_iteration_beta_sobolevs) + self.previous_iteration_beta_sobolevs = self.plasma_array.beta_sobolevs if self.tardis_config.plasma.line_interaction_type in ('downbranch', 'macroatom'): self.transition_probabilities = self.plasma_array.transition_probabilities diff --git a/tardis/plasma/base.py b/tardis/plasma/base.py index 5b58846c46b..e52abeaf1ec 100644 --- a/tardis/plasma/base.py +++ b/tardis/plasma/base.py @@ -205,6 +205,7 @@ class StandardPlasma(BasePlasma): def __init__(self, number_densities, atom_data, time_explosion, delta_treatment=None, nlte_config=None, ionization_mode='lte', excitation_mode='lte', w=None, - link_t_rad_t_electron=0.9): + link_t_rad_t_electron=0.9, nlte_species=None, + previous_beta_sobolevs=None): pass diff --git a/tardis/plasma/properties/atomic.py b/tardis/plasma/properties/atomic.py index 9dd8498c35c..d021bf979b3 100644 --- a/tardis/plasma/properties/atomic.py +++ b/tardis/plasma/properties/atomic.py @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) __all__ = ['Levels', 'Lines', 'LinesLowerLevelIndex', 'LinesUpperLevelIndex', - 'AtomicMass', 'IonizationData', 'ZetaData'] + 'AtomicMass', 'IonizationData', 'ZetaData', 'NLTEData'] class BaseAtomicDataProperty(ProcessingPlasmaProperty): __metaclass__ = ABCMeta @@ -173,3 +173,12 @@ def _filter_atomic_property(self, zeta_data, selected_atoms): def _set_index(self, zeta_data, atomic_data): return zeta_data.set_index(['atomic_number', 'ion_number']) +class NLTEData(ProcessingPlasmaProperty): + outputs = ('nlte_data',) + + def calculate(self, atomic_data): + if getattr(self, self.outputs[0]) is not None: + return (getattr(self, self.outputs[0]),) + else: + return atomic_data.nlte_data + diff --git a/tardis/plasma/properties/level_population.py b/tardis/plasma/properties/level_population.py index e31b056b2e6..b1831e6f253 100644 --- a/tardis/plasma/properties/level_population.py +++ b/tardis/plasma/properties/level_population.py @@ -30,73 +30,4 @@ class LevelNumberDensity(ProcessingPlasmaProperty): def calculate(level_population_fraction, ion_number_density): ion_number_density_broadcast = ion_number_density.ix[ level_population_fraction.index.droplevel(2)].values - return level_population_fraction * ion_number_density_broadcast - -class LevelPopulationNLTE(ProcessingPlasmaProperty): - @staticmethod - def calculate(self): - """ - Calculating the NLTE level populations for specific ions - - """ - - if not hasattr(self, 'beta_sobolevs'): - self.beta_sobolevs = np.zeros_like(self.tau_sobolevs.values) - - macro_atom.calculate_beta_sobolev(self.tau_sobolevs.values.ravel(order='F'), - self.beta_sobolevs.ravel(order='F')) - self.beta_sobolevs_precalculated = True - - if self.nlte_config.get('coronal_approximation', False): - beta_sobolevs = np.ones_like(self.beta_sobolevs) - j_blues = np.zeros_like(self.j_blues) - logger.info('using coronal approximation = setting beta_sobolevs to 1 AND j_blues to 0') - else: - beta_sobolevs = self.beta_sobolevs - j_blues = self.j_blues.values - - if self.nlte_config.get('classical_nebular', False): - logger.info('using Classical Nebular = setting beta_sobolevs to 1') - beta_sobolevs = np.ones_like(self.beta_sobolevs) - - for species in self.nlte_config.species: - logger.info('Calculating rates for species %s', species) - number_of_levels = self.atom_data.levels.energy.ix[species].count() - - level_population_proportionalitys = self.level_population_proportionalitys.ix[species].values - lnl = self.atom_data.nlte_data.lines_level_number_lower[species] - lnu = self.atom_data.nlte_data.lines_level_number_upper[species] - - lines_index = self.atom_data.nlte_data.lines_idx[species] - A_uls = self.atom_data.nlte_data.A_uls[species] - B_uls = self.atom_data.nlte_data.B_uls[species] - B_lus = self.atom_data.nlte_data.B_lus[species] - - r_lu_index = lnu * number_of_levels + lnl - r_ul_index = lnl * number_of_levels + lnu - - r_ul_matrix = np.zeros((number_of_levels, number_of_levels, len(self.t_rads)), dtype=np.float64) - r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, len(self.t_rads))) - r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + B_uls[np.newaxis].T * j_blues[lines_index] - r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] - - r_lu_matrix = np.zeros_like(r_ul_matrix) - r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, len(self.t_rads))) - r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * j_blues[lines_index] * beta_sobolevs[lines_index] - - collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ - self.electron_densities.values - - - rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix - - for i in xrange(number_of_levels): - rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) - - rates_matrix[0, :, :] = 1.0 - - x = np.zeros(rates_matrix.shape[0]) - x[0] = 1.0 - for i in xrange(len(self.t_rads)): - relative_level_population_proportionalitys = np.linalg.solve(rates_matrix[:, :, i], x) - self.level_population_proportionalitys[i].ix[species] = relative_level_population_proportionalitys * self.ion_populations[i].ix[species] \ No newline at end of file + return level_population_fraction * ion_number_density_broadcast \ No newline at end of file diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index a14789893a8..b058fe1f74a 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -8,7 +8,8 @@ logger = logging.getLogger(__name__) __all__ = ['LevelBoltzmannFactorLTE', 'LevelBoltzmannFactorDiluteLTE', - 'PartitionFunction'] + 'LevelBoltzmannFactorNoNLTE', 'LevelBoltzmannFactorNLTECoronal', + 'LevelBoltzmannFactorNLTEGeneral', 'PartitionFunction'] class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): """ @@ -19,7 +20,7 @@ class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): """ - outputs = ('level_boltzmann_factor',) + outputs = ('general_level_boltzmann_factor',) latex_formula = r'$g_{i, j, k} e^{E_{i, j, k} \times \beta_\textrm{rad}}$' @staticmethod @@ -36,7 +37,7 @@ def calculate(levels, beta_rad): class LevelBoltzmannFactorDiluteLTE(ProcessingPlasmaProperty): - outputs = ('level_boltzmann_factor',) + outputs = ('general_level_boltzmann_factor',) @staticmethod def calculate(levels, beta_rad, w): @@ -51,6 +52,151 @@ def calculate(levels, beta_rad, w): level_boltzmann_factor[~levels.metastable] *= w return level_boltzmann_factor +class LevelBoltzmannFactorNoNLTE(ProcessingPlasmaProperty): + + outputs = ('level_boltzmann_factor',) + + @staticmethod + def calculate(general_level_boltzmann_factor): + return general_level_boltzmann_factor + +class LevelBoltzmannFactorNLTECoronal(ProcessingPlasmaProperty): + + outputs = ('level_boltzmann_factor',) + + @staticmethod + def calculate(t_electron, lines, atomic_data, nlte_data, + general_level_boltzmann_factor, nlte_species): + """ + Calculating the NLTE level populations for specific ions + + """ + + beta_sobolevs = np.ones((len(lines), len(t_electron))) + j_blues = np.zeros((len(lines), len(t_electron))) + + for species in nlte_species: + logger.info('Calculating rates for species %s', species) + number_of_levels = atomic_data.levels.energy.ix[species].count() + + lnl = nlte_data.lines_level_number_lower[species] + lnu = nlte_data.lines_level_number_upper[species] + + lines_index = nlte_data.lines_idx[species] + A_uls = nlte_data.A_uls[species] + B_uls = nlte_data.B_uls[species] + B_lus = nlte_data.B_lus[species] + + r_lu_index = lnu * number_of_levels + lnl + r_ul_index = lnl * number_of_levels + lnu + + r_ul_matrix = np.zeros((number_of_levels, number_of_levels, + len(t_electron)), dtype=np.float64) + r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ + B_uls[np.newaxis].T * j_blues[lines_index] + r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] + + r_lu_matrix = np.zeros_like(r_ul_matrix) + r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ + j_blues[lines_index] * beta_sobolevs[lines_index] + +# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ +# self.electron_densities.values + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + + rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix + + for i in xrange(number_of_levels): + rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) + + rates_matrix[0, :, :] = 1.0 + + x = np.zeros(rates_matrix.shape[0]) + x[0] = 1.0 + for i in xrange(len(t_electron)): + level_boltzmann_factor = \ + np.linalg.solve(rates_matrix[:, :, i], x) + general_level_boltzmann_factor[i].ix[species] = \ + level_boltzmann_factor + return general_level_boltzmann_factor + +class LevelBoltzmannFactorNLTEClassicalNebular(ProcessingPlasmaProperty): + pass + +class LevelBoltzmannFactorNLTEGeneral(ProcessingPlasmaProperty): + outputs = ('level_boltzmann_factor',) + + @staticmethod + def calculate(t_electron, lines, atomic_data, nlte_data, + general_level_boltzmann_factor, nlte_species, j_blues, + previous_beta_sobolevs, lte_j_blues): + """ + Calculating the NLTE level populations for specific ions + + """ + if previous_beta_sobolevs is None: + beta_sobolevs = np.ones((len(lines), len(t_electron))) + else: + beta_sobolevs = previous_beta_sobolevs + + if len(j_blues)==0: + j_blues = lte_j_blues + + for species in nlte_species: + logger.info('Calculating rates for species %s', species) + number_of_levels = atomic_data.levels.energy.ix[species].count() + + lnl = nlte_data.lines_level_number_lower[species] + lnu = nlte_data.lines_level_number_upper[species] + + lines_index = nlte_data.lines_idx[species] + A_uls = nlte_data.A_uls[species] + B_uls = nlte_data.B_uls[species] + B_lus = nlte_data.B_lus[species] + + r_lu_index = lnu * number_of_levels + lnl + r_ul_index = lnl * number_of_levels + lnu + + r_ul_matrix = np.zeros((number_of_levels, number_of_levels, + len(t_electron)), dtype=np.float64) + r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ + B_uls[np.newaxis].T * j_blues[lines_index] + r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] + + r_lu_matrix = np.zeros_like(r_ul_matrix) + r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ + j_blues[lines_index] * beta_sobolevs[lines_index] + +# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ +# self.electron_densities.values + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + + rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix + + for i in xrange(number_of_levels): + rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) + + rates_matrix[0, :, :] = 1.0 + + x = np.zeros(rates_matrix.shape[0]) + x[0] = 1.0 + for i in xrange(len(t_electron)): + level_boltzmann_factor = \ + np.linalg.solve(rates_matrix[:, :, i], x) + general_level_boltzmann_factor[i].ix[species] = \ + level_boltzmann_factor + return general_level_boltzmann_factor + class PartitionFunction(ProcessingPlasmaProperty): outputs = ('partition_function',) latex_outputs = '$Z_{i, j}$' diff --git a/tardis/plasma/properties/plasma_input.py b/tardis/plasma/properties/plasma_input.py index b9fb4792fa8..17fd765e2bd 100644 --- a/tardis/plasma/properties/plasma_input.py +++ b/tardis/plasma/properties/plasma_input.py @@ -4,7 +4,8 @@ __all__ = ['TRadiative', 'DilutionFactor', 'AtomicData', 'Abundance', 'Density', 'TimeExplosion', 'JBlues', 'LinkTRadTElectron', - 'RadiationFieldCorrectionInput'] + 'RadiationFieldCorrectionInput', 'NLTESpecies', + 'PreviousIterationBetaSobolevs'] class Input(BasePlasmaProperty): @@ -65,3 +66,9 @@ class JBlues(DataFrameInput): class LinkTRadTElectron(StaticInput): outputs = ('link_t_rad_t_electron',) + +class NLTESpecies(StaticInput): + outputs = ('nlte_species',) + +class PreviousIterationBetaSobolevs(DynamicInput): + outputs = ('previous_beta_sobolevs',) diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index 6f5294d5cad..5d63a22b92b 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -7,14 +7,17 @@ TransitionProbabilities, StimulatedEmissionFactor, SelectedAtoms, PhiGeneral, PhiSahaNebular, LevelBoltzmannFactorDiluteLTE, DilutionFactor, ZetaData, ElectronTemperature, LinkTRadTElectron, BetaElectron, - RadiationFieldCorrection, RadiationFieldCorrectionInput) + RadiationFieldCorrection, RadiationFieldCorrectionInput, + LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTECoronal, NLTEData, + NLTESpecies, PreviousIterationBetaSobolevs, LTEJBlues, + LevelBoltzmannFactorNLTEGeneral) class PlasmaPropertyCollection(list): pass basic_inputs = PlasmaPropertyCollection([TRadiative, Abundance, Density, TimeExplosion, AtomicData, JBlues, DilutionFactor, LinkTRadTElectron, - RadiationFieldCorrectionInput]) + RadiationFieldCorrectionInput, NLTESpecies, PreviousIterationBetaSobolevs]) basic_properties = PlasmaPropertyCollection([BetaRadiation, Levels, Lines, AtomicMass, LevelPopulation, PartitionFunction, GElectron, IonizationData, NumberDensity, LinesLowerLevelIndex, @@ -28,3 +31,9 @@ class PlasmaPropertyCollection(list): ZetaData, ElectronTemperature, BetaElectron, RadiationFieldCorrection]) dilute_lte_excitation_properties = PlasmaPropertyCollection([ LevelBoltzmannFactorDiluteLTE]) +non_nlte_properties = PlasmaPropertyCollection([LevelBoltzmannFactorNoNLTE]) +nlte_coronal_properties = PlasmaPropertyCollection([ + LevelBoltzmannFactorNLTECoronal, NLTEData, NLTESpecies]) +nlte_general_properties = PlasmaPropertyCollection([ + LevelBoltzmannFactorNLTEGeneral, NLTEData, NLTESpecies, + PreviousIterationBetaSobolevs, LTEJBlues]) diff --git a/tardis/plasma/properties/radiative_properties.py b/tardis/plasma/properties/radiative_properties.py index d27dd406517..62b09cc1f0f 100644 --- a/tardis/plasma/properties/radiative_properties.py +++ b/tardis/plasma/properties/radiative_properties.py @@ -10,7 +10,7 @@ logger = logging.getLogger(__name__) __all__ = ['StimulatedEmissionFactor', 'TauSobolev', 'BetaSobolev', - 'TransitionProbabilities'] + 'TransitionProbabilities', 'LTEJBlues'] class StimulatedEmissionFactor(ProcessingPlasmaProperty): @@ -159,3 +159,20 @@ def calculate(self, atomic_data, beta_sobolev, j_blues, index=macro_atom_data.transition_line_id, columns=tau_sobolevs.columns) return transition_probabilities + +class LTEJBlues(ProcessingPlasmaProperty): + outputs = ('lte_j_blues',) + + @staticmethod + def calculate(lines, beta_rad): + beta_rad = pd.Series(beta_rad) + nu = pd.Series(lines.nu) + h = const.h.cgs.value + c = const.c.cgs.value + df = pd.DataFrame(1, index=nu.index, columns=beta_rad.index) + df = df.multiply(nu, axis='index') * beta_rad + exponential = (np.exp(h * df) - 1)**(-1) + remainder = (2 * (h * lines.nu.values ** 3) / + (c ** 2)) + j_blues = exponential.multiply(remainder, axis=0) + return pd.DataFrame(j_blues, index=lines.index, columns=beta_rad.index) \ No newline at end of file diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 603ce0a5108..5ea8c7fe5eb 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -6,7 +6,8 @@ from tardis.plasma.properties.property_collections import (basic_inputs, basic_properties, lte_excitation_properties, lte_ionization_properties, macro_atom_properties, dilute_lte_excitation_properties, - nebular_ionization_properties) + nebular_ionization_properties, non_nlte_properties, + nlte_coronal_properties, nlte_general_properties) logger = logging.getLogger(__name__) @@ -15,13 +16,15 @@ class LTEPlasma(BasePlasma): def __init__(self, t_rad, abundance, density, time_explosion, atomic_data, j_blues, link_t_rad_t_electron=0.9, delta_treatment=None): plasma_modules = basic_inputs + basic_properties + \ - lte_excitation_properties + lte_ionization_properties + lte_excitation_properties + lte_ionization_properties + \ + non_nlte_properties super(LTEPlasma, self).__init__(plasma_properties=plasma_modules, t_rad=t_rad, abundance=abundance, atomic_data=atomic_data, density=density, time_explosion=time_explosion, j_blues=j_blues, w=None, link_t_rad_t_electron=link_t_rad_t_electron, - delta_input=delta_treatment) + delta_input=delta_treatment, nlte_species=None, + previous_beta_sobolevs=None) class LegacyPlasmaArray(BasePlasma): @@ -43,8 +46,9 @@ def initial_w(self, number_densities): def update_radiationfield(self, t_rad, ws, j_blues, t_electrons=None, n_e_convergence_threshold=0.05, - initialize_nlte=False): - self.update(t_rad=t_rad, w=ws, j_blues=j_blues) + initialize_nlte=False, previous_iteration_beta_sobolevs=None): + self.update(t_rad=t_rad, w=ws, j_blues=j_blues, + previous_iteration_beta_sobolevs=previous_iteration_beta_sobolevs) def __init__(self, number_densities, atomic_data, time_explosion, t_rad=None, delta_treatment=None, nlte_config=None, @@ -68,9 +72,14 @@ def __init__(self, number_densities, atomic_data, time_explosion, else: raise NotImplementedError('Sorry ' + ionization_mode + ' not implemented yet.') + if nlte_config.species: - raise NotImplementedError('Sorry, NLTE treatment not implemented \ - yet.') + if nlte_config.coronal_approximation == True: + plasma_modules += nlte_coronal_properties + else: + plasma_modules += nlte_general_properties + else: + plasma_modules += non_nlte_properties if line_interaction_type in ('downbranch', 'macroatom'): plasma_modules += macro_atom_properties @@ -87,4 +96,5 @@ def __init__(self, number_densities, atomic_data, time_explosion, t_rad=t_rad, abundance=abundance, density=density, atomic_data=atomic_data, time_explosion=time_explosion, j_blues=None, w=w, link_t_rad_t_electron=link_t_rad_t_electron, - delta_input=delta_treatment) + delta_input=delta_treatment, nlte_species=nlte_config.species, + previous_beta_sobolevs=None) From 4c20e61a910c5c5d88a32e4b25f7dfccfb8e0425 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Tue, 7 Jul 2015 21:34:39 +0100 Subject: [PATCH 02/13] New NLTE system. --- tardis/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/model.py b/tardis/model.py index c9567ebf1e6..c8b13fb473d 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -135,7 +135,7 @@ def __init__(self, tardis_config): line_interaction_type=tardis_config.plasma.line_interaction_type, link_t_rad_t_electron=0.9) - self.previous_iteration_beta_sobolevs = np.ones((len(self.plasma_array.lines), len(self.plasma_array.t_electron))) + self.previous_iteration_beta_sobolevs = np.ones((len(self.plasma_array.lines), len(self.plasma_array.t_rad))) self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_virtual = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) From 4b5782e1d3b3a621cdc14fdd6e779ac2b8695831 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Tue, 7 Jul 2015 22:00:53 +0100 Subject: [PATCH 03/13] New NLTE system. --- tardis/model.py | 4 ++-- tardis/plasma/standard_plasmas.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tardis/model.py b/tardis/model.py index c8b13fb473d..1ad4785155e 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -239,9 +239,9 @@ def update_plasmas(self, initialize_nlte=False): self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, j_blues=self.j_blues, initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05, - previous_iteration_beta_sobolevs = self.previous_iteration_beta_sobolevs) + previous_beta_sobolevs = self.previous_iteration_beta_sobolevs) - self.previous_iteration_beta_sobolevs = self.plasma_array.beta_sobolevs + self.previous_iteration_beta_sobolevs = self.plasma_array.get_value('beta_sobolev') if self.tardis_config.plasma.line_interaction_type in ('downbranch', 'macroatom'): self.transition_probabilities = self.plasma_array.transition_probabilities diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 5ea8c7fe5eb..8734d58023e 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -46,9 +46,9 @@ def initial_w(self, number_densities): def update_radiationfield(self, t_rad, ws, j_blues, t_electrons=None, n_e_convergence_threshold=0.05, - initialize_nlte=False, previous_iteration_beta_sobolevs=None): + initialize_nlte=False, previous_beta_sobolevs=None): self.update(t_rad=t_rad, w=ws, j_blues=j_blues, - previous_iteration_beta_sobolevs=previous_iteration_beta_sobolevs) + previous_beta_sobolevs=previous_beta_sobolevs) def __init__(self, number_densities, atomic_data, time_explosion, t_rad=None, delta_treatment=None, nlte_config=None, From 1aaee4b3e12b87ae334f9d186bc07c11ecd569c1 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Wed, 8 Jul 2015 11:48:40 +0100 Subject: [PATCH 04/13] Inserted fix for lines/nlte_lines indexing discrepancy. --- tardis/plasma/properties/partition_function.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index b058fe1f74a..72818461da1 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -146,6 +146,9 @@ def calculate(t_electron, lines, atomic_data, nlte_data, if len(j_blues)==0: j_blues = lte_j_blues + else: + j_blues = pd.DataFrame(j_blues, index=lines.index, columns = + range(len(t_electron))) for species in nlte_species: logger.info('Calculating rates for species %s', species) @@ -159,6 +162,8 @@ def calculate(t_electron, lines, atomic_data, nlte_data, B_uls = nlte_data.B_uls[species] B_lus = nlte_data.B_lus[species] + j_blues_index = lines.index[lines_index] + r_lu_index = lnu * number_of_levels + lnl r_ul_index = lnl * number_of_levels + lnu @@ -167,14 +172,14 @@ def calculate(t_electron, lines, atomic_data, nlte_data, r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, len(t_electron))) r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ - B_uls[np.newaxis].T * j_blues[lines_index] + B_uls[np.newaxis].T * j_blues.ix[j_blues_index] r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] r_lu_matrix = np.zeros_like(r_ul_matrix) r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, len(t_electron))) r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ - j_blues[lines_index] * beta_sobolevs[lines_index] + j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] # collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ # self.electron_densities.values From 061e75efbbec57f1ed456f59ef7f0f930aa76324 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Wed, 8 Jul 2015 11:55:49 +0100 Subject: [PATCH 05/13] Added classical nebular NLTE. --- .../plasma/properties/partition_function.py | 71 ++++++++++++++++++- .../plasma/properties/property_collections.py | 5 +- tardis/plasma/standard_plasmas.py | 5 +- 3 files changed, 78 insertions(+), 3 deletions(-) diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index 72818461da1..ca42a7fc9e2 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -9,6 +9,7 @@ __all__ = ['LevelBoltzmannFactorLTE', 'LevelBoltzmannFactorDiluteLTE', 'LevelBoltzmannFactorNoNLTE', 'LevelBoltzmannFactorNLTECoronal', + 'LevelBoltzmannFactorNLTEClassicalNebular', 'LevelBoltzmannFactorNLTEGeneral', 'PartitionFunction'] class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): @@ -126,7 +127,75 @@ def calculate(t_electron, lines, atomic_data, nlte_data, return general_level_boltzmann_factor class LevelBoltzmannFactorNLTEClassicalNebular(ProcessingPlasmaProperty): - pass + outputs = ('level_boltzmann_factor',) + + @staticmethod + def calculate(t_electron, lines, atomic_data, nlte_data, + general_level_boltzmann_factor, nlte_species, j_blues, + previous_beta_sobolevs, lte_j_blues): + """ + Calculating the NLTE level populations for specific ions + + """ + beta_sobolevs = np.ones((len(lines), len(t_electron))) + + if len(j_blues)==0: + j_blues = lte_j_blues + else: + j_blues = pd.DataFrame(j_blues, index=lines.index, columns = + range(len(t_electron))) + + for species in nlte_species: + logger.info('Calculating rates for species %s', species) + number_of_levels = atomic_data.levels.energy.ix[species].count() + + lnl = nlte_data.lines_level_number_lower[species] + lnu = nlte_data.lines_level_number_upper[species] + + lines_index = nlte_data.lines_idx[species] + A_uls = nlte_data.A_uls[species] + B_uls = nlte_data.B_uls[species] + B_lus = nlte_data.B_lus[species] + + j_blues_index = lines.index[lines_index] + + r_lu_index = lnu * number_of_levels + lnl + r_ul_index = lnl * number_of_levels + lnu + + r_ul_matrix = np.zeros((number_of_levels, number_of_levels, + len(t_electron)), dtype=np.float64) + r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ + B_uls[np.newaxis].T * j_blues.ix[j_blues_index] + r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] + + r_lu_matrix = np.zeros_like(r_ul_matrix) + r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, + len(t_electron))) + r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ + j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] + +# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ +# self.electron_densities.values + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + + rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix + + for i in xrange(number_of_levels): + rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) + + rates_matrix[0, :, :] = 1.0 + + x = np.zeros(rates_matrix.shape[0]) + x[0] = 1.0 + for i in xrange(len(t_electron)): + level_boltzmann_factor = \ + np.linalg.solve(rates_matrix[:, :, i], x) + general_level_boltzmann_factor[i].ix[species] = \ + level_boltzmann_factor + return general_level_boltzmann_factor class LevelBoltzmannFactorNLTEGeneral(ProcessingPlasmaProperty): outputs = ('level_boltzmann_factor',) diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index 5d63a22b92b..d31f27d3a81 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -10,7 +10,7 @@ RadiationFieldCorrection, RadiationFieldCorrectionInput, LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTECoronal, NLTEData, NLTESpecies, PreviousIterationBetaSobolevs, LTEJBlues, - LevelBoltzmannFactorNLTEGeneral) + LevelBoltzmannFactorNLTEGeneral, LevelBoltzmannFactorNLTEClassicalNebular) class PlasmaPropertyCollection(list): pass @@ -34,6 +34,9 @@ class PlasmaPropertyCollection(list): non_nlte_properties = PlasmaPropertyCollection([LevelBoltzmannFactorNoNLTE]) nlte_coronal_properties = PlasmaPropertyCollection([ LevelBoltzmannFactorNLTECoronal, NLTEData, NLTESpecies]) +nlte_classical_nebular_properties = PlasmaPropertyCollection([ + LevelBoltzmannFactorNLTEClassicalNebular, NLTEData, NLTESpecies, + LTEJBlues]) nlte_general_properties = PlasmaPropertyCollection([ LevelBoltzmannFactorNLTEGeneral, NLTEData, NLTESpecies, PreviousIterationBetaSobolevs, LTEJBlues]) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 8734d58023e..7665102da4b 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -7,7 +7,8 @@ basic_properties, lte_excitation_properties, lte_ionization_properties, macro_atom_properties, dilute_lte_excitation_properties, nebular_ionization_properties, non_nlte_properties, - nlte_coronal_properties, nlte_general_properties) + nlte_coronal_properties, nlte_general_properties, + nlte_classical_nebular_properties) logger = logging.getLogger(__name__) @@ -76,6 +77,8 @@ def __init__(self, number_densities, atomic_data, time_explosion, if nlte_config.species: if nlte_config.coronal_approximation == True: plasma_modules += nlte_coronal_properties + elif nlte_config.classical_nebular == True: + plasma_modules += nlte_classical_nebular_properties else: plasma_modules += nlte_general_properties else: From 992778492c1014616d6465fa55c5f485dd24979f Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Wed, 8 Jul 2015 12:17:26 +0100 Subject: [PATCH 06/13] Re-added collision matrix. --- tardis/model.py | 5 +- .../plasma/properties/partition_function.py | 49 +++++++++++++------ tardis/plasma/properties/plasma_input.py | 6 ++- .../plasma/properties/property_collections.py | 9 ++-- tardis/plasma/standard_plasmas.py | 7 +-- 5 files changed, 52 insertions(+), 24 deletions(-) diff --git a/tardis/model.py b/tardis/model.py index 1ad4785155e..d5e8faf1f80 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -136,6 +136,7 @@ def __init__(self, tardis_config): link_t_rad_t_electron=0.9) self.previous_iteration_beta_sobolevs = np.ones((len(self.plasma_array.lines), len(self.plasma_array.t_rad))) + self.previous_iteration_electron_densities = self.plasma_array.get_value('electron_densities') self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_virtual = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) @@ -239,9 +240,11 @@ def update_plasmas(self, initialize_nlte=False): self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, j_blues=self.j_blues, initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05, - previous_beta_sobolevs = self.previous_iteration_beta_sobolevs) + previous_beta_sobolevs = self.previous_iteration_beta_sobolevs, + previous_electron_densities = self.previous_iteration_electron_densities) self.previous_iteration_beta_sobolevs = self.plasma_array.get_value('beta_sobolev') + self.previous_iteration_electron_densities = self.plasma_array.get_value('electron_densities') if self.tardis_config.plasma.line_interaction_type in ('downbranch', 'macroatom'): self.transition_probabilities = self.plasma_array.transition_probabilities diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index ca42a7fc9e2..dd1f6f16ade 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -67,7 +67,8 @@ class LevelBoltzmannFactorNLTECoronal(ProcessingPlasmaProperty): @staticmethod def calculate(t_electron, lines, atomic_data, nlte_data, - general_level_boltzmann_factor, nlte_species): + general_level_boltzmann_factor, nlte_species, + previous_electron_densities): """ Calculating the NLTE level populations for specific ions @@ -105,10 +106,16 @@ def calculate(t_electron, lines, atomic_data, nlte_data, r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ j_blues[lines_index] * beta_sobolevs[lines_index] -# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ -# self.electron_densities.values - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) + if atomic_data.has_collision_data: + if previous_electron_densities is None: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + else: + collision_matrix = nlte_data.get_collision_matrix(species, + t_electron) * previous_electron_densities.values + else: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix @@ -132,7 +139,7 @@ class LevelBoltzmannFactorNLTEClassicalNebular(ProcessingPlasmaProperty): @staticmethod def calculate(t_electron, lines, atomic_data, nlte_data, general_level_boltzmann_factor, nlte_species, j_blues, - previous_beta_sobolevs, lte_j_blues): + previous_beta_sobolevs, lte_j_blues, previous_electron_densities): """ Calculating the NLTE level populations for specific ions @@ -176,10 +183,16 @@ def calculate(t_electron, lines, atomic_data, nlte_data, r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] -# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ -# self.electron_densities.values - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) + if atomic_data.has_collision_data: + if previous_electron_densities is None: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + else: + collision_matrix = nlte_data.get_collision_matrix(species, + t_electron) * previous_electron_densities.values + else: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix @@ -203,7 +216,7 @@ class LevelBoltzmannFactorNLTEGeneral(ProcessingPlasmaProperty): @staticmethod def calculate(t_electron, lines, atomic_data, nlte_data, general_level_boltzmann_factor, nlte_species, j_blues, - previous_beta_sobolevs, lte_j_blues): + previous_beta_sobolevs, lte_j_blues, previous_electron_densities): """ Calculating the NLTE level populations for specific ions @@ -250,10 +263,16 @@ def calculate(t_electron, lines, atomic_data, nlte_data, r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] -# collision_matrix = self.atom_data.nlte_data.get_collision_matrix(species, self.t_electrons) * \ -# self.electron_densities.values - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) + if atomic_data.has_collision_data: + if previous_electron_densities is None: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) + else: + collision_matrix = nlte_data.get_collision_matrix(species, + t_electron) * previous_electron_densities.values + else: + collision_matrix = r_ul_matrix.copy() + collision_matrix.fill(0.0) rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix diff --git a/tardis/plasma/properties/plasma_input.py b/tardis/plasma/properties/plasma_input.py index 17fd765e2bd..74adbfb6f6e 100644 --- a/tardis/plasma/properties/plasma_input.py +++ b/tardis/plasma/properties/plasma_input.py @@ -5,7 +5,8 @@ __all__ = ['TRadiative', 'DilutionFactor', 'AtomicData', 'Abundance', 'Density', 'TimeExplosion', 'JBlues', 'LinkTRadTElectron', 'RadiationFieldCorrectionInput', 'NLTESpecies', - 'PreviousIterationBetaSobolevs'] + 'PreviousIterationBetaSobolevs', + 'PreviousIterationElectronDensities'] class Input(BasePlasmaProperty): @@ -72,3 +73,6 @@ class NLTESpecies(StaticInput): class PreviousIterationBetaSobolevs(DynamicInput): outputs = ('previous_beta_sobolevs',) + +class PreviousIterationElectronDensities(DynamicInput): + outputs = ('previous_electron_densities',) diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index d31f27d3a81..1e7768cc426 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -10,14 +10,16 @@ RadiationFieldCorrection, RadiationFieldCorrectionInput, LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTECoronal, NLTEData, NLTESpecies, PreviousIterationBetaSobolevs, LTEJBlues, - LevelBoltzmannFactorNLTEGeneral, LevelBoltzmannFactorNLTEClassicalNebular) + LevelBoltzmannFactorNLTEGeneral, LevelBoltzmannFactorNLTEClassicalNebular, + PreviousIterationElectronDensities) class PlasmaPropertyCollection(list): pass basic_inputs = PlasmaPropertyCollection([TRadiative, Abundance, Density, TimeExplosion, AtomicData, JBlues, DilutionFactor, LinkTRadTElectron, - RadiationFieldCorrectionInput, NLTESpecies, PreviousIterationBetaSobolevs]) + RadiationFieldCorrectionInput, NLTESpecies, PreviousIterationBetaSobolevs, + PreviousIterationElectronDensities]) basic_properties = PlasmaPropertyCollection([BetaRadiation, Levels, Lines, AtomicMass, LevelPopulation, PartitionFunction, GElectron, IonizationData, NumberDensity, LinesLowerLevelIndex, @@ -38,5 +40,4 @@ class PlasmaPropertyCollection(list): LevelBoltzmannFactorNLTEClassicalNebular, NLTEData, NLTESpecies, LTEJBlues]) nlte_general_properties = PlasmaPropertyCollection([ - LevelBoltzmannFactorNLTEGeneral, NLTEData, NLTESpecies, - PreviousIterationBetaSobolevs, LTEJBlues]) + LevelBoltzmannFactorNLTEGeneral, NLTEData, NLTESpecies, LTEJBlues]) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 7665102da4b..d34a925cd85 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -25,7 +25,7 @@ def __init__(self, t_rad, abundance, density, time_explosion, atomic_data, density=density, time_explosion=time_explosion, j_blues=j_blues, w=None, link_t_rad_t_electron=link_t_rad_t_electron, delta_input=delta_treatment, nlte_species=None, - previous_beta_sobolevs=None) + previous_beta_sobolevs=None, previous_electron_densities=None) class LegacyPlasmaArray(BasePlasma): @@ -49,7 +49,8 @@ def update_radiationfield(self, t_rad, ws, j_blues, t_electrons=None, n_e_convergence_threshold=0.05, initialize_nlte=False, previous_beta_sobolevs=None): self.update(t_rad=t_rad, w=ws, j_blues=j_blues, - previous_beta_sobolevs=previous_beta_sobolevs) + previous_beta_sobolevs=previous_beta_sobolevs, + previous_electron_densities=previous_electron_densities) def __init__(self, number_densities, atomic_data, time_explosion, t_rad=None, delta_treatment=None, nlte_config=None, @@ -100,4 +101,4 @@ def __init__(self, number_densities, atomic_data, time_explosion, atomic_data=atomic_data, time_explosion=time_explosion, j_blues=None, w=w, link_t_rad_t_electron=link_t_rad_t_electron, delta_input=delta_treatment, nlte_species=nlte_config.species, - previous_beta_sobolevs=None) + previous_beta_sobolevs=None, previous_electron_densities=None) From 8d9a0736b0ee829755d8d1b08001714bb9c3f710 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Wed, 8 Jul 2015 12:24:32 +0100 Subject: [PATCH 07/13] Fix. --- tardis/plasma/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tardis/plasma/base.py b/tardis/plasma/base.py index e52abeaf1ec..5e6a433fd86 100644 --- a/tardis/plasma/base.py +++ b/tardis/plasma/base.py @@ -206,6 +206,7 @@ def __init__(self, number_densities, atom_data, time_explosion, delta_treatment=None, nlte_config=None, ionization_mode='lte', excitation_mode='lte', w=None, link_t_rad_t_electron=0.9, nlte_species=None, - previous_beta_sobolevs=None): + previous_beta_sobolevs=None, + previous_electron_densities=None): pass From 7aa5c03354b1a8963de1fb879ac9faa3ec3cb44e Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Wed, 8 Jul 2015 12:30:14 +0100 Subject: [PATCH 08/13] Fix. --- tardis/plasma/standard_plasmas.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index d34a925cd85..5ce11925127 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -47,7 +47,8 @@ def initial_w(self, number_densities): def update_radiationfield(self, t_rad, ws, j_blues, t_electrons=None, n_e_convergence_threshold=0.05, - initialize_nlte=False, previous_beta_sobolevs=None): + initialize_nlte=False, previous_beta_sobolevs=None, + previous_electron_densities=None): self.update(t_rad=t_rad, w=ws, j_blues=j_blues, previous_beta_sobolevs=previous_beta_sobolevs, previous_electron_densities=previous_electron_densities) From 1070c25be6ebbdaf7b71b74450e0dad680ccab9c Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Thu, 30 Jul 2015 13:27:34 +0100 Subject: [PATCH 09/13] WIP: Will fail tests. Change to NLTE structure. --- tardis/model.py | 12 ++---------- tardis/plasma/base.py | 6 ++++++ tardis/plasma/properties/__init__.py | 1 + tardis/plasma/properties/nlte.py | 17 +++++++++++++++++ .../plasma/properties/property_collections.py | 8 ++++---- tardis/plasma/standard_plasmas.py | 18 +++++++++++------- 6 files changed, 41 insertions(+), 21 deletions(-) create mode 100644 tardis/plasma/properties/nlte.py diff --git a/tardis/model.py b/tardis/model.py index d5e8faf1f80..078fe14eee4 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -135,9 +135,6 @@ def __init__(self, tardis_config): line_interaction_type=tardis_config.plasma.line_interaction_type, link_t_rad_t_electron=0.9) - self.previous_iteration_beta_sobolevs = np.ones((len(self.plasma_array.lines), len(self.plasma_array.t_rad))) - self.previous_iteration_electron_densities = self.plasma_array.get_value('electron_densities') - self.spectrum = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_virtual = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) self.spectrum_reabsorbed = TARDISSpectrum(tardis_config.spectrum.frequency, tardis_config.supernova.distance) @@ -238,13 +235,8 @@ def calculate_updated_radiationfield(self, nubar_estimator, j_estimator): def update_plasmas(self, initialize_nlte=False): - self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, j_blues=self.j_blues, - initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05, - previous_beta_sobolevs = self.previous_iteration_beta_sobolevs, - previous_electron_densities = self.previous_iteration_electron_densities) - - self.previous_iteration_beta_sobolevs = self.plasma_array.get_value('beta_sobolev') - self.previous_iteration_electron_densities = self.plasma_array.get_value('electron_densities') + self.plasma_array.update_radiationfield(self.t_rads.value, self.ws, self.j_blues, + self.tardis_config.plasma.nlte, initialize_nlte=initialize_nlte, n_e_convergence_threshold=0.05) if self.tardis_config.plasma.line_interaction_type in ('downbranch', 'macroatom'): self.transition_probabilities = self.plasma_array.transition_probabilities diff --git a/tardis/plasma/base.py b/tardis/plasma/base.py index 5e6a433fd86..38ac549f0e4 100644 --- a/tardis/plasma/base.py +++ b/tardis/plasma/base.py @@ -119,6 +119,12 @@ def _init_properties(self, plasma_properties, **kwargs): plasma_property_objects.append(current_property_object) return plasma_property_objects + def store_previous_properties(self): + self.outputs_dict['previous_electron_densities'].set_value( + self.get_value('electron_densities')) + self.outputs_dict['previous_beta_sobolevs'].set_value( + self.get_value('beta_sobolev')) + def update(self, **kwargs): for key in kwargs: if key not in self.outputs_dict: diff --git a/tardis/plasma/properties/__init__.py b/tardis/plasma/properties/__init__.py index fc136911b8b..458c59f6a93 100644 --- a/tardis/plasma/properties/__init__.py +++ b/tardis/plasma/properties/__init__.py @@ -5,4 +5,5 @@ from tardis.plasma.properties.partition_function import * from tardis.plasma.properties.plasma_input import * from tardis.plasma.properties.radiative_properties import * +from tardis.plasma.properties.nlte import * diff --git a/tardis/plasma/properties/nlte.py b/tardis/plasma/properties/nlte.py new file mode 100644 index 00000000000..d0246fdcf4e --- /dev/null +++ b/tardis/plasma/properties/nlte.py @@ -0,0 +1,17 @@ +from tardis.plasma.properties.base import BasePlasmaProperty + +__all__ = ['PreviousElectronDensities', 'PreviousBetaSobolevs'] + +class PreviousIterationProperty(BasePlasmaProperty): + def _set_output_value(self, output, value): + setattr(self, output, value) + + def set_value(self, value): + assert len(self.outputs) == 1 + self._set_output_value(self.outputs[0], value) + +class PreviousElectronDensities(PreviousIterationProperty): + outputs = ('previous_electron_densities',) + +class PreviousBetaSobolevs(PreviousIterationProperty): + outputs = ('previous_beta_sobolevs',) \ No newline at end of file diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index 1e7768cc426..ebd939314dc 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -9,17 +9,17 @@ ZetaData, ElectronTemperature, LinkTRadTElectron, BetaElectron, RadiationFieldCorrection, RadiationFieldCorrectionInput, LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTECoronal, NLTEData, - NLTESpecies, PreviousIterationBetaSobolevs, LTEJBlues, + NLTESpecies, PreviousBetaSobolevs, LTEJBlues, LevelBoltzmannFactorNLTEGeneral, LevelBoltzmannFactorNLTEClassicalNebular, - PreviousIterationElectronDensities) + PreviousElectronDensities) class PlasmaPropertyCollection(list): pass basic_inputs = PlasmaPropertyCollection([TRadiative, Abundance, Density, TimeExplosion, AtomicData, JBlues, DilutionFactor, LinkTRadTElectron, - RadiationFieldCorrectionInput, NLTESpecies, PreviousIterationBetaSobolevs, - PreviousIterationElectronDensities]) + RadiationFieldCorrectionInput, NLTESpecies, PreviousBetaSobolevs, + PreviousElectronDensities]) basic_properties = PlasmaPropertyCollection([BetaRadiation, Levels, Lines, AtomicMass, LevelPopulation, PartitionFunction, GElectron, IonizationData, NumberDensity, LinesLowerLevelIndex, diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index 5ce11925127..cc4592fb743 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -45,13 +45,12 @@ def initial_t_rad(self, number_densities): def initial_w(self, number_densities): return np.ones(len(number_densities.columns)) * 0.5 - def update_radiationfield(self, t_rad, ws, j_blues, + def update_radiationfield(self, t_rad, ws, j_blues, nlte_config, t_electrons=None, n_e_convergence_threshold=0.05, - initialize_nlte=False, previous_beta_sobolevs=None, - previous_electron_densities=None): - self.update(t_rad=t_rad, w=ws, j_blues=j_blues, - previous_beta_sobolevs=previous_beta_sobolevs, - previous_electron_densities=previous_electron_densities) + initialize_nlte=False): + if nlte_config.species: + self.store_previous_properties() + self.update(t_rad=t_rad, w=ws, j_blues=j_blues) def __init__(self, number_densities, atomic_data, time_explosion, t_rad=None, delta_treatment=None, nlte_config=None, @@ -97,9 +96,14 @@ def __init__(self, number_densities, atomic_data, time_explosion, abundance, density = self.from_number_densities(number_densities, atomic_data) + initial_beta_sobolevs = np.ones((len(atomic_data.lines), + len(number_densities.columns))) + initial_electron_densities = number_densities.sum(axis=0) + super(LegacyPlasmaArray, self).__init__(plasma_properties=plasma_modules, t_rad=t_rad, abundance=abundance, density=density, atomic_data=atomic_data, time_explosion=time_explosion, j_blues=None, w=w, link_t_rad_t_electron=link_t_rad_t_electron, delta_input=delta_treatment, nlte_species=nlte_config.species, - previous_beta_sobolevs=None, previous_electron_densities=None) + previous_electron_densities=initial_electron_densities, + previous_beta_sobolevs=initial_beta_sobolevs) From ccc78d6109e7de84ff7ff62a5472835a6fb0e3c1 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Fri, 31 Jul 2015 00:02:45 +0100 Subject: [PATCH 10/13] Changing NLTE. --- .../plasma/properties/partition_function.py | 40 +++++++++++++------ .../plasma/properties/property_collections.py | 12 ++---- tardis/plasma/standard_plasmas.py | 12 ++---- 3 files changed, 34 insertions(+), 30 deletions(-) diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index dd1f6f16ade..982201d4cdc 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -8,17 +8,14 @@ logger = logging.getLogger(__name__) __all__ = ['LevelBoltzmannFactorLTE', 'LevelBoltzmannFactorDiluteLTE', - 'LevelBoltzmannFactorNoNLTE', 'LevelBoltzmannFactorNLTECoronal', - 'LevelBoltzmannFactorNLTEClassicalNebular', - 'LevelBoltzmannFactorNLTEGeneral', 'PartitionFunction'] + 'LevelBoltzmannFactorNoNLTE', 'LevelBoltzmannFactorNLTE', + 'PartitionFunction'] class LevelBoltzmannFactorLTE(ProcessingPlasmaProperty): """ Calculate the level population Boltzmann factor - .. math: {latex_formula} - """ outputs = ('general_level_boltzmann_factor',) @@ -60,7 +57,7 @@ class LevelBoltzmannFactorNoNLTE(ProcessingPlasmaProperty): @staticmethod def calculate(general_level_boltzmann_factor): return general_level_boltzmann_factor - +''' class LevelBoltzmannFactorNLTECoronal(ProcessingPlasmaProperty): outputs = ('level_boltzmann_factor',) @@ -71,7 +68,6 @@ def calculate(t_electron, lines, atomic_data, nlte_data, previous_electron_densities): """ Calculating the NLTE level populations for specific ions - """ beta_sobolevs = np.ones((len(lines), len(t_electron))) @@ -142,7 +138,6 @@ def calculate(t_electron, lines, atomic_data, nlte_data, previous_beta_sobolevs, lte_j_blues, previous_electron_densities): """ Calculating the NLTE level populations for specific ions - """ beta_sobolevs = np.ones((len(lines), len(t_electron))) @@ -209,17 +204,36 @@ def calculate(t_electron, lines, atomic_data, nlte_data, general_level_boltzmann_factor[i].ix[species] = \ level_boltzmann_factor return general_level_boltzmann_factor - -class LevelBoltzmannFactorNLTEGeneral(ProcessingPlasmaProperty): +''' +class LevelBoltzmannFactorNLTE(ProcessingPlasmaProperty): outputs = ('level_boltzmann_factor',) + def __init__(self, plasma_parent, classical_nebular=False, + coronal_approximation=False): + super(LevelBoltzmannFactorNLTE, self).__init__(plasma_parent) + if classical_nebular == True: + self.calculate = self._calculate_classical_nebular + elif coronal_approximation == True: + self.calculate = self._calculate_coronal_approximation + else: + self.calculate = self._calculate_general + self._update_inputs() + + def _calculate_classical_nebular(self, general_level_boltzmann_factor): + return general_level_boltzmann_factor + def _calculate_coronal_approximation(self, general_level_boltzmann_factor): + return general_level_boltzmann_factor + def _calculate_general(self, general_level_boltzmann_factor): + return general_level_boltzmann_factor + def _main_nlte_calculation(self): + pass +''' @staticmethod def calculate(t_electron, lines, atomic_data, nlte_data, general_level_boltzmann_factor, nlte_species, j_blues, previous_beta_sobolevs, lte_j_blues, previous_electron_densities): """ Calculating the NLTE level populations for specific ions - """ if previous_beta_sobolevs is None: beta_sobolevs = np.ones((len(lines), len(t_electron))) @@ -289,7 +303,7 @@ def calculate(t_electron, lines, atomic_data, nlte_data, general_level_boltzmann_factor[i].ix[species] = \ level_boltzmann_factor return general_level_boltzmann_factor - +''' class PartitionFunction(ProcessingPlasmaProperty): outputs = ('partition_function',) latex_outputs = '$Z_{i, j}$' @@ -300,4 +314,4 @@ class PartitionFunction(ProcessingPlasmaProperty): @staticmethod def calculate(levels, level_boltzmann_factor): return level_boltzmann_factor.groupby( - level=['atomic_number', 'ion_number']).sum() + level=['atomic_number', 'ion_number']).sum() \ No newline at end of file diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py index ebd939314dc..37e4cb5a5e5 100644 --- a/tardis/plasma/properties/property_collections.py +++ b/tardis/plasma/properties/property_collections.py @@ -8,9 +8,8 @@ PhiGeneral, PhiSahaNebular, LevelBoltzmannFactorDiluteLTE, DilutionFactor, ZetaData, ElectronTemperature, LinkTRadTElectron, BetaElectron, RadiationFieldCorrection, RadiationFieldCorrectionInput, - LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTECoronal, NLTEData, + LevelBoltzmannFactorNoNLTE, LevelBoltzmannFactorNLTE, NLTEData, NLTESpecies, PreviousBetaSobolevs, LTEJBlues, - LevelBoltzmannFactorNLTEGeneral, LevelBoltzmannFactorNLTEClassicalNebular, PreviousElectronDensities) class PlasmaPropertyCollection(list): @@ -34,10 +33,5 @@ class PlasmaPropertyCollection(list): dilute_lte_excitation_properties = PlasmaPropertyCollection([ LevelBoltzmannFactorDiluteLTE]) non_nlte_properties = PlasmaPropertyCollection([LevelBoltzmannFactorNoNLTE]) -nlte_coronal_properties = PlasmaPropertyCollection([ - LevelBoltzmannFactorNLTECoronal, NLTEData, NLTESpecies]) -nlte_classical_nebular_properties = PlasmaPropertyCollection([ - LevelBoltzmannFactorNLTEClassicalNebular, NLTEData, NLTESpecies, - LTEJBlues]) -nlte_general_properties = PlasmaPropertyCollection([ - LevelBoltzmannFactorNLTEGeneral, NLTEData, NLTESpecies, LTEJBlues]) +nlte_properties = PlasmaPropertyCollection([ + LevelBoltzmannFactorNLTE, NLTEData, NLTESpecies, LTEJBlues]) diff --git a/tardis/plasma/standard_plasmas.py b/tardis/plasma/standard_plasmas.py index cc4592fb743..15f2b567ae3 100644 --- a/tardis/plasma/standard_plasmas.py +++ b/tardis/plasma/standard_plasmas.py @@ -7,8 +7,7 @@ basic_properties, lte_excitation_properties, lte_ionization_properties, macro_atom_properties, dilute_lte_excitation_properties, nebular_ionization_properties, non_nlte_properties, - nlte_coronal_properties, nlte_general_properties, - nlte_classical_nebular_properties) + nlte_properties) logger = logging.getLogger(__name__) @@ -76,12 +75,7 @@ def __init__(self, number_densities, atomic_data, time_explosion, ' not implemented yet.') if nlte_config.species: - if nlte_config.coronal_approximation == True: - plasma_modules += nlte_coronal_properties - elif nlte_config.classical_nebular == True: - plasma_modules += nlte_classical_nebular_properties - else: - plasma_modules += nlte_general_properties + plasma_modules += nlte_properties else: plasma_modules += non_nlte_properties @@ -100,6 +94,8 @@ def __init__(self, number_densities, atomic_data, time_explosion, len(number_densities.columns))) initial_electron_densities = number_densities.sum(axis=0) + self.nlte_config = nlte_config + super(LegacyPlasmaArray, self).__init__(plasma_properties=plasma_modules, t_rad=t_rad, abundance=abundance, density=density, atomic_data=atomic_data, time_explosion=time_explosion, From d92cb1d3e006707c047bb9e87ef00e97bb136306 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Fri, 31 Jul 2015 10:33:00 +0100 Subject: [PATCH 11/13] Combined NLTE options into single property. --- .../plasma/properties/partition_function.py | 222 ++++-------------- 1 file changed, 43 insertions(+), 179 deletions(-) diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index 982201d4cdc..a84818a8e0d 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -57,51 +57,50 @@ class LevelBoltzmannFactorNoNLTE(ProcessingPlasmaProperty): @staticmethod def calculate(general_level_boltzmann_factor): return general_level_boltzmann_factor -''' -class LevelBoltzmannFactorNLTECoronal(ProcessingPlasmaProperty): +class LevelBoltzmannFactorNLTE(ProcessingPlasmaProperty): outputs = ('level_boltzmann_factor',) - @staticmethod - def calculate(t_electron, lines, atomic_data, nlte_data, - general_level_boltzmann_factor, nlte_species, - previous_electron_densities): - """ - Calculating the NLTE level populations for specific ions - """ + def calculate(self): + pass - beta_sobolevs = np.ones((len(lines), len(t_electron))) - j_blues = np.zeros((len(lines), len(t_electron))) + def __init__(self, plasma_parent, classical_nebular=False, + coronal_approximation=False): + super(LevelBoltzmannFactorNLTE, self).__init__(plasma_parent) + if classical_nebular == True: + self.calculate = self._calculate_classical_nebular + elif coronal_approximation == True: + self.calculate = self._calculate_coronal_approximation + else: + self.calculate = self._calculate_general + self._update_inputs() + def _main_nlte_calculation(self, nlte_species, atomic_data, nlte_data, + t_electron, j_blues, beta_sobolevs, general_level_boltzmann_factor, + previous_electron_densities): for species in nlte_species: logger.info('Calculating rates for species %s', species) number_of_levels = atomic_data.levels.energy.ix[species].count() - lnl = nlte_data.lines_level_number_lower[species] lnu = nlte_data.lines_level_number_upper[species] - lines_index = nlte_data.lines_idx[species] A_uls = nlte_data.A_uls[species] B_uls = nlte_data.B_uls[species] B_lus = nlte_data.B_lus[species] - r_lu_index = lnu * number_of_levels + lnl r_ul_index = lnl * number_of_levels + lnu - r_ul_matrix = np.zeros((number_of_levels, number_of_levels, len(t_electron)), dtype=np.float64) r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, len(t_electron))) r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ - B_uls[np.newaxis].T * j_blues[lines_index] + B_uls[np.newaxis].T * j_blues.ix[lines_index] r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] - r_lu_matrix = np.zeros_like(r_ul_matrix) r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, len(t_electron))) r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ - j_blues[lines_index] * beta_sobolevs[lines_index] - + j_blues.ix[lines_index] * beta_sobolevs[lines_index] if atomic_data.has_collision_data: if previous_electron_densities is None: collision_matrix = r_ul_matrix.copy() @@ -112,14 +111,10 @@ def calculate(t_electron, lines, atomic_data, nlte_data, else: collision_matrix = r_ul_matrix.copy() collision_matrix.fill(0.0) - rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix - for i in xrange(number_of_levels): rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) - rates_matrix[0, :, :] = 1.0 - x = np.zeros(rates_matrix.shape[0]) x[0] = 1.0 for i in xrange(len(t_electron)): @@ -129,181 +124,50 @@ def calculate(t_electron, lines, atomic_data, nlte_data, level_boltzmann_factor return general_level_boltzmann_factor -class LevelBoltzmannFactorNLTEClassicalNebular(ProcessingPlasmaProperty): - outputs = ('level_boltzmann_factor',) - - @staticmethod - def calculate(t_electron, lines, atomic_data, nlte_data, - general_level_boltzmann_factor, nlte_species, j_blues, - previous_beta_sobolevs, lte_j_blues, previous_electron_densities): - """ - Calculating the NLTE level populations for specific ions - """ + def _calculate_classical_nebular(self, t_electron, lines, atomic_data, + nlte_data, general_level_boltzmann_factor, nlte_species, j_blues, + previous_beta_sobolevs, lte_j_blues, previous_electron_densities): beta_sobolevs = np.ones((len(lines), len(t_electron))) - if len(j_blues)==0: j_blues = lte_j_blues else: j_blues = pd.DataFrame(j_blues, index=lines.index, columns = range(len(t_electron))) - - for species in nlte_species: - logger.info('Calculating rates for species %s', species) - number_of_levels = atomic_data.levels.energy.ix[species].count() - - lnl = nlte_data.lines_level_number_lower[species] - lnu = nlte_data.lines_level_number_upper[species] - - lines_index = nlte_data.lines_idx[species] - A_uls = nlte_data.A_uls[species] - B_uls = nlte_data.B_uls[species] - B_lus = nlte_data.B_lus[species] - - j_blues_index = lines.index[lines_index] - - r_lu_index = lnu * number_of_levels + lnl - r_ul_index = lnl * number_of_levels + lnu - - r_ul_matrix = np.zeros((number_of_levels, number_of_levels, - len(t_electron)), dtype=np.float64) - r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, - len(t_electron))) - r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ - B_uls[np.newaxis].T * j_blues.ix[j_blues_index] - r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] - - r_lu_matrix = np.zeros_like(r_ul_matrix) - r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, - len(t_electron))) - r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ - j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] - - if atomic_data.has_collision_data: - if previous_electron_densities is None: - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) - else: - collision_matrix = nlte_data.get_collision_matrix(species, - t_electron) * previous_electron_densities.values - else: - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) - - rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix - - for i in xrange(number_of_levels): - rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) - - rates_matrix[0, :, :] = 1.0 - - x = np.zeros(rates_matrix.shape[0]) - x[0] = 1.0 - for i in xrange(len(t_electron)): - level_boltzmann_factor = \ - np.linalg.solve(rates_matrix[:, :, i], x) - general_level_boltzmann_factor[i].ix[species] = \ - level_boltzmann_factor + general_level_boltzmann_factor = self._main_nlte_calculation( + nlte_species, atomic_data, nlte_data, t_electron, j_blues, + beta_sobolevs, general_level_boltzmann_factor, + previous_electron_densities) return general_level_boltzmann_factor -''' -class LevelBoltzmannFactorNLTE(ProcessingPlasmaProperty): - outputs = ('level_boltzmann_factor',) - - def __init__(self, plasma_parent, classical_nebular=False, - coronal_approximation=False): - super(LevelBoltzmannFactorNLTE, self).__init__(plasma_parent) - if classical_nebular == True: - self.calculate = self._calculate_classical_nebular - elif coronal_approximation == True: - self.calculate = self._calculate_coronal_approximation - else: - self.calculate = self._calculate_general - self._update_inputs() - def _calculate_classical_nebular(self, general_level_boltzmann_factor): - return general_level_boltzmann_factor - def _calculate_coronal_approximation(self, general_level_boltzmann_factor): - return general_level_boltzmann_factor - def _calculate_general(self, general_level_boltzmann_factor): + def _calculate_coronal_approximation(self, t_electron, lines, atomic_data, + nlte_data, general_level_boltzmann_factor, nlte_species, + previous_electron_densities): + beta_sobolevs = np.ones((len(lines), len(t_electron))) + j_blues = np.zeros((len(lines), len(t_electron))) + general_level_boltzmann_factor = self._main_nlte_calculation( + nlte_species, atomic_data, nlte_data, t_electron, j_blues, + beta_sobolevs, general_level_boltzmann_factor, + previous_electron_densities) return general_level_boltzmann_factor - def _main_nlte_calculation(self): - pass -''' - @staticmethod - def calculate(t_electron, lines, atomic_data, nlte_data, - general_level_boltzmann_factor, nlte_species, j_blues, - previous_beta_sobolevs, lte_j_blues, previous_electron_densities): - """ - Calculating the NLTE level populations for specific ions - """ + + def _calculate_general(self, t_electron, lines, atomic_data, nlte_data, + general_level_boltzmann_factor, nlte_species, j_blues, + previous_beta_sobolevs, lte_j_blues, previous_electron_densities): if previous_beta_sobolevs is None: beta_sobolevs = np.ones((len(lines), len(t_electron))) else: beta_sobolevs = previous_beta_sobolevs - if len(j_blues)==0: j_blues = lte_j_blues else: j_blues = pd.DataFrame(j_blues, index=lines.index, columns = range(len(t_electron))) - - for species in nlte_species: - logger.info('Calculating rates for species %s', species) - number_of_levels = atomic_data.levels.energy.ix[species].count() - - lnl = nlte_data.lines_level_number_lower[species] - lnu = nlte_data.lines_level_number_upper[species] - - lines_index = nlte_data.lines_idx[species] - A_uls = nlte_data.A_uls[species] - B_uls = nlte_data.B_uls[species] - B_lus = nlte_data.B_lus[species] - - j_blues_index = lines.index[lines_index] - - r_lu_index = lnu * number_of_levels + lnl - r_ul_index = lnl * number_of_levels + lnu - - r_ul_matrix = np.zeros((number_of_levels, number_of_levels, - len(t_electron)), dtype=np.float64) - r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, - len(t_electron))) - r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ - B_uls[np.newaxis].T * j_blues.ix[j_blues_index] - r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] - - r_lu_matrix = np.zeros_like(r_ul_matrix) - r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, - len(t_electron))) - r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ - j_blues.ix[j_blues_index] * beta_sobolevs[lines_index] - - if atomic_data.has_collision_data: - if previous_electron_densities is None: - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) - else: - collision_matrix = nlte_data.get_collision_matrix(species, - t_electron) * previous_electron_densities.values - else: - collision_matrix = r_ul_matrix.copy() - collision_matrix.fill(0.0) - - rates_matrix = r_lu_matrix + r_ul_matrix + collision_matrix - - for i in xrange(number_of_levels): - rates_matrix[i, i] = -rates_matrix[:, i].sum(axis=0) - - rates_matrix[0, :, :] = 1.0 - - x = np.zeros(rates_matrix.shape[0]) - x[0] = 1.0 - for i in xrange(len(t_electron)): - level_boltzmann_factor = \ - np.linalg.solve(rates_matrix[:, :, i], x) - general_level_boltzmann_factor[i].ix[species] = \ - level_boltzmann_factor + general_level_boltzmann_factor = self._main_nlte_calculation( + nlte_species, atomic_data, nlte_data, t_electron, j_blues, + beta_sobolevs, general_level_boltzmann_factor, + previous_electron_densities) return general_level_boltzmann_factor -''' + class PartitionFunction(ProcessingPlasmaProperty): outputs = ('partition_function',) latex_outputs = '$Z_{i, j}$' From affeab6c3b1b5383bcddbe0f676de0f92a9558a8 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Fri, 31 Jul 2015 16:02:05 +0100 Subject: [PATCH 12/13] Preventing nlte config contradiction. --- tardis/plasma/exceptions.py | 3 +++ tardis/plasma/properties/partition_function.py | 7 +++++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/tardis/plasma/exceptions.py b/tardis/plasma/exceptions.py index c450fa4755a..a0a1b728d8f 100644 --- a/tardis/plasma/exceptions.py +++ b/tardis/plasma/exceptions.py @@ -19,4 +19,7 @@ class NotInitializedModule(PlasmaException): pass class PlasmaIonizationError(PlasmaException): + pass + +class PlasmaConfigContradiction(PlasmaException): pass \ No newline at end of file diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index a84818a8e0d..2bcda0cf947 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -4,6 +4,7 @@ import pandas as pd from tardis.plasma.properties.base import ProcessingPlasmaProperty +from tardis.plasma.exceptions import PlasmaConfigContradiction logger = logging.getLogger(__name__) @@ -67,10 +68,12 @@ def calculate(self): def __init__(self, plasma_parent, classical_nebular=False, coronal_approximation=False): super(LevelBoltzmannFactorNLTE, self).__init__(plasma_parent) - if classical_nebular == True: + if classical_nebular == True and coronal_approximation == False: self.calculate = self._calculate_classical_nebular - elif coronal_approximation == True: + elif coronal_approximation == True and classical_nebular == False: self.calculate = self._calculate_coronal_approximation + elif coronal_approximation == True and classical_nebular == True: + raise PlasmaConfigContradiction else: self.calculate = self._calculate_general self._update_inputs() From e687f9f92da846391835df32eefcbb1c11ac78e7 Mon Sep 17 00:00:00 2001 From: Aoife Boyle Date: Sat, 1 Aug 2015 02:38:49 +0100 Subject: [PATCH 13/13] Small fix. --- tardis/plasma/properties/partition_function.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tardis/plasma/properties/partition_function.py b/tardis/plasma/properties/partition_function.py index 2bcda0cf947..c97d7823d36 100644 --- a/tardis/plasma/properties/partition_function.py +++ b/tardis/plasma/properties/partition_function.py @@ -82,6 +82,7 @@ def _main_nlte_calculation(self, nlte_species, atomic_data, nlte_data, t_electron, j_blues, beta_sobolevs, general_level_boltzmann_factor, previous_electron_densities): for species in nlte_species: + j_blues = j_blues.values logger.info('Calculating rates for species %s', species) number_of_levels = atomic_data.levels.energy.ix[species].count() lnl = nlte_data.lines_level_number_lower[species] @@ -97,13 +98,13 @@ def _main_nlte_calculation(self, nlte_species, atomic_data, nlte_data, r_ul_matrix_reshaped = r_ul_matrix.reshape((number_of_levels**2, len(t_electron))) r_ul_matrix_reshaped[r_ul_index] = A_uls[np.newaxis].T + \ - B_uls[np.newaxis].T * j_blues.ix[lines_index] + B_uls[np.newaxis].T * j_blues[lines_index] r_ul_matrix_reshaped[r_ul_index] *= beta_sobolevs[lines_index] r_lu_matrix = np.zeros_like(r_ul_matrix) r_lu_matrix_reshaped = r_lu_matrix.reshape((number_of_levels**2, len(t_electron))) r_lu_matrix_reshaped[r_lu_index] = B_lus[np.newaxis].T * \ - j_blues.ix[lines_index] * beta_sobolevs[lines_index] + j_blues[lines_index] * beta_sobolevs[lines_index] if atomic_data.has_collision_data: if previous_electron_densities is None: collision_matrix = r_ul_matrix.copy()