Skip to content

Commit

Permalink
nvm, changing assertation criteria in diluteBB is better than change …
Browse files Browse the repository at this point in the history
…radiationfield to be active shells only, this way one can change the v_inner on fly
  • Loading branch information
DeerWhale committed Feb 7, 2024
1 parent 372e618 commit 3969d9c
Show file tree
Hide file tree
Showing 3 changed files with 53 additions and 35 deletions.
16 changes: 12 additions & 4 deletions tardis/model/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -128,12 +128,16 @@ def t_inner(self, value):

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

@dilution_factor.setter
def dilution_factor(self, new_dilution_factor):
if len(new_dilution_factor) == self.no_of_shells:
self.radiation_field_state.dilution_factor = new_dilution_factor
self.radiation_field_state.dilution_factor[
self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
] = new_dilution_factor
else:
raise ValueError(
"Trying to set dilution_factor for unmatching number"
Expand All @@ -142,12 +146,16 @@ def dilution_factor(self, new_dilution_factor):

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

@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 = new_t_radiative
self.radiation_field_state.t_radiative[
self.geometry.v_inner_boundary_index : self.geometry.v_outer_boundary_index
] = new_t_radiative
else:
raise ValueError(
"Trying to set t_radiative for different number of shells."
Expand Down
43 changes: 15 additions & 28 deletions tardis/model/parse_input.py
Original file line number Diff line number Diff line change
Expand Up @@ -557,25 +557,23 @@ 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_active)
* config.plasma.initial_t_rad
np.ones(geometry.no_of_shells) * 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_active
assert len(t_radiative) == geometry.no_of_shells

if dilution_factor is None:
dilution_factor = calculate_geometric_dilution_factor(geometry)
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_active

return DiluteBlackBodyRadiationFieldState(t_radiative, dilution_factor)
assert len(dilution_factor) == geometry.no_of_shells

return DiluteBlackBodyRadiationFieldState(
t_radiative, dilution_factor, geometry
)


def initialize_packet_source(config, geometry, packet_source):
Expand Down Expand Up @@ -669,18 +667,11 @@ 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"][1:]
.iloc[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
.values
* t_rad_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_active) * config.plasma.initial_t_rad
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
)
else:
t_radiative = calculate_t_radiative_from_t_inner(
Expand All @@ -690,17 +681,13 @@ 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"][1:]
.iloc[
geometry.v_inner_boundary_index : geometry.v_outer_boundary_index
]
.values
)
dilution_factor = csvy_model_data["dilution_factor"].iloc[1:].values
else:
dilution_factor = calculate_geometric_dilution_factor(geometry)

return DiluteBlackBodyRadiationFieldState(t_radiative, dilution_factor)
return DiluteBlackBodyRadiationFieldState(
t_radiative, dilution_factor, geometry
)


def calculate_t_radiative_from_t_inner(geometry, packet_source):
Expand All @@ -722,7 +709,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_active - geometry.v_inner_boundary) / const.c)
* (1 + (geometry.v_middle - geometry.v_inner_boundary) / const.c)
)
return t_radiative

Expand All @@ -734,7 +721,7 @@ def calculate_geometric_dilution_factor(geometry):
1
- (
geometry.r_inner[geometry.v_inner_boundary_index] ** 2
/ geometry.r_middle_active**2
/ geometry.r_middle**2
)
.to(1)
.value
Expand Down
29 changes: 26 additions & 3 deletions tardis/model/radiation_field_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,37 @@ class DiluteBlackBodyRadiationFieldState:
Radiative temperature in each shell
dilution_factor : numpy.ndarray
Dilution Factors in each shell
geometry: tardis.model.Radial1DModel
The geometry of the model that uses to constrains the active shells
"""

def __init__(self, t_radiative: u.Quantity, dilution_factor: np.ndarray):
def __init__(
self,
t_radiative: u.Quantity,
dilution_factor: np.ndarray,
geometry=None,
):
# 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 > 0 * u.K)
assert np.all(dilution_factor > 0)
if (
geometry is not None
): # check the active shells only (this is used when setting up the radiation_field_state)
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
)
else:
assert np.all(t_radiative > 0 * u.K)
assert np.all(dilution_factor > 0)
self.t_radiative = t_radiative
self.dilution_factor = dilution_factor

Expand Down

0 comments on commit 3969d9c

Please sign in to comment.