Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

Gamma ray spectra dep #2793

Open
wants to merge 29 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
ddc060c
Changed GX_packet to accommodate positron_fraction
Aug 5, 2024
cc6a1e5
Added tardis like string
Aug 5, 2024
1be314c
Added positronium sampler
Aug 5, 2024
82dc7de
Added time evolution of decays as the composition changes
Aug 5, 2024
cf0258d
Added 3D Doppler function for all packets
Aug 5, 2024
05f09f5
changed name of root to discriminant-helps in understanding
Aug 5, 2024
6e5f459
Added sampling frequencies and positronium fraction
Aug 5, 2024
8d2179d
Modified energy deposition
Aug 5, 2024
0ab1571
Added tests and changed main loop
Aug 5, 2024
6dd1c54
Added changes to get deposition
Aug 13, 2024
ca8c054
Modified how times and effective times are used
Aug 19, 2024
9e8c19a
Changed mailmap
Aug 19, 2024
101c037
Black formatting
Aug 19, 2024
8484acd
changed GXpacket attribute
Aug 19, 2024
1d0398e
Fixed test
Aug 19, 2024
76115ad
Added a note to change appending to dataframes
Aug 19, 2024
906935e
Changed the positron fraction and updated with a function to extract …
Knights-Templars Aug 26, 2024
099f7a7
Black fixes
Knights-Templars Aug 26, 2024
bc5c0c6
Added black changes
Knights-Templars Aug 26, 2024
9f83760
Added ToDo notes
Knights-Templars Aug 26, 2024
8d39eb0
Merged updates into working branch
Knights-Templars Nov 25, 2024
0a5ac59
Added corrections to comments
Knights-Templars Dec 2, 2024
8fb6f48
Comments addressed
Knights-Templars Dec 2, 2024
cd6e099
Added review suggestions and small changes to the GX packet
Knights-Templars Dec 2, 2024
112f98a
Changes to conftest
Knights-Templars Dec 2, 2024
984b95e
Changes to pair_creation_packet by adding time_start
Knights-Templars Dec 2, 2024
0f9976e
Removed not necessary comments
Knights-Templars Dec 3, 2024
923b98c
Removed experimental features from this PR
Knights-Templars Dec 4, 2024
00a34a0
Changed name to time_start
Knights-Templars Dec 4, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .mailmap
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ Atharwa Kharkar <[email protected]> atharwa_24 <[email protected]
Atul Kumar <[email protected]>

Anirban Dutta <[email protected]>
Anirban Dutta <[email protected]>
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change mailmap!


Barnabás Barna <[email protected]>

Expand Down
5 changes: 5 additions & 0 deletions tardis/energy_input/GXPacket.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
("status", int64),
("shell", int64),
("time_current", float64),
("positron_energy", float64),
("tau", float64),
]

Expand All @@ -53,6 +54,7 @@
status,
shell,
time_current,
positron_energy,
):
self.location = location
self.direction = direction
Expand All @@ -64,6 +66,7 @@
self.shell = shell
self.time_current = time_current
# TODO: rename to tau_event
self.positron_energy = positron_energy

Check warning on line 69 in tardis/energy_input/GXPacket.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/GXPacket.py#L69

Added line #L69 was not covered by tests
self.tau = -np.log(np.random.random())

def get_location_r(self):
Expand Down Expand Up @@ -95,6 +98,7 @@
status,
shell,
time_current,
positron_energy,
):
self.location = location
self.direction = direction
Expand All @@ -106,6 +110,7 @@
self.shell = shell
self.time_current = time_current
self.number_of_packets = len(self.energy_rf)
self.positron_energy = positron_energy

Check warning on line 113 in tardis/energy_input/GXPacket.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/GXPacket.py#L113

Added line #L113 was not covered by tests
self.tau = -np.log(np.random.random(self.number_of_packets))


Expand Down
70 changes: 22 additions & 48 deletions tardis/energy_input/gamma_packet_loop.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import numpy as np
from numba import njit

from tardis.energy_input.gamma_ray_estimators import deposition_estimator_kasen
from tardis.energy_input.gamma_ray_grid import (
distance_trace,
move_packet,
Expand Down Expand Up @@ -38,17 +37,16 @@
pair_creation_opacity_type,
electron_number_density_time,
mass_density_time,
inv_volume_time,
iron_group_fraction_per_shell,
inner_velocities,
outer_velocities,
times,
dt_array,
times,
effective_time_array,
energy_bins,
energy_df_rows,
energy_plot_df_rows,
energy_out,
energy_deposited_gamma,
energy_deposited_positron,
packets_info_array,
):
"""Propagates packets through the simulation
Expand Down Expand Up @@ -103,21 +101,23 @@
escaped_packets = 0
scattered_packets = 0
packet_count = len(packets)
# Logging does not work with numba
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove the "print statement".

print("Entering gamma ray loop for " + str(packet_count) + " packets")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would much prefer using logging, and just doing it right above the function call rather than inside the function.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep now this is moving from being in-development, we should formalize some of these logs (which has already been done in some of the other major functions)


deposition_estimator = np.zeros_like(energy_df_rows)

for i in range(packet_count):
packet = packets[i]
time_index = get_index(packet.time_current, times)
energy_deposited_positron[

Check warning on line 110 in tardis/energy_input/gamma_packet_loop.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_packet_loop.py#L110

Added line #L110 was not covered by tests
packet.shell, time_index
] += packet.positron_energy

andrewfullard marked this conversation as resolved.
Show resolved Hide resolved
if time_index < 0:
print(packet.time_current, time_index)
raise ValueError("Packet time index less than 0!")

scattered = False

initial_energy = packet.energy_cmf
# Not used now. Useful for the deposition estimator.
# initial_energy = packet.energy_cmf

while packet.status == GXPacketStatus.IN_PROCESS:
# Get delta-time value for this step
Expand Down Expand Up @@ -214,17 +214,6 @@

packet = move_packet(packet, distance)

deposition_estimator[packet.shell, time_index] += (
(initial_energy * 1000)
* distance
* (packet.energy_cmf / initial_energy)
* deposition_estimator_kasen(
comoving_energy,
mass_density_time[packet.shell, time_index],
iron_group_fraction_per_shell[packet.shell],
)
)

if distance == distance_time:
time_index += 1

Expand All @@ -246,26 +235,9 @@

packet, ejecta_energy_gained = process_packet_path(packet)

# Save packets to dataframe rows
# convert KeV to eV / s / cm^3
energy_df_rows[packet.shell, time_index] += (
ejecta_energy_gained * 1000
)

energy_plot_df_rows[i] = np.array(
[
i,
ejecta_energy_gained * 1000
# * inv_volume_time[packet.shell, time_index]
/ dt,
packet.get_location_r(),
packet.time_current,
packet.shell,
compton_opacity,
photoabsorption_opacity,
pair_creation_opacity,
]
)
energy_deposited_gamma[

Check warning on line 238 in tardis/energy_input/gamma_packet_loop.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_packet_loop.py#L238

Added line #L238 was not covered by tests
packet.shell, time_index
] += ejecta_energy_gained

if packet.status == GXPacketStatus.PHOTOABSORPTION:
# Packet destroyed, go to the next packet
Expand All @@ -279,14 +251,18 @@

if packet.shell > len(mass_density_time[:, 0]) - 1:
rest_energy = packet.nu_rf * H_CGS_KEV
lum_rf = (packet.energy_rf * 1.6022e-9) / dt
bin_index = get_index(rest_energy, energy_bins)
bin_width = (
energy_bins[bin_index + 1] - energy_bins[bin_index]
)
energy_out[bin_index, time_index] += rest_energy / (
bin_width * dt
freq_bin_width = bin_width / H_CGS_KEV
energy_out[bin_index, time_index] += (

Check warning on line 259 in tardis/energy_input/gamma_packet_loop.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_packet_loop.py#L258-L259

Added lines #L258 - L259 were not covered by tests
packet.energy_rf
/ dt
/ freq_bin_width # Take light crossing time into account
)

luminosity = packet.energy_rf / dt

Check warning on line 265 in tardis/energy_input/gamma_packet_loop.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_packet_loop.py#L265

Added line #L265 was not covered by tests
packet.status = GXPacketStatus.ESCAPED
escaped_packets += 1
if scattered:
Expand All @@ -303,7 +279,7 @@
packet.nu_cmf,
packet.nu_rf,
packet.energy_cmf,
lum_rf,
luminosity,
packet.energy_rf,
packet.shell,
]
Expand All @@ -313,12 +289,10 @@
print("Scattered packets:", scattered_packets)

return (
energy_df_rows,
energy_plot_df_rows,
energy_out,
deposition_estimator,
bin_width,
packets_info_array,
energy_deposited_gamma,
energy_deposited_positron,
)


Expand Down
93 changes: 91 additions & 2 deletions tardis/energy_input/gamma_ray_channel.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import radioactivedecay as rd

from tardis.energy_input.util import KEV2ERG
from tardis.model.matter.decay import IsotopicMassFraction

logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
Expand Down Expand Up @@ -83,10 +84,10 @@
cumulative_decay_df : pd.DataFrame
total decays for x g of isotope for time 't'
"""
time_delta = u.Quantity(time_delta, u.s)
time_delta = u.Quantity(time_delta, u.d)
total_decays = {}
for shell, inventory in inventories.items():
total_decays[shell] = inventory.cumulative_decays(time_delta.value)
total_decays[shell] = inventory.cumulative_decays(time_delta.value, "d")

flattened_dict = {}

Expand Down Expand Up @@ -167,3 +168,91 @@
)

return isotope_decay_df


def time_evolve_mass_fraction(raw_isotope_mass_fraction, time_array):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe combine this with the time_evolve_cumilative_decay function below, and just have the combined function return both the new mass fractions and cumulative decays.

"""
Function to evolve the mass fraction of isotopes with time.

Parameters
----------
raw_isotope_mass_fraction : pd.DataFrame
isotope mass fraction in mass fractions.
time_array : np.array
array of time in days.

Returns
-------
time_evolved_isotope_mass_fraction : pd.DataFrame
time evolved mass fraction of isotopes.
"""

initial_isotope_mass_fraction = raw_isotope_mass_fraction
isotope_mass_fraction_list = []

Check warning on line 191 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L190-L191

Added lines #L190 - L191 were not covered by tests

for time in time_array:

Check warning on line 193 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L193

Added line #L193 was not covered by tests

decayed_isotope_mass_fraction = IsotopicMassFraction(

Check warning on line 195 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L195

Added line #L195 was not covered by tests
initial_isotope_mass_fraction
).decay(time)
isotope_mass_fraction_list.append(decayed_isotope_mass_fraction)
initial_isotope_mass_fraction = decayed_isotope_mass_fraction

Check warning on line 199 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L198-L199

Added lines #L198 - L199 were not covered by tests

time_evolved_isotope_mass_fraction = pd.concat(

Check warning on line 201 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L201

Added line #L201 was not covered by tests
isotope_mass_fraction_list, keys=time_array, names=["time"]
)

return time_evolved_isotope_mass_fraction

Check warning on line 205 in tardis/energy_input/gamma_ray_channel.py

View check run for this annotation

Codecov / codecov/patch

tardis/energy_input/gamma_ray_channel.py#L205

Added line #L205 was not covered by tests


def time_evolve_cumulative_decay(
raw_isotope_mass_fraction, shell_masses, gamma_ray_lines, time_array
):
"""
Function to calculate the total decays for each isotope for each shell at each time step.

Parameters
----------
raw_isotope_mass_fraction : pd.DataFrame
isotope abundance in mass fractions.
shell_masses : numpy.ndarray
shell masses in units of g
gamma_ray_lines : pd.DataFrame
gamma ray lines from nndc stored as a pandas dataframe.
time_array : numpy.ndarray
array of time steps in days.

Returns
-------
time_evolve_decay_df : pd.DataFrame
dataframe of isotopes for each shell with their decay mode, number of decays, radiation type,
radiation energy and radiation intensity at each time step.

"""
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function duplicates some code from the function above that could be factored out into a different function


isotope_decay_df_list = []
initial_isotope_mass_fraction = raw_isotope_mass_fraction

dt = np.diff(time_array)
# Do not append dataframes to empty lists
for time in dt:
#
isotope_dict = create_isotope_dicts(
initial_isotope_mass_fraction, shell_masses
)
inventories = create_inventories_dict(isotope_dict)
total_decays = calculate_total_decays(inventories, time)
isotope_df_time = create_isotope_decay_df(total_decays, gamma_ray_lines)
isotope_decay_df_list.append(isotope_df_time)

decayed_isotope_mass_fraction = IsotopicMassFraction(
initial_isotope_mass_fraction
).decay(time)

initial_isotope_mass_fraction = decayed_isotope_mass_fraction

time_evolved_decay_df = pd.concat(
isotope_decay_df_list, keys=time_array, names=["time"]
)

return time_evolved_decay_df
Loading
Loading