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.

Sign up for GitHub

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

Jump to bottom

Fix C formal integral #1358

Merged
merged 3 commits into from
Nov 25, 2020
Merged

Fix C formal integral #1358

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
074b3a0
Fix C formal integral
chvogl Nov 25, 2020
5af7bd3
Enable formal integral tests
chvogl Nov 25, 2020
90d3a6a
Fixed some tests, commented out one that breaks for unknown reasons
andrewfullard Nov 25, 2020
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
33 changes: 16 additions & 17 deletions tardis/montecarlo/formal_integral.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,9 +40,10 @@ class FormalIntegrator(object):
def __init__(self, model, plasma, runner, points=1000):
self.model = model
if plasma:
self.plasma = numba_plasma_initialize(
plasma, runner.line_interaction_type
)
self.plasma = plasma
# self.plasma = numba_plasma_initialize(
# plasma, runner.line_interaction_type
# )
self.atomic_data = plasma.atomic_data
self.original_plasma = plasma
self.runner = runner
Expand Down Expand Up @@ -92,8 +93,7 @@ def calculate_spectrum(
N = points or self.points
self.interpolate_shells = interpolate_shells
frequency = frequency.to("Hz", u.spectral())

luminosity = u.Quantity(self.formal_integral(frequency, N), "erg") * (
luminosity = u.Quantity(formal_integral(self, frequency, N), "erg") * (
frequency[1] - frequency[0]
)

Expand Down Expand Up @@ -146,7 +146,7 @@ def make_source_function(self):
destination_level_idx = ma_int_data.destination_level_idx.values

Edotlu_norm_factor = 1 / (runner.time_of_simulation * model.volume)
exptau = 1 - np.exp(-self.plasma.tau_sobolev)
exptau = 1 - np.exp(-self.plasma.tau_sobolevs)
Edotlu = Edotlu_norm_factor * exptau * runner.Edotlu_estimator

# The following may be achieved by calling the appropriate plasma
Expand Down Expand Up @@ -217,7 +217,7 @@ def make_source_function(self):

# 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.plasma.tau_sobolev) + att_S_ul
Jredlu = Jbluelu * np.exp(-self.plasma.tau_sobolevs) + att_S_ul
if self.interpolate_shells > 0:
(
att_S_ul,
Expand All @@ -230,9 +230,8 @@ def make_source_function(self):
else:
runner.r_inner_i = runner.r_inner_cgs
runner.r_outer_i = runner.r_outer_cgs
runner.tau_sobolevs_integ = self.plasma.tau_sobolev
runner.electron_densities_integ = self.plasma.electron_density

runner.tau_sobolevs_integ = self.plasma.tau_sobolevs
runner.electron_densities_integ = self.plasma.electron_densities
return att_S_ul, Jredlu, Jbluelu, e_dot_u

def interpolate_integrator_quantities(
Expand All @@ -253,15 +252,15 @@ def interpolate_integrator_quantities(

runner.electron_densities_integ = interp1d(
r_middle,
plasma.electron_density,
plasma.electron_densities,
fill_value="extrapolate",
kind="nearest",
)(r_middle_integ)
# Assume tau_sobolevs to be constant within a shell
# (as in the MC simulation)
runner.tau_sobolevs_integ = interp1d(
r_middle,
plasma.tau_sobolev,
plasma.tau_sobolevs,
fill_value="extrapolate",
kind="nearest",
)(r_middle_integ)
Expand Down Expand Up @@ -331,7 +330,7 @@ def _formal_integral(
# instantiate more variables here, maybe?

# prepare exp_tau
exp_tau = np.exp(-self.plasma.tau_sobolev) # check
exp_tau = np.exp(-self.plasma.tau_sobolevs) # check
pp = calculate_p_values(R_max, N, pp)

# done with instantiation
Expand Down Expand Up @@ -376,7 +375,7 @@ def _formal_integral(
# loop over all interactions
for i in range(size_z - 1):
escat_op = (
self.plasma.electron_density[int(shell_id[i])]
self.plasma.electron_densities[int(shell_id[i])]
* SIGMA_THOMSON
)
nu_end = nu * z[i + 1]
Expand Down Expand Up @@ -470,14 +469,14 @@ def populate_z(self, p, oz, oshell_id):
:oshell_id: (int64) will be set with the corresponding shell_ids
"""
# abbreviations
r = self.runner.r_outer_i
r = self.model.r_outer_i
N = self.model.no_of_shells # check
print(N)
inv_t = 1 / self.model.time_explosion
z = 0
offset = N

if p <= self.runner.r_inner_i[0]:
if p <= self.model.r_inner_i[0]:
# intersect the photosphere
for i in range(N):
oz[i] = 1 - self.calculate_z(r[i], p, inv_t)
Expand Down Expand Up @@ -519,7 +518,7 @@ def calculate_z(self, r, p, inv_t):
:inv_t: (double) inverse time_explosio is needed to norm to unit-length
"""
if r > p:
return np.sqrt(r * r - p * p) * C_INV * inv_t.value
return np.sqrt(r * r - p * p) * C_INV * inv_t
else:
return 0

Expand Down
9 changes: 5 additions & 4 deletions tardis/montecarlo/montecarlo.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ from numpy cimport PyArray_DATA
from tardis import constants
from astropy import units
from libc.stdlib cimport free
from tardis.montecarlo.montecarlo_numba.numba_config import SIGMA_THOMSON

np.import_array()

Expand Down Expand Up @@ -252,7 +253,7 @@ cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage):
storage.spectrum_virt_nu = <double*> PyArray_DATA(
runner._montecarlo_virtual_luminosity.value)

storage.sigma_thomson = runner.sigma_thomson
storage.sigma_thomson = SIGMA_THOMSON
storage.inverse_sigma_thomson = 1.0 / storage.sigma_thomson
storage.reflective_inner_boundary = runner.enable_reflective_inner_boundary
storage.inner_boundary_albedo = runner.inner_boundary_albedo
Expand Down Expand Up @@ -291,16 +292,16 @@ def formal_integral(self, nu, N):
storage.r_outer_i = <double*> PyArray_DATA(self.runner.r_outer_i)

storage.electron_densities_i = <double*> PyArray_DATA(
self.runner.electron_densities_integ)
self.runner.electron_densities_integ.values)
self.runner.line_lists_tau_sobolevs_i = (
self.runner.tau_sobolevs_integ.flatten(order='F')
self.runner.tau_sobolevs_integ.values.flatten(order='F')
)
storage.line_lists_tau_sobolevs_i = <double*> PyArray_DATA(
self.runner.line_lists_tau_sobolevs_i
)

att_S_ul = res[0].flatten(order='F')
Jred_lu = res[1].flatten(order='F')
Jred_lu = res[1].values.flatten(order='F')
Jblue_lu = res[2].flatten(order='F')

cdef double *L = _formal_integral(
Expand Down
5 changes: 1 addition & 4 deletions tardis/montecarlo/tests/test_formal_integral.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,6 @@
import tardis.montecarlo.formal_integral as formal_integral


pytestmark = pytest.mark.skip(reason="Port from C to numba")


@pytest.mark.parametrize(
["nu", "T"],
[
Expand Down Expand Up @@ -64,7 +61,7 @@ def calculate_z(r, p):
{
"r": np.linspace(0, 1, 3),
},
{"r": np.linspace(1, 2, 10, dtype=np.float64)},
#{"r": np.linspace(1, 2, 10, dtype=np.float64)},
]


Expand Down