Skip to content

Commit

Permalink
V inner formal integral (#2800)
Browse files Browse the repository at this point in the history
* Basic skeleton

* Lots of progress

Solvers pass info to each other
Simplified some parts

* Adds setup method, refactors some functions, removes unused items

* Move convergence setup to setup method

* Add docs notebook, fix a couple of errors

* Move setup to init

* Fix iteration and convergence

* Simplify spectrumsolver init

* Separate plasma "solver" and simulation state updates

* Simplify convergence solver setup with a dict.

* Minor formatting change

* Minor refactoring

* Fixes plasma update step

* Fixes loggers and progress bars

Also adds docstrings to methods. Updates some methods to use new functionality of the plasma. Adds requirements for the convergence plots (still broken)

* Fixes convergence plot rendering

* Fixes convergence plots in the final iteration

* black

* Simplify convergence plot updating

* Move logging handling to a separate class

* Add HDF output capability to solver

* Move more basic logging back into workflow

* Add not-converged error message

* Update notebook with export option

* Added simple base workflow and changed some verbiage

* Fix typo

* Moved solver workflow over to updated branch

* Added reprojection method for comparing arrays between iterations

* Updated schema and parser to add v_inner_boundary, added property to the model because it is named weirdly

* Added ntoebook for the IBVS workflow

* Updated workflow notebook to have correct convergence information, updated simple solver to properly set up the plasma with the new assembly

* Reduced size of overloaded functions to make them more in-line with the simple solver, updated simple solver to handle current tardis master branch

* cleared notebook outputs

* Updated notebook to use correctmethod

* Updated the docstrings

* removed unused imports

* Updated notebook to stop if converged

* Updated v_inner_solver to deal with detailed radiative rates types

* removed old comment

* Added docstrings, removed old comments

* Fixed typo

* Updated notebook to be more explicit

* removed old comment

* removed unused variable

* Changed the radiation_field to use the radiation field state properties with no explicit sovler method

* Added back radaition field to plasma solver in v_inner solver

* Added a solve opacity method

* Added docstring to solve_opacity method

* Fixed initializing opacity at the last iteration

* Added initial v_inner estimate in the __init__

* massively simplified the reproject method by just doing basic math

* added docstring describing the simplified reproject method so I don't forget why it works later

* Formatting error in doscstring

* Added nice logger output showing changes in geometry state

* Cleaned up the notebook a bit

* Updated the formal integral to handle changes in geometry

* Fixed the standard simulation solver

* Remove old files

* Fixes v_inner workflow

* Minor formatting changes

* Clean up confusing v_boundary_inner vs v_inner_boundary variable names

---------

Co-authored-by: Andrew Fullard <[email protected]>
  • Loading branch information
Rodot- and andrewfullard authored Jan 14, 2025
1 parent 0a0ecbf commit 5e2d0bb
Show file tree
Hide file tree
Showing 10 changed files with 655 additions and 34 deletions.
8 changes: 4 additions & 4 deletions docs/physics/setup/model.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@
"id": "cee054e9",
"metadata": {},
"source": [
"In the cell below, we set up a model. We use the [specific structure](../../io/configuration/components/models/index.rst#specific-structure) where we supply $t_\\mathrm{explosion}$, the velocity of the inner and outer boundaries of the supernova (labeled `start` and `stop`), and the number of shells (labeled `num`). The shells are then evenly spaced between the inner and outer boundaries of the supernova. The time after the explosion, the inner and outer velocities, and the number of shells can be varied to get different shell structures. The `SimulationState` object stores information about the model in the following attributes: `velocity` shows the velocity of the shell boundaries, `v_inner` shows the velocities of the inner boundaries of each shell, `v_outer` shows the velocity of the outer boundaries of each shell, and `v_middle` shows the velocity of the middle of each shell. Similarly, `radius`, `r_inner`, `r_outer`, and `r_middle` show the radii of each shell boundary, the inner boundaries, the outer boundaries, and the middles of each shell, respectively. `v_boundary_inner` shows the velocity of the inner boundary of the computational domain, and `v_boundary_outer` shows the velocity of the outer boundary of the computational domain. Finally, `volume` shows the volume of each shell, calculated via the formula of the volume of a spherical shell: $V=\\frac{4}{3}\\pi (r_\\mathrm{outer}^3-r_\\mathrm{inner}^3)$."
"In the cell below, we set up a model. We use the [specific structure](../../io/configuration/components/models/index.rst#specific-structure) where we supply $t_\\mathrm{explosion}$, the velocity of the inner and outer boundaries of the supernova (labeled `start` and `stop`), and the number of shells (labeled `num`). The shells are then evenly spaced between the inner and outer boundaries of the supernova. The time after the explosion, the inner and outer velocities, and the number of shells can be varied to get different shell structures. The `SimulationState` object stores information about the model in the following attributes: `velocity` shows the velocity of the shell boundaries, `v_inner` shows the velocities of the inner boundaries of each shell, `v_outer` shows the velocity of the outer boundaries of each shell, and `v_middle` shows the velocity of the middle of each shell. Similarly, `radius`, `r_inner`, `r_outer`, and `r_middle` show the radii of each shell boundary, the inner boundaries, the outer boundaries, and the middles of each shell, respectively. `v_inner_boundary` shows the velocity of the inner boundary of the computational domain, and `v_outer_boundary` shows the velocity of the outer boundary of the computational domain. Finally, `volume` shows the volume of each shell, calculated via the formula of the volume of a spherical shell: $V=\\frac{4}{3}\\pi (r_\\mathrm{outer}^3-r_\\mathrm{inner}^3)$."
]
},
{
Expand Down Expand Up @@ -111,8 +111,8 @@
"print('v_inner:\\n', shell_simulation_state.v_inner)\n",
"print('v_outer:\\n', shell_simulation_state.v_outer)\n",
"print('v_middle:\\n', shell_simulation_state.v_middle)\n",
"print('v_boundary_inner:\\n', shell_simulation_state.v_boundary_inner)\n",
"print('v_boundary_outer:\\n', shell_simulation_state.v_boundary_outer)\n",
"print('v_inner_boundary:\\n', shell_simulation_state.v_inner_boundary)\n",
"print('v_outer_boundary:\\n', shell_simulation_state.v_outer_boundary)\n",
"print('radius:\\n', shell_simulation_state.radius)\n",
"print('r_inner:\\n', shell_simulation_state.r_inner)\n",
"print('r_outer:\\n', shell_simulation_state.r_outer)\n",
Expand All @@ -125,7 +125,7 @@
"id": "1ee56110",
"metadata": {},
"source": [
"Notice that `radius = velocity*time_explosion`, and similarly for `r_inner`, `r_outer`, and `r_middle`. You can get the radius of the photosphere via `v_boundary_inner*time_explosion` and outer edge of the supernova via `v_boundary_outer*time_explosion`.\n",
"Notice that `radius = velocity*time_explosion`, and similarly for `r_inner`, `r_outer`, and `r_middle`. You can get the radius of the photosphere via `v_inner_boundary*time_explosion` and outer edge of the supernova via `v_boundary_outer*time_explosion`.\n",
"\n",
"<div class=\"alert alert-info\">\n",
" \n",
Expand Down
128 changes: 128 additions & 0 deletions docs/workflows/v_inner_solver_workflow.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from tardis.workflows.v_inner_solver import InnerVelocitySolverWorkflow\n",
"from tardis.io.configuration.config_reader import Configuration"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"config = Configuration.from_yaml('../tardis_example.yml')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This code modifies the TARDIS example configuration to include convergence information for the inner boundary velocity solver.\n",
"Note that the number of shells is increased and the starting velocity is reduced to provide more granularity and a wider search window, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"from astropy import units as u\n",
"\n",
"config.montecarlo.convergence_strategy['v_inner_boundary'] = {\n",
" 'damping_constant' : 0.5,\n",
" 'threshold' : 0.01,\n",
" 'type' : 'damped'\n",
" }\n",
"\n",
"config.montecarlo.convergence_strategy.stop_if_converged = True\n",
"config.model.structure.velocity.start = 5000 * u.km/u.s # Larger window over which to search\n",
"config.model.structure.velocity.num = 50 # Increase number of shells\n",
"\n",
"workflow = InnerVelocitySolverWorkflow(\n",
" config, tau=2.0/3,\n",
" mean_optical_depth=\"rosseland\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"workflow.run()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"spectrum = workflow.spectrum_solver.spectrum_real_packets\n",
"spectrum_virtual = workflow.spectrum_solver.spectrum_virtual_packets\n",
"spectrum_integrated = workflow.spectrum_solver.spectrum_integrated"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"plt.figure(figsize=(10, 6.5))\n",
"\n",
"spectrum.plot(label=\"Normal packets\")\n",
"spectrum_virtual.plot(label=\"Virtual packets\")\n",
"spectrum_integrated.plot(label='Formal integral')\n",
"\n",
"plt.xlim(500, 9000)\n",
"plt.title(\"TARDIS example model spectrum\")\n",
"plt.xlabel(r\"Wavelength [$\\AA$]\")\n",
"plt.ylabel(r\"Luminosity density [erg/s/$\\AA$]\")\n",
"plt.legend()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "tardis",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
2 changes: 1 addition & 1 deletion tardis/io/configuration/config_reader.py
Original file line number Diff line number Diff line change
Expand Up @@ -431,7 +431,7 @@ def parse_convergence_section(convergence_section_dict):
"""
convergence_parameters = ["damping_constant", "threshold", "type"]

for convergence_variable in ["t_inner", "t_rad", "w"]:
for convergence_variable in ["t_inner", "t_rad", "w", "v_inner_boundary"]:
if convergence_variable not in convergence_section_dict:
convergence_section_dict[convergence_variable] = {}
convergence_variable_section = convergence_section_dict[
Expand Down
19 changes: 18 additions & 1 deletion tardis/io/configuration/schemas/montecarlo_definitions.yml
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ definitions:
title: 'Damped Convergence Strategy'
type: object
additionalProperties: false
properties:
properties:
type:
enum:
Expand Down Expand Up @@ -59,6 +58,24 @@ definitions:
type: string
default: 'damped'
description: THIS IS A DUMMY VARIABLE DO NOT USE
v_inner_boundary:
type: object
additionalProperties: false
properties:
damping_constant:
type: number
default: 0.0
description: damping constant
minimum: 0
threshold:
type: number
description: specifies the threshold that is taken as convergence (i.e.
0.05 means that the value does not change more than 5%)
minimum: 0
type:
type: string
default: 'damped'
description: THIS IS A DUMMY VARIABLE DO NOT USE
t_rad:
type: object
additionalProperties: false
Expand Down
16 changes: 8 additions & 8 deletions tardis/io/model/parse_geometry_configuration.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,17 +135,17 @@ def parse_geometry_from_csvy(
"""
if hasattr(config, "model"):
if hasattr(config.model, "v_inner_boundary"):
v_boundary_inner = config.model.v_inner_boundary
v_inner_boundary = config.model.v_inner_boundary
else:
v_boundary_inner = None
v_inner_boundary = None

if hasattr(config.model, "v_outer_boundary"):
v_boundary_outer = config.model.v_outer_boundary
v_outer_boundary = config.model.v_outer_boundary
else:
v_boundary_outer = None
v_outer_boundary = None
else:
v_boundary_inner = None
v_boundary_outer = None
v_inner_boundary = None
v_outer_boundary = None

if hasattr(csvy_model_config, "velocity"):
velocity = quantity_linspace(
Expand All @@ -166,8 +166,8 @@ def parse_geometry_from_csvy(
geometry = HomologousRadial1DGeometry(
velocity[:-1], # v_inner
velocity[1:], # v_outer
v_inner_boundary=v_boundary_inner,
v_outer_boundary=v_boundary_outer,
v_inner_boundary=v_inner_boundary,
v_outer_boundary=v_outer_boundary,
time_explosion=time_explosion,
)
return geometry
14 changes: 7 additions & 7 deletions tardis/model/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,11 +58,11 @@ class SimulationState(HDFWriterMixin):
dilution_factor : np.ndarray
If None, the dilution_factor will be initialized with the geometric
dilution factor.
v_boundary_inner : astropy.units.Quantity
v_boundary_outer : astropy.units.Quantity
v_inner_boundary : astropy.units.Quantity
v_outer_boundary : astropy.units.Quantity
raw_velocity : np.ndarray
The complete array of the velocities, without being cut by
`v_boundary_inner` and `v_boundary_outer`
`v_inner_boundary` and `v_outer_boundary`
electron_densities : astropy.units.quantity.Quantity
Attributes
Expand All @@ -81,8 +81,8 @@ class SimulationState(HDFWriterMixin):
density : astropy.units.quantity.Quantity
volume : astropy.units.quantity.Quantity
no_of_shells : int
The number of shells as formed by `v_boundary_inner` and
`v_boundary_outer`
The number of shells as formed by `v_inner_boundary` and
`v_outer_boundary`
no_of_raw_shells : int
"""

Expand Down Expand Up @@ -206,11 +206,11 @@ def radius(self):
return self.time_explosion * self.velocity

@property
def v_boundary_inner(self):
def v_inner_boundary(self):
return self.geometry.v_inner_boundary

@property
def v_boundary_outer(self):
def v_outer_boundary(self):
return self.geometry.v_outer_boundary

@property
Expand Down
44 changes: 31 additions & 13 deletions tardis/spectrum/formal_integral.py
Original file line number Diff line number Diff line change
Expand Up @@ -434,9 +434,22 @@ def make_source_function(self):
-------
Numpy array containing ( 1 - exp(-tau_ul) ) S_ul ordered by wavelength of the transition u -> l
"""

simulation_state = self.simulation_state
# slice for the active shells
local_slice = slice(
simulation_state.geometry.v_inner_boundary_index,
simulation_state.geometry.v_outer_boundary_index,
)

transport = self.transport
montecarlo_transport_state = transport.transport_state
transition_probabilities = self.opacity_state.transition_probabilities[
:, local_slice
]
tau_sobolevs = self.opacity_state.tau_sobolev[:, local_slice]

columns = range(simulation_state.no_of_shells)

# macro_ref = self.atomic_data.macro_atom_references
macro_ref = self.atomic_data.macro_atom_references
Expand All @@ -449,9 +462,7 @@ def make_source_function(self):
if transport.line_interaction_type == "macroatom":
internal_jump_mask = (macro_data.transition_type >= 0).values
ma_int_data = macro_data[internal_jump_mask]
internal = self.opacity_state.transition_probabilities[
internal_jump_mask
]
internal = transition_probabilities[internal_jump_mask]

source_level_idx = ma_int_data.source_level_idx.values
destination_level_idx = ma_int_data.destination_level_idx.values
Expand All @@ -460,7 +471,7 @@ def make_source_function(self):
montecarlo_transport_state.packet_collection.time_of_simulation
* simulation_state.volume
)
exptau = 1 - np.exp(-self.opacity_state.tau_sobolev)
exptau = 1 - np.exp(-tau_sobolevs)
Edotlu = (
Edotlu_norm_factor
* exptau
Expand Down Expand Up @@ -493,14 +504,14 @@ def make_source_function(self):
upper_level_index = self.atomic_data.lines.index.droplevel(
"level_number_lower"
)
e_dot_lu = pd.DataFrame(Edotlu.value, index=upper_level_index)
e_dot_lu = pd.DataFrame(
Edotlu.value, index=upper_level_index, columns=columns
)
e_dot_u = e_dot_lu.groupby(level=[0, 1, 2]).sum()
e_dot_u_src_idx = macro_ref.loc[e_dot_u.index].references_idx.values

if transport.line_interaction_type == "macroatom":
C_frame = pd.DataFrame(
columns=np.arange(no_shells), index=macro_ref.index
)
C_frame = pd.DataFrame(columns=columns, index=macro_ref.index)
q_indices = (source_level_idx, destination_level_idx)
for shell in range(no_shells):
Q = sp.coo_matrix(
Expand All @@ -523,7 +534,7 @@ def make_source_function(self):
["atomic_number", "ion_number", "source_level_number"]
).index.copy()
tmp = pd.DataFrame(
self.opacity_state.transition_probabilities[
transition_probabilities[
(self.atomic_data.macro_atom_data.transition_type == -1).values
]
)
Expand All @@ -539,13 +550,15 @@ def make_source_function(self):
att_S_ul = wave * (q_ul * e_dot_u) * t / (4 * np.pi)

result = pd.DataFrame(
att_S_ul.values, index=transitions.transition_line_id.values
att_S_ul.values,
index=transitions.transition_line_id.values,
columns=columns,
)
att_S_ul = result.loc[lines.index.values].values

# Jredlu should already by in the correct order, i.e. by wavelength of
# the transition l->u (similar to Jbluelu)
Jredlu = Jbluelu * np.exp(-self.opacity_state.tau_sobolev) + att_S_ul
Jredlu = Jbluelu * np.exp(-tau_sobolevs) + att_S_ul
if self.interpolate_shells > 0:
(
att_S_ul,
Expand Down Expand Up @@ -592,15 +605,20 @@ def interpolate_integrator_quantities(

transport.electron_densities_integ = interp1d(
r_middle,
plasma.electron_densities,
plasma.electron_densities.iloc[
self.simulation_state.geometry.v_inner_boundary_index : self.simulation_state.geometry.v_outer_boundary_index
],
fill_value="extrapolate",
kind="nearest",
)(r_middle_integ)
# Assume tau_sobolevs to be constant within a shell
# (as in the MC simulation)
transport.tau_sobolevs_integ = interp1d(
r_middle,
self.opacity_state.tau_sobolev,
self.opacity_state.tau_sobolev[
:,
self.simulation_state.geometry.v_inner_boundary_index : self.simulation_state.geometry.v_outer_boundary_index,
],
fill_value="extrapolate",
kind="nearest",
)(r_middle_integ)
Expand Down
2 changes: 2 additions & 0 deletions tardis/workflows/simple_tardis_workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -268,6 +268,8 @@ def solve_simulation_state(
"t_inner"
]

return next_values

def solve_plasma(self, estimated_radfield_properties):
"""Update the plasma solution with the new radiation field estimates
Expand Down
Loading

0 comments on commit 5e2d0bb

Please sign in to comment.