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 packet source refactor #2546

Merged
merged 25 commits into from
Apr 24, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
5b3351a
First steps
andrewfullard Mar 11, 2024
14f81c8
Add more methods to the packet source
andrewfullard Mar 13, 2024
c135284
Documentation, cleanup, positronium handling
andrewfullard Mar 13, 2024
ddd2c23
Black formatting, cleanup
andrewfullard Mar 13, 2024
4055329
Fix kappa_calculation import loop
andrewfullard Mar 18, 2024
f1edb89
Working packet source
andrewfullard Mar 18, 2024
4a0624d
Positron deposition and packet tracker arrays
andrewfullard Mar 18, 2024
72148c1
Fix parents dict
andrewfullard Mar 20, 2024
70a0dfe
Result now mostly consistent with master
andrewfullard Mar 25, 2024
2ee34f5
Fix kappa calculation test
andrewfullard Mar 25, 2024
485157a
Adds missing docstring for calculate_energy_factors
andrewfullard Mar 25, 2024
2ead93b
Fix util tests and black
andrewfullard Mar 25, 2024
cdb3b10
Testing skeleton for packet source
andrewfullard Mar 27, 2024
ac403be
black tests
andrewfullard Mar 27, 2024
88e9023
Address comments
andrewfullard Apr 3, 2024
45b0c34
Beginnings of an alternate packet source for the dataframe
andrewfullard Apr 3, 2024
0e9149a
Possible deposition fix, import fixes
andrewfullard Apr 9, 2024
0555b22
Added old functions
Knights-Templars Apr 9, 2024
4f66494
Cleanup, comments, and more dataframe sampling setup
andrewfullard Apr 9, 2024
33213d6
Corrected packet time computation
andrewfullard Apr 9, 2024
cc27c05
Modified sampling to use positron information
andrewfullard Apr 17, 2024
4dd6bc6
Remove packet sampling call and replace with dataframe method
andrewfullard Apr 22, 2024
0576977
Update method call to match current usage
andrewfullard Apr 22, 2024
25804e3
Comment responses
andrewfullard Apr 24, 2024
d4ff4b8
More comments, improve packet dataframe index handling
andrewfullard Apr 24, 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
30 changes: 30 additions & 0 deletions tardis/energy_input/GXPacket.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,36 @@ def get_location_r(self):
)


class GXPacketCollection:
"""
Gamma-ray packet collection
"""

def __init__(
self,
location,
direction,
energy_rf,
energy_cmf,
nu_rf,
nu_cmf,
status,
shell,
time_current,
):
self.location = location
self.direction = direction
self.energy_rf = energy_rf
self.energy_cmf = energy_cmf
self.nu_rf = nu_rf
self.nu_cmf = nu_cmf
self.status = status
self.shell = shell
self.time_current = time_current
Copy link
Member

Choose a reason for hiding this comment

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

What is time_current?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The current local time value of the packet, if there was a clock attached to it

self.number_of_packets = len(self.energy_rf)
self.tau = -np.log(np.random.random(self.number_of_packets))
Copy link
Member

Choose a reason for hiding this comment

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

Is this tau same as the mean life of isotopes? Shall we give a separate name for this variable if this the optical depth?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Given that we are in a radiative transfer code, I think we should rename the mean lifetime instead



# @njit(**njit_dict_no_parallel)
def initialize_packet_properties(
isotope_energy,
Expand Down
2 changes: 0 additions & 2 deletions tardis/energy_input/__init__.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
"""
Contains classes and functions to handle energy deposition and transport.
"""

from tardis.energy_input.util import *
7 changes: 1 addition & 6 deletions tardis/energy_input/gamma_packet_loop.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
photoabsorption_opacity_calculation,
pair_creation_opacity_calculation,
photoabsorption_opacity_calculation_kasen,
kappa_calculation,
pair_creation_opacity_artis,
SIGMA_T,
)
Expand All @@ -18,7 +19,6 @@
doppler_factor_3d,
C_CGS,
H_CGS_KEV,
kappa_calculation,
get_index,
)
from tardis.energy_input.GXPacket import GXPacketStatus
Expand Down Expand Up @@ -50,7 +50,6 @@ def gamma_packet_loop(
energy_df_rows,
energy_plot_df_rows,
energy_out,
packets_out,
packets_info_array,
):
"""Propagates packets through the simulation
Expand Down Expand Up @@ -289,9 +288,6 @@ def gamma_packet_loop(
energy_out[bin_index, time_index] += rest_energy / (
Copy link
Member

Choose a reason for hiding this comment

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

Instead of the energy_out, shall we dump the information into an array which which record the "escaping" packets information. That will help in getting the gamma-ray spectra in a post-processing step?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

energy_out is such an array, but we can record more packet information if you want. I consider this out of scope for this PR though.

bin_width * dt
)
packets_out[bin_index] = np.array(
[bin_index, i, rest_energy, packet.Z, packet.A]
)
packet.status = GXPacketStatus.ESCAPED
escaped_packets += 1
if scattered:
Expand Down Expand Up @@ -322,7 +318,6 @@ def gamma_packet_loop(
energy_plot_df_rows,
energy_out,
deposition_estimator,
packets_out,
bin_width,
Copy link
Member

Choose a reason for hiding this comment

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

Do we need the bin_width as an output? WE can do that later in a post-processing step.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Probably not, but the blame said you added it. I can remove it, but it may be out of scope for this PR.

packets_info_array,
)
Expand Down
2 changes: 1 addition & 1 deletion tardis/energy_input/gamma_ray_estimators.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,13 @@
compton_opacity_calculation,
SIGMA_T,
photoabsorption_opacity_calculation,
kappa_calculation,
)
from tardis.energy_input.util import (
angle_aberration_gamma,
doppler_factor_3d,
H_CGS_KEV,
ELECTRON_MASS_ENERGY_KEV,
kappa_calculation,
)


Expand Down
7 changes: 4 additions & 3 deletions tardis/energy_input/gamma_ray_interactions.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,12 @@
from numba import njit

from tardis.montecarlo.montecarlo_numba import njit_dict_no_parallel
from tardis.montecarlo.montecarlo_numba.opacities import compton_opacity_partial
from tardis.montecarlo.montecarlo_numba.opacities import (
compton_opacity_partial,
kappa_calculation,
)
from tardis.energy_input.util import (
get_random_unit_vector,
kappa_calculation,
euler_rodrigues,
compton_theta_distribution,
get_perpendicular_vector,
Expand Down Expand Up @@ -157,7 +159,6 @@ def get_compton_fraction_urilight(energy):

accept = False
while not accept:

z = np.random.random(3)
alpha1 = np.log(1.0 / x0)
alpha2 = (1.0 - x0**2.0) / 2.0
Expand Down
Loading
Loading