From 55a190dbfa304d80ce1c63dc6408ad5e6050503b Mon Sep 17 00:00:00 2001 From: Swayam Shah Date: Thu, 2 Jan 2025 23:26:22 +0530 Subject: [PATCH 1/3] for test --- docs/tardis_example.yml | 13 +++++- docs/workflows/standard_workflow.ipynb | 59 +++++++++++++++++++++--- tardis/io/configuration/schemas/base.yml | 12 ++--- tardis/io/logger/logger.py | 49 ++++++++------------ 4 files changed, 89 insertions(+), 44 deletions(-) diff --git a/docs/tardis_example.yml b/docs/tardis_example.yml index 08d25b44eca..9224ddaa3f7 100644 --- a/docs/tardis_example.yml +++ b/docs/tardis_example.yml @@ -16,6 +16,9 @@ model: num: 20 density: type: branch85_w7 + w7_v_0: 2.5e9 cm/s + w7_rho_0: 1e-5 g/cm^3 + w7_time_0: 0.1 day abundances: type: uniform @@ -41,7 +44,7 @@ montecarlo: last_no_of_packets: 1.e+5 no_of_virtual_packets: 10 - + enable_full_relativity: false convergence_strategy: type: damped damping_constant: 1.0 @@ -51,7 +54,15 @@ montecarlo: t_inner: damping_constant: 0.5 + tracking: + initial_array_length: 1000 + track_rpacket: true + virtual_spectrum_spawn_range: 0.01 + debug_packets: false + logger_buffer: 1000 spectrum: start: 500 angstrom stop: 20000 angstrom num: 10000 + + diff --git a/docs/workflows/standard_workflow.ipynb b/docs/workflows/standard_workflow.ipynb index efdc3539c06..d698417431b 100644 --- a/docs/workflows/standard_workflow.ipynb +++ b/docs/workflows/standard_workflow.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -21,11 +21,58 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "workflow = StandardTARDISWorkflow(config, show_convergence_plots=True,show_progress_bars=True,convergence_plots_kwargs={\"export_convergence_plots\":True})" + "# from tardis.io.atom_data import download_atom_data\n", + "# download_atom_data('kurucz_cd23_chianti_H_He')" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[\u001b[1mtardis.io.model.parse_atom_data\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\t\n", + "\tReading Atomic Data from ..\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mparse_atom_data.py\u001b[0m:40)\n", + "[\u001b[1mtardis.io.atom_data.util\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\t\n", + "\tAtom Data kurucz_cd23_chianti_H_He.h5 not found in local path.\n", + "\tExists in TARDIS Data repo C:\\Users\\admin\\Downloads\\tardis-data\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mutil.py\u001b[0m:34)\n", + "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tReading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5 = 864f1753714343c41f99cb065710cace (\u001b[1mbase.py\u001b[0m:258)\n", + "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tNon provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules (\u001b[1mbase.py\u001b[0m:262)\n", + "[\u001b[1mtardis.io.model.parse_density_configuration\u001b[0m][\u001b[1;33mWARNING\u001b[0m] \n", + "\tNumber of density points larger than number of shells. Assuming inner point irrelevant (\u001b[1mparse_density_configuration.py\u001b[0m:114)\n", + "[\u001b[1mtardis.model.matter.decay\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tDecaying abundances for 1123200.0 seconds (\u001b[1mdecay.py\u001b[0m:101)\n", + "[\u001b[1mtardis.io.model.parse_atom_data\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\t\n", + "\tReading Atomic Data from ..\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mparse_atom_data.py\u001b[0m:40)\n", + "[\u001b[1mtardis.io.atom_data.util\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\t\n", + "\tAtom Data kurucz_cd23_chianti_H_He.h5 not found in local path.\n", + "\tExists in TARDIS Data repo C:\\Users\\admin\\Downloads\\tardis-data\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mutil.py\u001b[0m:34)\n", + "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tReading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5 = 864f1753714343c41f99cb065710cace (\u001b[1mbase.py\u001b[0m:258)\n", + "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tNon provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules (\u001b[1mbase.py\u001b[0m:262)\n", + "[\u001b[1mtardis.io.model.parse_density_configuration\u001b[0m][\u001b[1;33mWARNING\u001b[0m] \n", + "\tNumber of density points larger than number of shells. Assuming inner point irrelevant (\u001b[1mparse_density_configuration.py\u001b[0m:114)\n", + "[\u001b[1mtardis.model.matter.decay\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", + "\tDecaying abundances for 1123200.0 seconds (\u001b[1mdecay.py\u001b[0m:101)\n" + ] + } + ], + "source": [ + "workflow = StandardTARDISWorkflow(config, show_convergence_plots=False,show_progress_bars=False,convergence_plots_kwargs={\"export_convergence_plots\":False})" ] }, { @@ -88,7 +135,7 @@ ], "metadata": { "kernelspec": { - "display_name": "tardis", + "display_name": ".venv", "language": "python", "name": "python3" }, @@ -102,7 +149,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.4" + "version": "3.10.11" } }, "nbformat": 4, diff --git a/tardis/io/configuration/schemas/base.yml b/tardis/io/configuration/schemas/base.yml index 492379aa322..9250a1d8a0f 100644 --- a/tardis/io/configuration/schemas/base.yml +++ b/tardis/io/configuration/schemas/base.yml @@ -6,28 +6,28 @@ properties: type: string description: Version of the configuration file. The current version is 1.0 and no other versions are allowed supernova: - $ref: supernova.yml + $ref: "..\\supernova.yml" description: a section pertaining to observations of the supernova atom_data: type: string description: path or filename to the Atomic Data HDF5 file plasma: - $ref: plasma.yml + $ref: "..\\plasma.yml" description: configuration of the plasma microphysics csvy_model: type: string description: defining the model using csvy format (see :ref:`csvy-models`) model: - $ref: model.yml + $ref: "..\\model_definitions.yml" description: defining the model in the config yaml file montecarlo: - $ref: montecarlo.yml + $ref: "..\\montecarlo_definitions.yml" description: configuring the physics of the monte carlo radiative transfer spectrum: - $ref: spectrum.yml + $ref: "..\\spectrum.yml" description: Final spectrum sampling debug: - $ref: debug.yml + $ref: "..\\debug.yml" description: Debugging setup for the simulation required: - tardis_config_version diff --git a/tardis/io/logger/logger.py b/tardis/io/logger/logger.py index 44da5fa187f..f9dda8a6b33 100644 --- a/tardis/io/logger/logger.py +++ b/tardis/io/logger/logger.py @@ -77,41 +77,28 @@ def logging_state(log_level, tardis_config, specific_log_level): specific_log_level: boolean Allows to set specific logging levels. Logs of the `log_level` level would be output. """ - if "debug" in tardis_config: - specific_log_level = ( - tardis_config["debug"]["specific_log_level"] - if specific_log_level is None - else specific_log_level - ) - - logging_level = ( - log_level if log_level else tardis_config["debug"]["log_level"] - ) + # Ensure the debug section exists + if "debug" not in tardis_config: + tardis_config["debug"] = {} - # Displays a message when both log_level & tardis["debug"]["log_level"] are specified - if log_level and tardis_config["debug"]["log_level"]: - print( - "log_level is defined both in Functional Argument & YAML Configuration {debug section}" - ) - print( - f"log_level = {log_level.upper()} will be used for Log Level Determination\n" - ) + # Set specific_log_level if not provided + if specific_log_level is None: + specific_log_level = tardis_config["debug"].get("specific_log_level", DEFAULT_SPECIFIC_STATE) + # Set logging_level if not provided + if log_level is None: + logging_level = tardis_config["debug"].get("log_level", DEFAULT_LOG_LEVEL) else: - # Adds empty `debug` section for the YAML - tardis_config["debug"] = {} + logging_level = log_level - if log_level: - logging_level = log_level - else: - tardis_config["debug"]["log_level"] = DEFAULT_LOG_LEVEL - logging_level = tardis_config["debug"]["log_level"] - - if not specific_log_level: - tardis_config["debug"][ - "specific_log_level" - ] = DEFAULT_SPECIFIC_STATE - specific_log_level = tardis_config["debug"]["specific_log_level"] + # Displays a message when both log_level & tardis["debug"]["log_level"] are specified + if log_level and tardis_config["debug"].get("log_level"): + print( + "log_level is defined both in Functional Argument & YAML Configuration {debug section}" + ) + print( + f"log_level = {log_level.upper()} will be used for Log Level Determination\n" + ) logging_level = logging_level.upper() if logging_level not in LOGGING_LEVELS: From 6d3b92c0b75a1e27d5dead373cec65ae5ae902d3 Mon Sep 17 00:00:00 2001 From: Swayam Shah Date: Sun, 12 Jan 2025 01:52:02 +0530 Subject: [PATCH 2/3] Revert "for test" This reverts commit 55a190dbfa304d80ce1c63dc6408ad5e6050503b. --- docs/tardis_example.yml | 13 +----- docs/workflows/standard_workflow.ipynb | 59 +++--------------------- tardis/io/configuration/schemas/base.yml | 12 ++--- tardis/io/logger/logger.py | 49 ++++++++++++-------- 4 files changed, 44 insertions(+), 89 deletions(-) diff --git a/docs/tardis_example.yml b/docs/tardis_example.yml index 9224ddaa3f7..08d25b44eca 100644 --- a/docs/tardis_example.yml +++ b/docs/tardis_example.yml @@ -16,9 +16,6 @@ model: num: 20 density: type: branch85_w7 - w7_v_0: 2.5e9 cm/s - w7_rho_0: 1e-5 g/cm^3 - w7_time_0: 0.1 day abundances: type: uniform @@ -44,7 +41,7 @@ montecarlo: last_no_of_packets: 1.e+5 no_of_virtual_packets: 10 - enable_full_relativity: false + convergence_strategy: type: damped damping_constant: 1.0 @@ -54,15 +51,7 @@ montecarlo: t_inner: damping_constant: 0.5 - tracking: - initial_array_length: 1000 - track_rpacket: true - virtual_spectrum_spawn_range: 0.01 - debug_packets: false - logger_buffer: 1000 spectrum: start: 500 angstrom stop: 20000 angstrom num: 10000 - - diff --git a/docs/workflows/standard_workflow.ipynb b/docs/workflows/standard_workflow.ipynb index d698417431b..efdc3539c06 100644 --- a/docs/workflows/standard_workflow.ipynb +++ b/docs/workflows/standard_workflow.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -21,58 +21,11 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# from tardis.io.atom_data import download_atom_data\n", - "# download_atom_data('kurucz_cd23_chianti_H_He')" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[\u001b[1mtardis.io.model.parse_atom_data\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\t\n", - "\tReading Atomic Data from ..\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mparse_atom_data.py\u001b[0m:40)\n", - "[\u001b[1mtardis.io.atom_data.util\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\t\n", - "\tAtom Data kurucz_cd23_chianti_H_He.h5 not found in local path.\n", - "\tExists in TARDIS Data repo C:\\Users\\admin\\Downloads\\tardis-data\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mutil.py\u001b[0m:34)\n", - "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tReading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5 = 864f1753714343c41f99cb065710cace (\u001b[1mbase.py\u001b[0m:258)\n", - "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tNon provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules (\u001b[1mbase.py\u001b[0m:262)\n", - "[\u001b[1mtardis.io.model.parse_density_configuration\u001b[0m][\u001b[1;33mWARNING\u001b[0m] \n", - "\tNumber of density points larger than number of shells. Assuming inner point irrelevant (\u001b[1mparse_density_configuration.py\u001b[0m:114)\n", - "[\u001b[1mtardis.model.matter.decay\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tDecaying abundances for 1123200.0 seconds (\u001b[1mdecay.py\u001b[0m:101)\n", - "[\u001b[1mtardis.io.model.parse_atom_data\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\t\n", - "\tReading Atomic Data from ..\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mparse_atom_data.py\u001b[0m:40)\n", - "[\u001b[1mtardis.io.atom_data.util\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\t\n", - "\tAtom Data kurucz_cd23_chianti_H_He.h5 not found in local path.\n", - "\tExists in TARDIS Data repo C:\\Users\\admin\\Downloads\\tardis-data\\kurucz_cd23_chianti_H_He.h5 (\u001b[1mutil.py\u001b[0m:34)\n", - "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tReading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5 = 864f1753714343c41f99cb065710cace (\u001b[1mbase.py\u001b[0m:258)\n", - "[\u001b[1mtardis.io.atom_data.base\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tNon provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules (\u001b[1mbase.py\u001b[0m:262)\n", - "[\u001b[1mtardis.io.model.parse_density_configuration\u001b[0m][\u001b[1;33mWARNING\u001b[0m] \n", - "\tNumber of density points larger than number of shells. Assuming inner point irrelevant (\u001b[1mparse_density_configuration.py\u001b[0m:114)\n", - "[\u001b[1mtardis.model.matter.decay\u001b[0m][\u001b[1;37mINFO\u001b[0m ] \n", - "\tDecaying abundances for 1123200.0 seconds (\u001b[1mdecay.py\u001b[0m:101)\n" - ] - } - ], - "source": [ - "workflow = StandardTARDISWorkflow(config, show_convergence_plots=False,show_progress_bars=False,convergence_plots_kwargs={\"export_convergence_plots\":False})" + "workflow = StandardTARDISWorkflow(config, show_convergence_plots=True,show_progress_bars=True,convergence_plots_kwargs={\"export_convergence_plots\":True})" ] }, { @@ -135,7 +88,7 @@ ], "metadata": { "kernelspec": { - "display_name": ".venv", + "display_name": "tardis", "language": "python", "name": "python3" }, @@ -149,7 +102,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.11" + "version": "3.12.4" } }, "nbformat": 4, diff --git a/tardis/io/configuration/schemas/base.yml b/tardis/io/configuration/schemas/base.yml index 9250a1d8a0f..492379aa322 100644 --- a/tardis/io/configuration/schemas/base.yml +++ b/tardis/io/configuration/schemas/base.yml @@ -6,28 +6,28 @@ properties: type: string description: Version of the configuration file. The current version is 1.0 and no other versions are allowed supernova: - $ref: "..\\supernova.yml" + $ref: supernova.yml description: a section pertaining to observations of the supernova atom_data: type: string description: path or filename to the Atomic Data HDF5 file plasma: - $ref: "..\\plasma.yml" + $ref: plasma.yml description: configuration of the plasma microphysics csvy_model: type: string description: defining the model using csvy format (see :ref:`csvy-models`) model: - $ref: "..\\model_definitions.yml" + $ref: model.yml description: defining the model in the config yaml file montecarlo: - $ref: "..\\montecarlo_definitions.yml" + $ref: montecarlo.yml description: configuring the physics of the monte carlo radiative transfer spectrum: - $ref: "..\\spectrum.yml" + $ref: spectrum.yml description: Final spectrum sampling debug: - $ref: "..\\debug.yml" + $ref: debug.yml description: Debugging setup for the simulation required: - tardis_config_version diff --git a/tardis/io/logger/logger.py b/tardis/io/logger/logger.py index f9dda8a6b33..44da5fa187f 100644 --- a/tardis/io/logger/logger.py +++ b/tardis/io/logger/logger.py @@ -77,28 +77,41 @@ def logging_state(log_level, tardis_config, specific_log_level): specific_log_level: boolean Allows to set specific logging levels. Logs of the `log_level` level would be output. """ - # Ensure the debug section exists - if "debug" not in tardis_config: - tardis_config["debug"] = {} + if "debug" in tardis_config: + specific_log_level = ( + tardis_config["debug"]["specific_log_level"] + if specific_log_level is None + else specific_log_level + ) + + logging_level = ( + log_level if log_level else tardis_config["debug"]["log_level"] + ) - # Set specific_log_level if not provided - if specific_log_level is None: - specific_log_level = tardis_config["debug"].get("specific_log_level", DEFAULT_SPECIFIC_STATE) + # Displays a message when both log_level & tardis["debug"]["log_level"] are specified + if log_level and tardis_config["debug"]["log_level"]: + print( + "log_level is defined both in Functional Argument & YAML Configuration {debug section}" + ) + print( + f"log_level = {log_level.upper()} will be used for Log Level Determination\n" + ) - # Set logging_level if not provided - if log_level is None: - logging_level = tardis_config["debug"].get("log_level", DEFAULT_LOG_LEVEL) else: - logging_level = log_level + # Adds empty `debug` section for the YAML + tardis_config["debug"] = {} - # Displays a message when both log_level & tardis["debug"]["log_level"] are specified - if log_level and tardis_config["debug"].get("log_level"): - print( - "log_level is defined both in Functional Argument & YAML Configuration {debug section}" - ) - print( - f"log_level = {log_level.upper()} will be used for Log Level Determination\n" - ) + if log_level: + logging_level = log_level + else: + tardis_config["debug"]["log_level"] = DEFAULT_LOG_LEVEL + logging_level = tardis_config["debug"]["log_level"] + + if not specific_log_level: + tardis_config["debug"][ + "specific_log_level" + ] = DEFAULT_SPECIFIC_STATE + specific_log_level = tardis_config["debug"]["specific_log_level"] logging_level = logging_level.upper() if logging_level not in LOGGING_LEVELS: From 8ff2fa7e2a87dd0862c616ff1b10c39ab383e523 Mon Sep 17 00:00:00 2001 From: Swayam Shah Date: Sun, 26 Jan 2025 19:13:56 +0530 Subject: [PATCH 3/3] moved get_tau_integ --- tardis/opacities/opacity_state.py | 92 ++++++++++++++++++++++++++++++ tardis/workflows/util.py | 92 ------------------------------ tardis/workflows/v_inner_solver.py | 5 +- 3 files changed, 94 insertions(+), 95 deletions(-) diff --git a/tardis/opacities/opacity_state.py b/tardis/opacities/opacity_state.py index 6db3b5791b9..5fee61deb3d 100644 --- a/tardis/opacities/opacity_state.py +++ b/tardis/opacities/opacity_state.py @@ -6,6 +6,7 @@ from tardis.opacities.macro_atom.macroatom_state import MacroAtomState from tardis.opacities.tau_sobolev import calculate_sobolev_line_opacity from tardis.transport.montecarlo.configuration import montecarlo_globals +from tardis.workflows.util import * opacity_state_spec = [ ("electron_density", float64[:]), @@ -364,6 +365,97 @@ def to_numba( k_packet_idx, ) + def get_tau_integ(self, plasma, simulation_state, bin_size=10): + """Estimate the integrated mean optical depth at each velocity bin + + Parameters + ---------- + plasma : tardis.plasma.BasePlasma + The tardis legacy plasma + simulation_state : tardis.model.base.SimulationState + the current simulation state + bin_size : int, optional. Default : 10 + bin size for the aggregation of line opacities + + Returns + ------- + dict + rosseland : np.ndarray + Roassland Mean Optical Depth + planck : np.ndarray + Planck Mean Optical Depth + """ + index = plasma.atomic_data.lines.nu.index + freqs = plasma.atomic_data.lines.nu.values * u.Hz + order = np.argsort(freqs) + freqs = freqs[order] + taus = self.tau_sobolev.values[order] + + check_bin_size = True + while check_bin_size: + extra = bin_size - len(freqs) % bin_size + extra_freqs = (np.arange(extra + 1) + 1) * u.Hz + extra_taus = np.zeros((extra + 1, taus.shape[1])) + freqs = np.hstack((extra_freqs, freqs)) + taus = np.vstack((extra_taus, taus)) + + bins_low = freqs[:-bin_size:bin_size] + bins_high = freqs[bin_size::bin_size] + delta_nu = bins_high - bins_low + n_bins = len(delta_nu) + + if np.any(delta_nu == 0): + bin_size += 2 + else: + check_bin_size = False + + taus = taus[1 : n_bins * bin_size + 1] + freqs = freqs[1 : n_bins * bin_size + 1] + + ct = simulation_state.time_explosion * const.c + t_rad = simulation_state.radiation_field_state.temperature + + def B(nu, T): + return ( + 2 + * const.h + * nu**3 + / const.c**2 + / (np.exp(const.h * nu / (const.k_B * T)) - 1) + ) + + def U(nu, T): + return ( + B(nu, T) ** 2 * (const.c / nu) ** 2 * (2 * const.k_B * T**2) ** -1 + ) + + kappa_exp = ( + (bins_low / delta_nu).reshape(-1, 1) + / ct + * (1 - np.exp(-taus.reshape(n_bins, bin_size, -1))).sum(axis=1) + ) + kappa_thom = plasma.electron_densities.values * u.cm ** (-3) * const.sigma_T + Bdnu = B(bins_low.reshape(-1, 1), t_rad.reshape(1, -1)) * delta_nu.reshape( + -1, 1 + ) + kappa_planck = kappa_thom + (Bdnu * kappa_exp).sum(axis=0) / ( + Bdnu.sum(axis=0) + ) + + udnu = U(bins_low.reshape(-1, 1), t_rad.reshape(1, -1)) * delta_nu.reshape( + -1, 1 + ) + kappa_tot = kappa_thom + kappa_exp + kappa_rosseland = ( + (udnu * kappa_tot**-1).sum(axis=0) / (udnu.sum(axis=0)) + ) ** -1 + + dr = simulation_state.geometry.r_outer - simulation_state.geometry.r_inner + dtau = kappa_planck * dr + planck_integ_tau = np.cumsum(dtau[::-1])[::-1] + rosseland_integ_tau = np.cumsum((kappa_rosseland * dr)[::-1])[::-1] + + return {"rosseland": rosseland_integ_tau, "planck": planck_integ_tau} def opacity_state_initialize( plasma, diff --git a/tardis/workflows/util.py b/tardis/workflows/util.py index 43aee66989d..bb6f4ed0411 100644 --- a/tardis/workflows/util.py +++ b/tardis/workflows/util.py @@ -2,95 +2,3 @@ from astropy import constants as const from astropy import units as u - -def get_tau_integ(plasma, opacity_state, simulation_state, bin_size=10): - """Estimate the integrated mean optical depth at each velocity bin - - Parameters - ---------- - plasma : tardis.plasma.BasePlasma - The tardis legacy plasma - simulation_state : tardis.model.base.SimulationState - the current simulation state - bin_size : int, optional. Default : 10 - bin size for the aggregation of line opacities - - Returns - ------- - dict - rosseland : np.ndarray - Roassland Mean Optical Depth - planck : np.ndarray - Planck Mean Optical Depth - """ - index = plasma.atomic_data.lines.nu.index - freqs = plasma.atomic_data.lines.nu.values * u.Hz - order = np.argsort(freqs) - freqs = freqs[order] - taus = opacity_state.tau_sobolev.values[order] - - check_bin_size = True - while check_bin_size: - extra = bin_size - len(freqs) % bin_size - extra_freqs = (np.arange(extra + 1) + 1) * u.Hz - extra_taus = np.zeros((extra + 1, taus.shape[1])) - freqs = np.hstack((extra_freqs, freqs)) - taus = np.vstack((extra_taus, taus)) - - bins_low = freqs[:-bin_size:bin_size] - bins_high = freqs[bin_size::bin_size] - delta_nu = bins_high - bins_low - n_bins = len(delta_nu) - - if np.any(delta_nu == 0): - bin_size += 2 - else: - check_bin_size = False - - taus = taus[1 : n_bins * bin_size + 1] - freqs = freqs[1 : n_bins * bin_size + 1] - - ct = simulation_state.time_explosion * const.c - t_rad = simulation_state.radiation_field_state.temperature - - def B(nu, T): - return ( - 2 - * const.h - * nu**3 - / const.c**2 - / (np.exp(const.h * nu / (const.k_B * T)) - 1) - ) - - def U(nu, T): - return ( - B(nu, T) ** 2 * (const.c / nu) ** 2 * (2 * const.k_B * T**2) ** -1 - ) - - kappa_exp = ( - (bins_low / delta_nu).reshape(-1, 1) - / ct - * (1 - np.exp(-taus.reshape(n_bins, bin_size, -1))).sum(axis=1) - ) - kappa_thom = plasma.electron_densities.values * u.cm ** (-3) * const.sigma_T - Bdnu = B(bins_low.reshape(-1, 1), t_rad.reshape(1, -1)) * delta_nu.reshape( - -1, 1 - ) - kappa_planck = kappa_thom + (Bdnu * kappa_exp).sum(axis=0) / ( - Bdnu.sum(axis=0) - ) - - udnu = U(bins_low.reshape(-1, 1), t_rad.reshape(1, -1)) * delta_nu.reshape( - -1, 1 - ) - kappa_tot = kappa_thom + kappa_exp - kappa_rosseland = ( - (udnu * kappa_tot**-1).sum(axis=0) / (udnu.sum(axis=0)) - ) ** -1 - - dr = simulation_state.geometry.r_outer - simulation_state.geometry.r_inner - dtau = kappa_planck * dr - planck_integ_tau = np.cumsum(dtau[::-1])[::-1] - rosseland_integ_tau = np.cumsum((kappa_rosseland * dr)[::-1])[::-1] - - return {"rosseland": rosseland_integ_tau, "planck": planck_integ_tau} diff --git a/tardis/workflows/v_inner_solver.py b/tardis/workflows/v_inner_solver.py index 20edd13c8c9..5cc88231edf 100644 --- a/tardis/workflows/v_inner_solver.py +++ b/tardis/workflows/v_inner_solver.py @@ -5,10 +5,10 @@ from astropy import units as u from scipy.interpolate import interp1d +from tardis.opacities.opacity_state import OpacityState from tardis.plasma.radiation_field import DilutePlanckianRadiationField from tardis.simulation.convergence import ConvergenceSolver from tardis.workflows.simple_tardis_workflow import SimpleTARDISWorkflow -from tardis.workflows.util import get_tau_integ # logging support logger = logging.getLogger(__name__) @@ -54,9 +54,8 @@ def estimate_v_inner(self): Need some way to return and inspect the optical depths for later logging """ tau_integ = np.log( - get_tau_integ( + self.opacity_states["opacity_state"].get_tau_integ( self.plasma_solver, - self.opacity_states["opacity_state"], self.simulation_state, )[self.mean_optical_depth] )