Skip to content

Commit

Permalink
modify to a neater way - shape the radiationfield to carry the t_rad …
Browse files Browse the repository at this point in the history
…and W of the active shells only
  • Loading branch information
DeerWhale committed Feb 6, 2024
1 parent dafdc28 commit 7e9bac7
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 39 deletions.
8 changes: 2 additions & 6 deletions tardis/model/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -144,16 +144,12 @@ def dilution_factor(self, new_dilution_factor):

@property
def t_radiative(self):
return self.radiation_field_state.t_radiative[
self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
]
return self.radiation_field_state.t_radiative

@t_radiative.setter
def t_radiative(self, new_t_radiative):
if len(new_t_radiative) == self.no_of_shells:
self.radiation_field_state.t_radiative[
self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
] = new_t_radiative
self.radiation_field_state.t_radiative = new_t_radiative
else:
raise ValueError(
"Trying to set t_radiative for different number of shells."
Expand Down
43 changes: 25 additions & 18 deletions tardis/model/parse_input.py
Original file line number Diff line number Diff line change
Expand Up @@ -555,26 +555,25 @@ def parse_radiation_field_state(
if t_radiative is None:
if config.plasma.initial_t_rad > 0 * u.K:
t_radiative = (
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
np.ones(geometry.no_of_shells_active)
* config.plasma.initial_t_rad
)
else:
t_radiative = calculate_t_radiative_from_t_inner(
geometry, packet_source
)

assert len(t_radiative) == geometry.no_of_shells
assert len(t_radiative) == geometry.no_of_shells_active

if dilution_factor is None:
dilution_factor = calculate_geometric_dilution_factor(geometry)
elif len(dilution_factor) != geometry.no_of_shells:
elif len(dilution_factor) != geometry.no_of_shells_active:
dilution_factor = dilution_factor[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
assert len(dilution_factor) == geometry.no_of_shells
assert len(dilution_factor) == geometry.no_of_shells_active

return DiluteThermalRadiationFieldState(
t_radiative, dilution_factor, geometry
)
return DiluteThermalRadiationFieldState(t_radiative, dilution_factor)


def initialize_packet_source(config, geometry, packet_source):
Expand Down Expand Up @@ -660,14 +659,18 @@ def parse_csvy_radiation_field_state(
t_rad_unit = u.Unit(
csvy_model_config.datatype.fields[t_rad_field_index]["unit"]
)
t_radiative = csvy_model_data["t_rad"].iloc[1:].values * t_rad_unit

elif config.plasma.initial_t_rad > 0 * u.K:
t_radiative = (
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
csvy_model_data["t_rad"][1:]
.iloc[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
.values
* t_rad_unit
)

elif config.plasma.initial_t_rad > 0 * u.K:
t_radiative = (
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
np.ones(geometry.no_of_shells_active) * config.plasma.initial_t_rad
)
else:
t_radiative = calculate_t_radiative_from_t_inner(
Expand All @@ -677,13 +680,17 @@ def parse_csvy_radiation_field_state(
if hasattr(csvy_model_data, "columns") and (
"dilution_factor" in csvy_model_data.columns
):
dilution_factor = csvy_model_data["dilution_factor"].iloc[1:].values
dilution_factor = (
csvy_model_data["dilution_factor"][1:]
.iloc[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
.values
)
else:
dilution_factor = calculate_geometric_dilution_factor(geometry)

return DiluteThermalRadiationFieldState(
t_radiative, dilution_factor, geometry
)
return DiluteThermalRadiationFieldState(t_radiative, dilution_factor)


def calculate_t_radiative_from_t_inner(geometry, packet_source):
Expand All @@ -705,7 +712,7 @@ def calculate_t_radiative_from_t_inner(geometry, packet_source):
lambda_wien_inner = const.b_wien / packet_source.temperature
t_radiative = const.b_wien / (
lambda_wien_inner
* (1 + (geometry.v_middle - geometry.v_inner_boundary) / const.c)
* (1 + (geometry.v_middle_active - geometry.v_inner_boundary) / const.c)
)
return t_radiative

Expand All @@ -715,7 +722,7 @@ def calculate_geometric_dilution_factor(geometry):
1
- np.sqrt(
1
- (geometry.r_inner_active[0] ** 2 / geometry.r_middle**2)
- (geometry.r_inner_active[0] ** 2 / geometry.r_middle_active**2)
.to(1)
.value
)
Expand Down
18 changes: 3 additions & 15 deletions tardis/model/radiation_field_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,23 +18,11 @@ class DiluteThermalRadiationFieldState:
Dilution Factors in each shell
"""

def __init__(
self, t_radiative: u.Quantity, dilution_factor: np.ndarray, geometry
):
def __init__(self, t_radiative: u.Quantity, dilution_factor: np.ndarray):
# ensuring that the radiation_field has both
# dilution_factor and t_radiative equal length
assert len(t_radiative) == len(dilution_factor)
assert np.all(
t_radiative[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
> 0 * u.K
)
assert np.all(
dilution_factor[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
> 0
)
assert np.all(t_radiative > 0 * u.K)
assert np.all(dilution_factor > 0)
self.t_radiative = t_radiative
self.dilution_factor = dilution_factor

0 comments on commit 7e9bac7

Please sign in to comment.