Skip to content

Commit

Permalink
Gamma ray spectra dep (#2793)
Browse files Browse the repository at this point in the history 
* Changed GX_packet to accommodate positron_fraction

* Added tardis like string

* Added positronium sampler

* Added time evolution of decays as the composition changes

* Added 3D Doppler function for all packets

* changed name of root to discriminant-helps in understanding

* Added sampling frequencies and positronium fraction

* Modified energy deposition

* Added tests and changed main loop

* Added changes to get deposition

* Modified how times and effective times are used

* Changed mailmap

* Black formatting

* changed GXpacket attribute

* Fixed test

* Added a note to change appending to dataframes

* Changed the positron fraction and updated with a function to extract packet info

* Black fixes

* Added black changes

* Added ToDo notes

* Added corrections to comments

* Comments addressed

* Added review suggestions and small changes to the GX packet

* Changes to conftest

* Changes to pair_creation_packet by adding time_start

* Removed not necessary comments

* Removed experimental features from this PR

* Changed name to time_start

* Refurbished the positronium sampling

* added a new array for getting energy in photons by area by time

---------

Co-authored-by: Anirban Dutta <duttaan2@moria.egr.msu.edu>
  • Loading branch information
@Knights-Templars
Knights-Templars and Anirban Dutta authored Feb 26, 2025
1 parent e0297c6 commit a51af90
Show file tree
Hide file tree
Showing 14 changed files with 769 additions and 390 deletions.
1 change: 1 addition & 0 deletions .mailmap
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ Atharwa Kharkar <atharwa.kharkar19@vit.edu> atharwa_24 <atharwa.kharkar19@vit.ed
Atul Kumar <kr.atul.atk@gmail.com>

Anirban Dutta <anirbaniamdutta@gmail.com>
Anirban Dutta <duttaan2@msu.edu>

Barnabás Barna <bbarna@titan.physx.u-szeged.hu>

Expand Down
18 changes: 11 additions & 7 deletions tardis/energy_input/GXPacket.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ class GXPacketStatus(IntEnum):
("nu_cmf", float64),
("status", int64),
("shell", int64),
("time_current", float64),
("time_start", float64),
("time_index", int64),
("tau", float64),
]

Expand All @@ -52,7 +53,8 @@ def __init__(
nu_cmf,
status,
shell,
time_current,
time_start,
time_index,
):
self.location = location
self.direction = direction
Expand All @@ -62,7 +64,8 @@ def __init__(
self.nu_cmf = nu_cmf
self.status = status
self.shell = shell
self.time_current = time_current
self.time_start = time_start
self.time_index = time_index
# TODO: rename to tau_event
self.tau = -np.log(np.random.random())

Expand Down Expand Up @@ -94,7 +97,8 @@ def __init__(
nu_cmf,
status,
shell,
time_current,
time_start,
time_index,
):
self.location = location
self.direction = direction
Expand All @@ -104,9 +108,9 @@ def __init__(
self.nu_cmf = nu_cmf
self.status = status
self.shell = shell
self.time_current = time_current
self.number_of_packets = len(self.energy_rf)
self.tau = -np.log(np.random.random(self.number_of_packets))
self.time_start = time_start
self.time_index = time_index
self.tau = -np.log(np.random.random())


# @njit(**njit_dict_no_parallel)
Expand Down
100 changes: 41 additions & 59 deletions tardis/energy_input/gamma_packet_loop.py
View file
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,17 @@ def gamma_packet_loop(
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_out_cosi,
total_energy,
energy_deposited_gamma,
packets_info_array,
):
"""Propagates packets through the simulation
Expand Down Expand Up @@ -103,21 +102,20 @@ def gamma_packet_loop(
escaped_packets = 0
scattered_packets = 0
packet_count = len(packets)
# Logging does not work with numba. Using print instead.
print("Entering gamma ray loop for " + str(packet_count) + " packets")

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)
time_index = packet.time_index

if time_index < 0:
print(packet.time_current, time_index)
print(packet.time_start, 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 All @@ -129,7 +127,7 @@ def gamma_packet_loop(
doppler_factor = doppler_factor_3d(
packet.direction,
packet.location,
effective_time_array[time_index],
times[time_index],
)

kappa = kappa_calculation(comoving_energy)
Expand Down Expand Up @@ -210,21 +208,10 @@ def gamma_packet_loop(
distance_interaction, distance_boundary, distance_time
)

packet.time_current += distance / C_CGS
packet.time_start += distance / C_CGS

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 @@ -234,7 +221,7 @@ def gamma_packet_loop(
else:
packet.shell = get_index(
packet.get_location_r(),
inner_velocities * effective_time_array[time_index],
inner_velocities * times[time_index],
)

elif distance == distance_interaction:
Expand All @@ -246,47 +233,43 @@ def gamma_packet_loop(

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,
]
)
# Ejecta gains energy from the packets (gamma-rays)
energy_deposited_gamma[
packet.shell, time_index
] += ejecta_energy_gained
# Ejecta gains energy from both gamma-rays and positrons
total_energy[packet.shell, time_index] += ejecta_energy_gained

if packet.status == GXPacketStatus.PHOTOABSORPTION:
# Packet destroyed, go to the next packet
break
else:
packet.status = GXPacketStatus.IN_PROCESS
scattered = True
packet.status = GXPacketStatus.IN_PROCESS
scattered = True

else:
packet.shell += shell_change

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

# get energy out in ergs per second per keV
energy_out[bin_index, time_index] += (
packet.energy_rf
/ dt
/ freq_bin_width # Take light crossing time into account
)
# get energy out in photons per second per keV
energy_out_cosi[bin_index, time_index] += (1
/ dt
/ bin_width
)

luminosity = packet.energy_rf / dt
packet.status = GXPacketStatus.ESCAPED
escaped_packets += 1
if scattered:
Expand All @@ -303,22 +286,21 @@ def gamma_packet_loop(
packet.nu_cmf,
packet.nu_rf,
packet.energy_cmf,
lum_rf,
luminosity,
packet.energy_rf,
packet.shell,
]
)

print("Escaped packets:", escaped_packets)
print("Scattered packets:", scattered_packets)
print("Number of escaped packets:", escaped_packets)
print("Number of scattered packets:", scattered_packets)

return (
energy_df_rows,
energy_plot_df_rows,
energy_out,
deposition_estimator,
bin_width,
energy_out_cosi,
packets_info_array,
energy_deposited_gamma,
total_energy,
)


Expand Down Expand Up @@ -355,7 +337,7 @@ def process_packet_path(packet):
doppler_factor = doppler_factor_3d(
packet.direction,
packet.location,
packet.time_current,
packet.time_start,
)

packet.nu_rf = packet.nu_cmf / doppler_factor
Expand Down
Loading

0 comments on commit a51af90

Please sign in to comment.