From 77fa4657663e6d1cb77bd8c309e7e0c445ba46e0 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 10 May 2022 15:13:58 -0400 Subject: [PATCH 01/23] added v_outer and v_inner to NumbaModel --- tardis/montecarlo/montecarlo_numba/numba_interface.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index 4971db407a6..fc88f5650ea 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -17,12 +17,14 @@ ("r_inner", float64[:]), ("r_outer", float64[:]), ("time_explosion", float64), + ("v_inner", float64[:]), + ("v_outer", float64[:]), ] @jitclass(numba_model_spec) class NumbaModel(object): - def __init__(self, r_inner, r_outer, time_explosion): + def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): """ Model for the Numba mode @@ -30,11 +32,16 @@ def __init__(self, r_inner, r_outer, time_explosion): ---------- r_inner : numpy.ndarray r_outer : numpy.ndarray + v_inner : numpy.ndarray + v_outer : numpy.ndarray time_explosion : float """ self.r_inner = r_inner self.r_outer = r_outer + self.v_inner = v_inner + self.v_outer = v_outer self.time_explosion = time_explosion + numba_plasma_spec = [ From a7da24749f42dca92514535bf15bd80e9755a96b Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 10 May 2022 16:02:28 -0400 Subject: [PATCH 02/23] defined default value for v_inner and v_outer in NumbaModel --- tardis/montecarlo/montecarlo_numba/numba_interface.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index fc88f5650ea..8dd55995ff8 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -24,7 +24,8 @@ @jitclass(numba_model_spec) class NumbaModel(object): - def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): + def __init__(self, r_inner, r_outer, time_explosion): + # how to have a function that has default value for an object, and then gets assigned value afterwards? """ Model for the Numba mode @@ -38,8 +39,8 @@ def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): """ self.r_inner = r_inner self.r_outer = r_outer - self.v_inner = v_inner - self.v_outer = v_outer + self.v_inner = np.zeros(r_inner.shape) + self.v_outer = np.zeros(r_outer.shape) self.time_explosion = time_explosion From f6b4e3153cd6e44d6f1ddcf643f6bdb6f977e60c Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 10 May 2022 16:04:45 -0400 Subject: [PATCH 03/23] added enable_nonhomologous_expansion to io/schemas/montecarlo.yml --- tardis/io/schemas/montecarlo.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tardis/io/schemas/montecarlo.yml b/tardis/io/schemas/montecarlo.yml index 4a31f5228c5..d3811194e07 100644 --- a/tardis/io/schemas/montecarlo.yml +++ b/tardis/io/schemas/montecarlo.yml @@ -64,6 +64,10 @@ properties: default: false description: Enables a more complete treatment of relativitic effects. This includes angle aberration as well as use of the fully general Doppler formula. + enable_nonhomologous_expansion: + type: boolean + default: false + description: Enables nonhomologous expansion. Treats shells as piece-wise homologous areas for velocity-radius dependence. tracking: type: object default: {} From 816a1e800783a91649003c04d33797104f98293a Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Fri, 13 May 2022 13:32:38 -0400 Subject: [PATCH 04/23] added v_inner and v_outer to NumbaPlasma --- tardis/montecarlo/montecarlo_numba/numba_interface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index 8dd55995ff8..db2382840e7 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -24,7 +24,7 @@ @jitclass(numba_model_spec) class NumbaModel(object): - def __init__(self, r_inner, r_outer, time_explosion): + def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): # how to have a function that has default value for an object, and then gets assigned value afterwards? """ Model for the Numba mode @@ -39,8 +39,8 @@ def __init__(self, r_inner, r_outer, time_explosion): """ self.r_inner = r_inner self.r_outer = r_outer - self.v_inner = np.zeros(r_inner.shape) - self.v_outer = np.zeros(r_outer.shape) + self.v_inner = v_inner + self.v_outer = v_outer self.time_explosion = time_explosion From 13195eaf68597723256955670abcad5526308aac Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Fri, 13 May 2022 13:36:33 -0400 Subject: [PATCH 05/23] added v_inner and v_outer as arguments for NumbaModel in run() in montecarlo_numba --- tardis/montecarlo/montecarlo_numba/base.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py index be5b8f06e1e..e46aa547014 100644 --- a/tardis/montecarlo/montecarlo_numba/base.py +++ b/tardis/montecarlo/montecarlo_numba/base.py @@ -48,6 +48,8 @@ def montecarlo_radial1d( numba_model = NumbaModel( runner.r_inner_cgs, runner.r_outer_cgs, + runner.v_inner_cgs, + runner.v_outer_cgs, model.time_explosion.to("s").value, ) numba_plasma = numba_plasma_initialize(plasma, runner.line_interaction_type) From e6ad725881464bac0a0b36657b63ac9746941e90 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Fri, 13 May 2022 13:46:19 -0400 Subject: [PATCH 06/23] added v_outer_cgs to MonteCarloRunner --- tardis/montecarlo/base.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index ce3664ead89..641f1ae11a4 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -202,6 +202,8 @@ def _initialize_geometry_arrays(self, model): self.r_inner_cgs = model.r_inner.to("cm").value self.r_outer_cgs = model.r_outer.to("cm").value self.v_inner_cgs = model.v_inner.to("cm/s").value + self.v_outer_cgs = model.v_outer.to("cm/s").value + def _initialize_packets(self, T, no_of_packets, iteration, radius): # the iteration is added each time to preserve randomness From 8a1872761faca2ad73d552a2cd605eacf2296c34 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Fri, 13 May 2022 13:50:19 -0400 Subject: [PATCH 07/23] removed a comment --- tardis/montecarlo/montecarlo_numba/numba_interface.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index db2382840e7..fc88f5650ea 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -25,7 +25,6 @@ @jitclass(numba_model_spec) class NumbaModel(object): def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): - # how to have a function that has default value for an object, and then gets assigned value afterwards? """ Model for the Numba mode From 8d28731d227facfd133d7c692438be211983bcd0 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Fri, 13 May 2022 14:29:59 -0400 Subject: [PATCH 08/23] ran black on the changed files --- tardis/montecarlo/base.py | 13 ++++++------- .../montecarlo_numba/numba_interface.py | 17 ++++++----------- 2 files changed, 12 insertions(+), 18 deletions(-) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index 641f1ae11a4..cb9ef9607f6 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -29,7 +29,7 @@ logger = logging.getLogger(__name__) -MAX_SEED_VAL = 2**32 - 1 +MAX_SEED_VAL = 2 ** 32 - 1 # MAX_SEED_VAL must be multiple orders of magnitude larger than no_of_packets; # otherwise, each packet would not have its own seed. Here, we set the max @@ -73,14 +73,14 @@ class MontecarloRunner(HDFWriterMixin): hdf_name = "runner" w_estimator_constant = ( - (const.c**2 / (2 * const.h)) - * (15 / np.pi**4) + (const.c ** 2 / (2 * const.h)) + * (15 / np.pi ** 4) * (const.h / const.k_B) ** 4 / (4 * np.pi) ).cgs.value t_rad_estimator_constant = ( - (np.pi**4 / (15 * 24 * zeta(5, 1))) * (const.h / const.k_B) + (np.pi ** 4 / (15 * 24 * zeta(5, 1))) * (const.h / const.k_B) ).cgs.value def __init__( @@ -204,7 +204,6 @@ def _initialize_geometry_arrays(self, model): self.v_inner_cgs = model.v_inner.to("cm/s").value self.v_outer_cgs = model.v_outer.to("cm/s").value - def _initialize_packets(self, T, no_of_packets, iteration, radius): # the iteration is added each time to preserve randomness # across different simulations with the same temperature, @@ -563,7 +562,7 @@ def calculate_radiationfield_properties(self): w = self.j_estimator / ( 4 * const.sigma_sb.cgs.value - * t_rad**4 + * t_rad ** 4 * self.time_of_simulation.value * self.volume.value ) @@ -587,7 +586,7 @@ def calculate_luminosity_inner(self, model): * np.pi * const.sigma_sb.cgs * model.r_inner[0] ** 2 - * model.t_inner**4 + * model.t_inner ** 4 ).to("erg/s") def calculate_time_of_simulation(self, model): diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index fc88f5650ea..7a3d44440f9 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -41,7 +41,6 @@ def __init__(self, r_inner, r_outer, v_inner, v_outer, time_explosion): self.v_inner = v_inner self.v_outer = v_outer self.time_explosion = time_explosion - numba_plasma_spec = [ @@ -586,16 +585,12 @@ def configuration_initialize(runner, number_of_vpackets): montecarlo_configuration.full_relativity = runner.enable_full_relativity montecarlo_configuration.montecarlo_seed = runner.seed montecarlo_configuration.single_packet_seed = runner.single_packet_seed - montecarlo_configuration.v_packet_spawn_start_frequency = ( - runner.virtual_spectrum_spawn_range.end.to( - u.Hz, equivalencies=u.spectral() - ).value - ) - montecarlo_configuration.v_packet_spawn_end_frequency = ( - runner.virtual_spectrum_spawn_range.start.to( - u.Hz, equivalencies=u.spectral() - ).value - ) + montecarlo_configuration.v_packet_spawn_start_frequency = runner.virtual_spectrum_spawn_range.end.to( + u.Hz, equivalencies=u.spectral() + ).value + montecarlo_configuration.v_packet_spawn_end_frequency = runner.virtual_spectrum_spawn_range.start.to( + u.Hz, equivalencies=u.spectral() + ).value montecarlo_configuration.VPACKET_LOGGING = runner.virt_logging From cad480e5eda957dd3719055c47f1b4c7e3dd5200 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Mon, 16 May 2022 13:46:18 -0400 Subject: [PATCH 09/23] added right arguments in test --- tardis/montecarlo/montecarlo_numba/tests/conftest.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tardis/montecarlo/montecarlo_numba/tests/conftest.py b/tardis/montecarlo/montecarlo_numba/tests/conftest.py index 6b9be84d2ba..91d1ae63cfd 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/conftest.py +++ b/tardis/montecarlo/montecarlo_numba/tests/conftest.py @@ -52,6 +52,8 @@ def verysimple_numba_model(nb_simulation_verysimple): return NumbaModel( runner.r_inner_cgs, runner.r_outer_cgs, + runner.v_inner_cgs, + runner.v_outer_cgs, model.time_explosion.to("s").value, ) From 1649f22c45200e50a9864c631e5868f462ff8ff0 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Mon, 16 May 2022 15:30:32 -0400 Subject: [PATCH 10/23] added v_inner, v_outer to tests --- .../montecarlo_numba/r_packet_transport.py | 8 +++---- .../tests/test_cuda_formal_integral.py | 24 ++++--------------- .../tests/test_numba_formal_integral.py | 24 ++++--------------- .../montecarlo_numba/tests/test_packet.py | 11 ++++----- 4 files changed, 17 insertions(+), 50 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py index 6767a7e4771..cc57397c43e 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py @@ -44,10 +44,9 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators): r_inner = numba_model.r_inner[r_packet.current_shell_id] r_outer = numba_model.r_outer[r_packet.current_shell_id] - ( - distance_boundary, - delta_shell, - ) = calculate_distance_boundary(r_packet.r, r_packet.mu, r_inner, r_outer) + (distance_boundary, delta_shell,) = calculate_distance_boundary( + r_packet.r, r_packet.mu, r_inner, r_outer + ) # defining start for line interaction start_line_id = r_packet.next_line_id @@ -57,7 +56,6 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators): tau_trace_line_combined = 0.0 # e scattering initialization - cur_electron_density = numba_plasma.electron_density[ r_packet.current_shell_id ] diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py index 90a36010d91..0ce96db2128 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py @@ -31,12 +31,7 @@ not GPUs_available, reason="No GPU is available to test CUDA function" ) @pytest.mark.parametrize( - ["nu", "T"], - [ - (1e14, 1e4), - (0, 1), - (1, 1), - ], + ["nu", "T"], [(1e14, 1e4), (0, 1), (1, 1),], ) def test_intensity_black_body_cuda(nu, T): """ @@ -101,12 +96,8 @@ def trapezoid_integration_caller(data, h, actual): TESTDATA_model = [ - { - "r": np.linspace(1, 2, 3, dtype=np.float64), - }, - { - "r": np.linspace(0, 1, 3), - }, + {"r": np.linspace(1, 2, 3, dtype=np.float64),}, + {"r": np.linspace(0, 1, 3),}, # {"r": np.linspace(1, 2, 10, dtype=np.float64)}, ] @@ -117,7 +108,7 @@ def formal_integral_model(request): This gets the Numba model to be used in later tests """ r = request.param["r"] - model = NumbaModel(r[:-1], r[1:], 1 / c.c.cgs.value) + model = NumbaModel(r[:-1], r[1:], r[:-1], r[1:], 1 / c.c.cgs.value) return model @@ -207,12 +198,7 @@ def populate_z_caller( @pytest.mark.parametrize( - "N", - [ - 100, - 1000, - 10000, - ], + "N", [100, 1000, 10000,], ) def test_calculate_p_values(N): """ diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py index b019c787620..673e4d8e99c 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py @@ -11,12 +11,7 @@ @pytest.mark.parametrize( - ["nu", "T"], - [ - (1e14, 1e4), - (0, 1), - (1, 1), - ], + ["nu", "T"], [(1e14, 1e4), (0, 1), (1, 1),], ) def test_intensity_black_body(nu, T): func = formal_integral.intensity_black_body @@ -44,12 +39,8 @@ def calculate_z(r, p): TESTDATA = [ - { - "r": np.linspace(1, 2, 3, dtype=np.float64), - }, - { - "r": np.linspace(0, 1, 3), - }, + {"r": np.linspace(1, 2, 3, dtype=np.float64),}, + {"r": np.linspace(0, 1, 3),}, # {"r": np.linspace(1, 2, 10, dtype=np.float64)}, ] @@ -57,7 +48,7 @@ def calculate_z(r, p): @pytest.fixture(scope="function", params=TESTDATA) def formal_integral_model(request): r = request.param["r"] - model = NumbaModel(r[:-1], r[1:], 1 / c.c.cgs.value) + model = NumbaModel(r[:-1], r[1:], r[:-1], r[1:], 1 / c.c.cgs.value) return model @@ -152,12 +143,7 @@ def test_populate_z_shells(formal_integral_model, p): @pytest.mark.parametrize( - "N", - [ - 100, - 1000, - 10000, - ], + "N", [100, 1000, 10000,], ) def test_calculate_p_values(N): r = 1.0 diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index c227a87feff..e6f272feeb7 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -6,9 +6,7 @@ import tardis.montecarlo.montecarlo_numba.frame_transformations as frame_transformations import tardis.montecarlo.montecarlo_numba.opacities as opacities import tardis.montecarlo.montecarlo_numba.r_packet_transport as r_packet_transport -from tardis.montecarlo.montecarlo_numba.estimators import ( - update_line_estimators, -) +from tardis.montecarlo.montecarlo_numba.estimators import update_line_estimators import tardis.montecarlo.montecarlo_numba.utils as utils import tardis.montecarlo.montecarlo_numba.numba_interface as numba_interface @@ -31,6 +29,8 @@ def model(): return numba_interface.NumbaModel( r_inner=np.array([6.912e14, 8.64e14], dtype=np.float64), r_outer=np.array([8.64e14, 1.0368e15], dtype=np.float64), + v_inner=np.array([6.912e14, 8.64e14], dtype=np.float64), + v_outer=np.array([8.64e14, 1.0368e15], dtype=np.float64), time_explosion=5.2e7, ) @@ -224,10 +224,7 @@ def test_get_random_mu(set_seed_fixture): 0, 5.2e7, [[2.249675256109242, 0.0, 0.0], [0.0, 0.0, 0.0]], - [ - [2.249675256109242 * 0.4, 0.0, 1.0], - [0.0, 0.0, 1.0], - ], + [[2.249675256109242 * 0.4, 0.0, 1.0], [0.0, 0.0, 1.0],], ), ( 1, From b998d2643aad6421ed496f8aec2423d915e8e6ec Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Mon, 16 May 2022 15:59:17 -0400 Subject: [PATCH 11/23] added temporary fix to formal_integral.py --- tardis/montecarlo/montecarlo_numba/formal_integral.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tardis/montecarlo/montecarlo_numba/formal_integral.py b/tardis/montecarlo/montecarlo_numba/formal_integral.py index 404daec5249..f6e2f2494c1 100644 --- a/tardis/montecarlo/montecarlo_numba/formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/formal_integral.py @@ -292,6 +292,8 @@ def generate_numba_objects(self): self.numba_model = NumbaModel( self.runner.r_inner_i, self.runner.r_outer_i, + self.runner.r_inner_i / self.model.time_explosion.to("s").value, + self.runner.r_outer_i / self.model.time_explosion.to("s").value, self.model.time_explosion.to("s").value, ) self.numba_plasma = numba_plasma_initialize( From 32f8e18568fb9dc99588580d539aa6b04a82c26b Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Mon, 16 May 2022 16:17:48 -0400 Subject: [PATCH 12/23] ran black on modified files --- tardis/montecarlo/base.py | 12 +++++----- .../montecarlo_numba/numba_interface.py | 16 +++++++++----- .../montecarlo_numba/r_packet_transport.py | 7 +++--- .../tests/test_cuda_formal_integral.py | 22 +++++++++++++++---- .../tests/test_numba_formal_integral.py | 22 +++++++++++++++---- .../montecarlo_numba/tests/test_packet.py | 5 ++++- 6 files changed, 60 insertions(+), 24 deletions(-) diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py index cb9ef9607f6..b5f9d502440 100644 --- a/tardis/montecarlo/base.py +++ b/tardis/montecarlo/base.py @@ -29,7 +29,7 @@ logger = logging.getLogger(__name__) -MAX_SEED_VAL = 2 ** 32 - 1 +MAX_SEED_VAL = 2**32 - 1 # MAX_SEED_VAL must be multiple orders of magnitude larger than no_of_packets; # otherwise, each packet would not have its own seed. Here, we set the max @@ -73,14 +73,14 @@ class MontecarloRunner(HDFWriterMixin): hdf_name = "runner" w_estimator_constant = ( - (const.c ** 2 / (2 * const.h)) - * (15 / np.pi ** 4) + (const.c**2 / (2 * const.h)) + * (15 / np.pi**4) * (const.h / const.k_B) ** 4 / (4 * np.pi) ).cgs.value t_rad_estimator_constant = ( - (np.pi ** 4 / (15 * 24 * zeta(5, 1))) * (const.h / const.k_B) + (np.pi**4 / (15 * 24 * zeta(5, 1))) * (const.h / const.k_B) ).cgs.value def __init__( @@ -562,7 +562,7 @@ def calculate_radiationfield_properties(self): w = self.j_estimator / ( 4 * const.sigma_sb.cgs.value - * t_rad ** 4 + * t_rad**4 * self.time_of_simulation.value * self.volume.value ) @@ -586,7 +586,7 @@ def calculate_luminosity_inner(self, model): * np.pi * const.sigma_sb.cgs * model.r_inner[0] ** 2 - * model.t_inner ** 4 + * model.t_inner**4 ).to("erg/s") def calculate_time_of_simulation(self, model): diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py index 7a3d44440f9..3a9efd2c0ba 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_interface.py +++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py @@ -585,12 +585,16 @@ def configuration_initialize(runner, number_of_vpackets): montecarlo_configuration.full_relativity = runner.enable_full_relativity montecarlo_configuration.montecarlo_seed = runner.seed montecarlo_configuration.single_packet_seed = runner.single_packet_seed - montecarlo_configuration.v_packet_spawn_start_frequency = runner.virtual_spectrum_spawn_range.end.to( - u.Hz, equivalencies=u.spectral() - ).value - montecarlo_configuration.v_packet_spawn_end_frequency = runner.virtual_spectrum_spawn_range.start.to( - u.Hz, equivalencies=u.spectral() - ).value + montecarlo_configuration.v_packet_spawn_start_frequency = ( + runner.virtual_spectrum_spawn_range.end.to( + u.Hz, equivalencies=u.spectral() + ).value + ) + montecarlo_configuration.v_packet_spawn_end_frequency = ( + runner.virtual_spectrum_spawn_range.start.to( + u.Hz, equivalencies=u.spectral() + ).value + ) montecarlo_configuration.VPACKET_LOGGING = runner.virt_logging diff --git a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py index cc57397c43e..e12c8b803de 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py @@ -44,9 +44,10 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators): r_inner = numba_model.r_inner[r_packet.current_shell_id] r_outer = numba_model.r_outer[r_packet.current_shell_id] - (distance_boundary, delta_shell,) = calculate_distance_boundary( - r_packet.r, r_packet.mu, r_inner, r_outer - ) + ( + distance_boundary, + delta_shell, + ) = calculate_distance_boundary(r_packet.r, r_packet.mu, r_inner, r_outer) # defining start for line interaction start_line_id = r_packet.next_line_id diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py index 0ce96db2128..3541288537b 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py @@ -31,7 +31,12 @@ not GPUs_available, reason="No GPU is available to test CUDA function" ) @pytest.mark.parametrize( - ["nu", "T"], [(1e14, 1e4), (0, 1), (1, 1),], + ["nu", "T"], + [ + (1e14, 1e4), + (0, 1), + (1, 1), + ], ) def test_intensity_black_body_cuda(nu, T): """ @@ -96,8 +101,12 @@ def trapezoid_integration_caller(data, h, actual): TESTDATA_model = [ - {"r": np.linspace(1, 2, 3, dtype=np.float64),}, - {"r": np.linspace(0, 1, 3),}, + { + "r": np.linspace(1, 2, 3, dtype=np.float64), + }, + { + "r": np.linspace(0, 1, 3), + }, # {"r": np.linspace(1, 2, 10, dtype=np.float64)}, ] @@ -198,7 +207,12 @@ def populate_z_caller( @pytest.mark.parametrize( - "N", [100, 1000, 10000,], + "N", + [ + 100, + 1000, + 10000, + ], ) def test_calculate_p_values(N): """ diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py index 673e4d8e99c..bfe49b4cfa5 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py @@ -11,7 +11,12 @@ @pytest.mark.parametrize( - ["nu", "T"], [(1e14, 1e4), (0, 1), (1, 1),], + ["nu", "T"], + [ + (1e14, 1e4), + (0, 1), + (1, 1), + ], ) def test_intensity_black_body(nu, T): func = formal_integral.intensity_black_body @@ -39,8 +44,12 @@ def calculate_z(r, p): TESTDATA = [ - {"r": np.linspace(1, 2, 3, dtype=np.float64),}, - {"r": np.linspace(0, 1, 3),}, + { + "r": np.linspace(1, 2, 3, dtype=np.float64), + }, + { + "r": np.linspace(0, 1, 3), + }, # {"r": np.linspace(1, 2, 10, dtype=np.float64)}, ] @@ -143,7 +152,12 @@ def test_populate_z_shells(formal_integral_model, p): @pytest.mark.parametrize( - "N", [100, 1000, 10000,], + "N", + [ + 100, + 1000, + 10000, + ], ) def test_calculate_p_values(N): r = 1.0 diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index e6f272feeb7..6330dfaaf8e 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -224,7 +224,10 @@ def test_get_random_mu(set_seed_fixture): 0, 5.2e7, [[2.249675256109242, 0.0, 0.0], [0.0, 0.0, 0.0]], - [[2.249675256109242 * 0.4, 0.0, 1.0], [0.0, 0.0, 1.0],], + [ + [2.249675256109242 * 0.4, 0.0, 1.0], + [0.0, 0.0, 1.0], + ], ), ( 1, From 1b3b8d06af36f5e88b08a60c6626d56b9f925353 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 19 May 2022 22:46:26 -0400 Subject: [PATCH 13/23] changed dummy values in the test --- tardis/montecarlo/montecarlo_numba/tests/test_packet.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py index 6330dfaaf8e..0ad3c0e2e4a 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_packet.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_packet.py @@ -29,8 +29,8 @@ def model(): return numba_interface.NumbaModel( r_inner=np.array([6.912e14, 8.64e14], dtype=np.float64), r_outer=np.array([8.64e14, 1.0368e15], dtype=np.float64), - v_inner=np.array([6.912e14, 8.64e14], dtype=np.float64), - v_outer=np.array([8.64e14, 1.0368e15], dtype=np.float64), + v_inner=np.array([-1, -1], dtype=np.float64), + v_outer=np.array([-1, -1], dtype=np.float64), time_explosion=5.2e7, ) From 179937e628f351fb9c2d4307e83f2f6650e38763 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 19 May 2022 22:51:01 -0400 Subject: [PATCH 14/23] Reverted content of file. --- tardis/montecarlo/montecarlo_numba/r_packet_transport.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py index e12c8b803de..6767a7e4771 100644 --- a/tardis/montecarlo/montecarlo_numba/r_packet_transport.py +++ b/tardis/montecarlo/montecarlo_numba/r_packet_transport.py @@ -57,6 +57,7 @@ def trace_packet(r_packet, numba_model, numba_plasma, estimators): tau_trace_line_combined = 0.0 # e scattering initialization + cur_electron_density = numba_plasma.electron_density[ r_packet.current_shell_id ] From 0d85b6f890680df6dbe0e6eea4bebeb13e38d074 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 24 May 2022 14:17:14 -0400 Subject: [PATCH 15/23] made needed change in test_cuda_formal_integral.py --- .../montecarlo_numba/tests/test_cuda_formal_integral.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py index 3541288537b..9eaa077c22b 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py @@ -117,7 +117,13 @@ def formal_integral_model(request): This gets the Numba model to be used in later tests """ r = request.param["r"] - model = NumbaModel(r[:-1], r[1:], r[:-1], r[1:], 1 / c.c.cgs.value) + model = NumbaModel( + r[:-1], + r[1:], + r[:-1] * c.c.cgs.value, + r[1:] * c.c.cgs.value, + 1 / c.c.cgs.value, + ) return model From aa4049403bc617fc23fd96c1c55b1c10fd4f4f43 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 24 May 2022 14:20:49 -0400 Subject: [PATCH 16/23] made needed change in test_numba_formal_integral.py --- .../montecarlo_numba/tests/test_numba_formal_integral.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py index bfe49b4cfa5..2a83aa0e550 100644 --- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py +++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_formal_integral.py @@ -57,7 +57,13 @@ def calculate_z(r, p): @pytest.fixture(scope="function", params=TESTDATA) def formal_integral_model(request): r = request.param["r"] - model = NumbaModel(r[:-1], r[1:], r[:-1], r[1:], 1 / c.c.cgs.value) + model = NumbaModel( + r[:-1], + r[1:], + r[:-1] * c.c.cgs.value, + r[1:] * c.c.cgs.value, + 1 / c.c.cgs.value, + ) return model From 3927a7b1e99f607b50138887e4fc3a6b299c4fbf Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 26 May 2022 14:59:30 -0400 Subject: [PATCH 17/23] added nonhomologous_grid.py and enable_nonhomologous_expansion --- tardis/montecarlo/montecarlo_numba/numba_config.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tardis/montecarlo/montecarlo_numba/numba_config.py b/tardis/montecarlo/montecarlo_numba/numba_config.py index 4bceebf975c..132d206f2d1 100644 --- a/tardis/montecarlo/montecarlo_numba/numba_config.py +++ b/tardis/montecarlo/montecarlo_numba/numba_config.py @@ -8,3 +8,4 @@ H = const.h.cgs.value ENABLE_FULL_RELATIVITY = False +ENABLE_NONHOMOLOGOUS_EXPANSION = False From fdf63ed3823640408f1b2d51e2bedc4d55d068ff Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 26 May 2022 15:12:45 -0400 Subject: [PATCH 18/23] added nonhom_grip.py --- .../montecarlo_numba/nonhomologous_grid.py | 66 +++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py diff --git a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py new file mode 100644 index 00000000000..109bfc716ac --- /dev/null +++ b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py @@ -0,0 +1,66 @@ +import numpy as np +from numba import njit + +from tardis.montecarlo.montecarlo_numba import ( + njit_dict_no_parallel, +) + +from tardis.montecarlo.montecarlo_numba.numba_config import ENABLE_NONHOMOLOGOUS_EXPANSION + +@njit(**njit_dict_no_parallel) +def velocity(r_packet, shell_id, numba_model): + """ + Velocity at radius r + + Parameters + ---------- + r_packet: RPacket + shell_id: int + numba_model: NumbaModel + + Returns + ----------- + v: float, current velocity + """ + v_inner = numba_model.v_inner[shell_id] + v_outer = numba_model.v_outer[shell_id] + r_inner = numba_model.r_inner[shell_id] + r_outer = numba_model.r_outer[shell_id] + r = r_packet.r + frac = (v_outer-v_inner)/(r_outer-r_inner) + return v_inner + frac * (r - r_inner) + +@njit(**njit_dict_no_parallel) +def dvdr(shell_id, numba_model): + """ + dv/dr for the current shell + + Parameters + ---------- + shell_id: int + numba_model: NumbaModel + + Returns + ----------- + dvdr: float, dv/dr of the current shell + """ + v_inner = numba_model.v_inner[shell_id] + v_outer = numba_model.v_outer[shell_id] + r_inner = numba_model.r_inner[shell_id] + r_outer = numba_model.r_outer[shell_id] + return (v_outer-v_inner)/(r_outer-r_inner) + + +def tau_sobolev_factor(r_packet, shell_id, numba_model, numba_plasma): + if ENABLE_NONHOMOLOGOUS_EXPANSION: + v = velocity(r_packet, shell_id, numba_model) + r = r_packet.r + dvdr = dvdr(shell_id, numba_model) + mu = r_packet.mu + factor = 1.0/((1 - mu * mu)*v/r + mu * mu * dvdr) + + + else: + factor = np.ones(numba_plasma.tau_sobolev.shape) + return factor * numba_plasma.tau_sobolev + From 568551a945776ee58d7581cf24027a73eb1d09bc Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 31 May 2022 12:05:51 -0400 Subject: [PATCH 19/23] minor changes to nonhomologous_grip --- .../montecarlo_numba/nonhomologous_grid.py | 40 ++++++++++++------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py index 109bfc716ac..950c287c83b 100644 --- a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py +++ b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py @@ -8,20 +8,20 @@ from tardis.montecarlo.montecarlo_numba.numba_config import ENABLE_NONHOMOLOGOUS_EXPANSION @njit(**njit_dict_no_parallel) -def velocity(r_packet, shell_id, numba_model): +def velocity(r_packet, numba_model): """ Velocity at radius r Parameters ---------- r_packet: RPacket - shell_id: int numba_model: NumbaModel Returns ----------- v: float, current velocity """ + shell_id = r_packet.current_shell_id v_inner = numba_model.v_inner[shell_id] v_outer = numba_model.v_outer[shell_id] r_inner = numba_model.r_inner[shell_id] @@ -31,19 +31,20 @@ def velocity(r_packet, shell_id, numba_model): return v_inner + frac * (r - r_inner) @njit(**njit_dict_no_parallel) -def dvdr(shell_id, numba_model): +def dvdr(r_packet, numba_model): """ dv/dr for the current shell Parameters ---------- - shell_id: int + r_packet: RPacket numba_model: NumbaModel Returns ----------- dvdr: float, dv/dr of the current shell """ + shell_id = r_packet.current_shell_id v_inner = numba_model.v_inner[shell_id] v_outer = numba_model.v_outer[shell_id] r_inner = numba_model.r_inner[shell_id] @@ -51,16 +52,27 @@ def dvdr(shell_id, numba_model): return (v_outer-v_inner)/(r_outer-r_inner) -def tau_sobolev_factor(r_packet, shell_id, numba_model, numba_plasma): - if ENABLE_NONHOMOLOGOUS_EXPANSION: - v = velocity(r_packet, shell_id, numba_model) - r = r_packet.r - dvdr = dvdr(shell_id, numba_model) - mu = r_packet.mu - factor = 1.0/((1 - mu * mu)*v/r + mu * mu * dvdr) +@njit(**njit_dict_no_parallel) +def tau_sobolev_factor(r_packet, numba_model, numba_plasma): + """ + The angle and velocity dependent Tau Sobolev factor. Is called when ENABLE_NONHOMOLOGOUS_EXPANSION is set to True. + + Parameters + ---------- + r_packet: RPacket + numba_model: NumbaModel + numba_plasma: NumbaPlasma + Returns + ----------- + factor = //put the equation here + """ + shell_id = r_packet.current_shell_id - else: - factor = np.ones(numba_plasma.tau_sobolev.shape) - return factor * numba_plasma.tau_sobolev + v = velocity(r_packet, shell_id, numba_model) + r = r_packet.r + dvdr = dvdr(shell_id, numba_model) + mu = r_packet.mu + factor = 1.0/((1 - mu * mu)*v/r + mu * mu * dvdr) + return factor From 98c8536c109418badc0395e5e577a0a9a8301404 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Tue, 31 May 2022 12:13:17 -0400 Subject: [PATCH 20/23] made minor changes to nonhomogolous_grid.py --- tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py index 950c287c83b..1d5c7b83def 100644 --- a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py +++ b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py @@ -5,7 +5,6 @@ njit_dict_no_parallel, ) -from tardis.montecarlo.montecarlo_numba.numba_config import ENABLE_NONHOMOLOGOUS_EXPANSION @njit(**njit_dict_no_parallel) def velocity(r_packet, numba_model): From efd9170ba0612eaf04b533f1f4046590295a8c9d Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 9 Jun 2022 15:09:11 -0400 Subject: [PATCH 21/23] added root-finding function --- .../montecarlo_numba/nonhomologous_grid.py | 117 ++++++++++++++++++ 1 file changed, 117 insertions(+) diff --git a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py index 1d5c7b83def..d317a4e5fda 100644 --- a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py +++ b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py @@ -75,3 +75,120 @@ def tau_sobolev_factor(r_packet, numba_model, numba_plasma): factor = 1.0/((1 - mu * mu)*v/r + mu * mu * dvdr) return factor + +@njit(**njit_dict_no_parallel) +def roots(a, b, c, d, e, threshold): + ''' + Solves ax^4 + bx^3 + cx^2 + dx + e = 0, for the real roots greater than the threshold returns (x - threshold). + + Parameters + ----------- + a, b, c, d, e: coefficients of the equations ax^4 + bx^3 + cx^2 + dx + e = 0, float + threshold: lower needed limit on roots, float + + ''' + delta = 256*a**3*e**3 - 192*a**2*b*d*e**2 - 128*a**2*c**2*e**2 + 144*a**2*c*d**2*e - 27*a**2*d**4 + 144*a*b**2*c*e**2 - 6*a*b**2*d**2*e - 80*a*b*c**2*d*e + 18*a*b*c*d**3 + 16*a*c**4*e - 4*a*c**3*d**2 - 27*b**4*e**2 + 18*b**3*c*d*e - 4*b**3*d**3 - 4*b**2*c**3*e + b**2*c**2*d**2 + P = 8*a*c - 3*b**2 + R = b**3 + 8*d*a**2 - 4*a*b*c + delta_0 = c**2 - 3*b*d + 12*a*e + delta_1 = 2*c**3 - 9*b*c*d + 27*a*d**2 + 27*b**2*e -72*a*c*e + D = 64*a**3*e - 16*a**2*c**2 + 16*a*b**2*c - 16*a**2*b*d - 3*b**4 + p = (8*a*c - 3*b**2)/(8*a**2) + q = (b**3 - 4*a*b*c + 8*a**2*d)/(8*a**3) + # print("delta: ", delta) + # print("P: ", P) + # print("R: ", R) + # print("delta0: ", delta_0) + # print("D: ", D) + # print("p: ", p) + # print("q: ", q) + # print("delta1: ", delta_1) + x_1, x_2, x_3, x_4 = None, None, None, None + if (delta<0):#do care + + Q = np.cbrt((delta_1 + np.sqrt(delta_1**2 - 4 * delta_0**3))/2) + S = np.sqrt(-2/3 * p + 1/(3*a) * (Q + delta_0/Q))/2 + + if -4*S**2 - 2*p - q/S >= 0: + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_3 = None + x_4 = None + if -4*S**2 - 2*p + q/S>=0: + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_1 = None + x_2 = None + # print("two distinct real roots") + elif (delta>0): + if (P<0 and D<0): #do care + phi = np.arccos(delta_1/(2 * np.sqrt(delta_0**3))) + S = np.sqrt(-2/3 * p + 2/(3*a) *np.sqrt(delta_0) * np.cos(phi/3) )/2 + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + # print("four distinct real roots") + elif(P>0 and D>0): #do not care + + print("all complex") + else: + if (P<0 and D<0 and delta_0!= 0): #do care + phi = np.arccos(delta_1/(2 * np.sqrt(delta_0**3))) + S = np.sqrt(-2/3 * p + 2/(3*a) *np.sqrt(delta_0) * np.cos(phi/3) )/2 + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + # print("one real double root and two real simple roots") + elif(D>0) or (P>0 and (D!=0 or R!=0)): #do care + phi = np.arccos(delta_1/(2 * np.sqrt(delta_0**3))) + S = np.sqrt(-2/3 * p + 2/(3*a) *np.sqrt(delta_0) * np.cos(phi/3) )/2 + if -4*S**2 - 2*p - q/S >= 0: + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + if -4*S**2 - 2*p + q/S>=0: + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + # print("one real double root and two complex") + elif (delta_0 == 0 and D!=0): #do care + # print("one triple real root and one simple real root") + phi = np.arccos(delta_1/(2 * np.sqrt(delta_0**3))) + S = np.sqrt(-2/3 * p + 2/(3*a) *np.sqrt(delta_0) * np.cos(phi/3) )/2 + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + elif (D == 0): + if (P < 0): + print("two double real roots") #do care + phi = np.arccos(delta_1/(2 * np.sqrt(delta_0**3))) + S = np.sqrt(-2/3 * p + 2/(3*a) *np.sqrt(delta_0) * np.cos(phi/3) )/2 + x_1 = -b/(4*a) + S + 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_2 = -b/(4*a) + S - 1/2 * np.sqrt(-4*S**2 - 2*p - q/S) + x_3 = -b/(4*a) - S + 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + x_4 = -b/(4*a) - S - 1/2 * np.sqrt(-4*S**2 - 2*p + q/S) + elif (P > 0 and R == 0): + print("all complex") #don't care + elif (delta_0 == 0): + print("all four equal to: ", -b/(4*a)) #do care + x_1 = -b/(4*a) + roots = [] + if x_1 != None: + x_1 -= threshold + if x_1 > 0: + roots.append(x_1) + if x_2 != None: + x_2 -= threshold + if x_2 > 0: + roots.append(x_2) + if x_3 != None: + x_3 -= threshold + if x_3 > 0: + roots.append(x_3) + if x_4 != None: + x_4 -= threshold + if x_4 > 0: + roots.append(x_4) + + return roots From 1304fd124ac03cec99dfafdf29bfc2243ff5463d Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 9 Jun 2022 16:41:34 -0400 Subject: [PATCH 22/23] Added calculate_distance_line, still not tested --- .../montecarlo_numba/calculate_distances.py | 29 +++++++++++++++++++ .../montecarlo_numba/nonhomologous_grid.py | 6 ++-- 2 files changed, 32 insertions(+), 3 deletions(-) diff --git a/tardis/montecarlo/montecarlo_numba/calculate_distances.py b/tardis/montecarlo/montecarlo_numba/calculate_distances.py index 44f35a506e9..10ab28e4458 100644 --- a/tardis/montecarlo/montecarlo_numba/calculate_distances.py +++ b/tardis/montecarlo/montecarlo_numba/calculate_distances.py @@ -18,6 +18,11 @@ from tardis.montecarlo.montecarlo_numba.r_packet import ( print_r_packet_properties, ) +from tardis.montecarlo.montecarlo_numba import ( + velocity, + dv_dr, + roots, +) @njit(**njit_dict_no_parallel) @@ -145,3 +150,27 @@ def calculate_distance_electron(electron_density, tau_event): """ # add full_relativity here return tau_event / (electron_density * SIGMA_THOMSON) + + + +@njit(**njit_dict_no_parallel) +def calculate_distance_line_nonhomologous( + r_packet, comov_nu, is_last_line, nu_line, time_explosion, numba_model +): + dvdr = dv_dr(r_packet, numba_model) + nu = r_packet.nu + nu_ratio = (nu - nu_line)/nu + r_i = r_packet.r + mu_i = r_packet.mu + v_i = velocity(r_packet, numba_model) + + #coefficients for the quartic equation + a = dvdr**2 + b = -2 * C_SPEED_OF_LIGHT * nu_ratio * dvdr + c = C_SPEED_OF_LIGHT**2 * nu_ratio**2 + dvdr**2 * r_i**2 * (1 - mu_i**2) - (v_i - dvdr * r_i)**2 + d = -2 * C_SPEED_OF_LIGHT * nu_ratio * dvdr * r_i**2 * (1 - mu_i**2) + e = C_SPEED_OF_LIGHT**2 * nu_ratio**2 * r_i**2 * (1 - mu_i**2) + correction = r_i * mu_i + possible_l = roots(a, b, c, d, e, correction) + distance_line = min(possible_l) + return distance_line diff --git a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py index d317a4e5fda..e5809f36fc5 100644 --- a/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py +++ b/tardis/montecarlo/montecarlo_numba/nonhomologous_grid.py @@ -30,7 +30,7 @@ def velocity(r_packet, numba_model): return v_inner + frac * (r - r_inner) @njit(**njit_dict_no_parallel) -def dvdr(r_packet, numba_model): +def dv_dr(r_packet, numba_model): """ dv/dr for the current shell @@ -68,9 +68,9 @@ def tau_sobolev_factor(r_packet, numba_model, numba_plasma): """ shell_id = r_packet.current_shell_id - v = velocity(r_packet, shell_id, numba_model) + v = velocity(r_packet, numba_model) r = r_packet.r - dvdr = dvdr(shell_id, numba_model) + dvdr = dvdr(r_packet, numba_model) mu = r_packet.mu factor = 1.0/((1 - mu * mu)*v/r + mu * mu * dvdr) return factor From aafe50dc719e24ac648a103b701fe6d8e30454c7 Mon Sep 17 00:00:00 2001 From: Sona Chitchyan Date: Thu, 9 Jun 2022 17:07:37 -0400 Subject: [PATCH 23/23] minor bug --- tardis/montecarlo/montecarlo_numba/calculate_distances.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/montecarlo/montecarlo_numba/calculate_distances.py b/tardis/montecarlo/montecarlo_numba/calculate_distances.py index 10ab28e4458..1ff8972a152 100644 --- a/tardis/montecarlo/montecarlo_numba/calculate_distances.py +++ b/tardis/montecarlo/montecarlo_numba/calculate_distances.py @@ -18,7 +18,7 @@ from tardis.montecarlo.montecarlo_numba.r_packet import ( print_r_packet_properties, ) -from tardis.montecarlo.montecarlo_numba import ( +from tardis.montecarlo.montecarlo_numba.nonhomologous_grid import ( velocity, dv_dr, roots,