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

Update number of shells in SimulationState #2504

Merged
merged 14 commits into from
Feb 7, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
19 changes: 7 additions & 12 deletions tardis/model/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,9 @@ def dilution_factor(self):
@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 @@ -151,7 +153,9 @@ def t_radiative(self):
@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 Expand Up @@ -224,7 +228,7 @@ def volume(self):

@property
def no_of_shells(self):
return self.geometry.no_of_shells
return self.geometry.no_of_shells_active

@property
def no_of_raw_shells(self):
Expand Down Expand Up @@ -262,15 +266,6 @@ def from_config(cls, config, atom_data):
density, nuclide_mass_fraction, atom_data.atom_data.mass.copy()
)

nuclide_mass_fraction = parse_abundance_config(
config, geometry, time_explosion
)

# using atom_data.mass.copy() to ensure that the original atom_data is not modified
composition = Composition(
density, nuclide_mass_fraction, atom_data.atom_data.mass.copy()
)

Copy link
Contributor Author

@DeerWhale DeerWhale Jan 24, 2024

Choose a reason for hiding this comment

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

This modification is just to delete the duplicated codes

packet_source = parse_packet_source(config, geometry)
radiation_field_state = parse_radiation_field_state(
config,
Expand Down
30 changes: 17 additions & 13 deletions tardis/model/parse_input.py
Original file line number Diff line number Diff line change
Expand Up @@ -568,13 +568,12 @@ def parse_radiation_field_state(

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

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 @@ -610,13 +609,13 @@ def initialize_packet_source(config, geometry, packet_source):

luminosity_requested = config.supernova.luminosity_requested
if config.plasma.initial_t_inner > 0.0 * u.K:
packet_source.radius = geometry.r_inner[0]
packet_source.radius = geometry.r_inner_active[0]
packet_source.temperature = config.plasma.initial_t_inner

elif (config.plasma.initial_t_inner < 0.0 * u.K) and (
luminosity_requested is not None
):
packet_source.radius = geometry.r_inner[0]
packet_source.radius = geometry.r_inner_active[0]
packet_source.set_temperature_from_luminosity(luminosity_requested)
else:
raise ValueError(
Expand Down Expand Up @@ -674,9 +673,6 @@ def parse_csvy_radiation_field_state(
t_radiative = (
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
)
t_radiative = (
np.ones(geometry.no_of_shells) * config.plasma.initial_t_rad
)
else:
t_radiative = calculate_t_radiative_from_t_inner(
geometry, packet_source
Expand All @@ -689,7 +685,9 @@ def parse_csvy_radiation_field_state(
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 Down Expand Up @@ -720,6 +718,12 @@ def calculate_geometric_dilution_factor(geometry):
return 0.5 * (
1
- np.sqrt(
1 - (geometry.r_inner[0] ** 2 / geometry.r_middle**2).to(1).value
1
- (
geometry.r_inner[geometry.v_inner_boundary_index] ** 2
Copy link
Contributor

Choose a reason for hiding this comment

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

Could this be geometry.r_inner_active[0]?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Tried that but since r_inner_active[0] is overwrite by the v_inner_boundary value (which can be values in between the discrete shells), sometimes it can be a bigger value than the first element in r_middle, which cause the first w to be nan, and fails the assertion in the diluteBBradiationField.

/ geometry.r_middle**2
Copy link
Contributor

Choose a reason for hiding this comment

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

Double check that r_middle is correct here (matching the active radii)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

r_middle is the raw shells, but if we would like to keep the raw shape of the shells in radiationfield, then we should use r_middle I think.

)
.to(1)
.value
)
)
23 changes: 21 additions & 2 deletions tardis/model/radiation_field_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,18 +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,
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[
Copy link
Contributor

Choose a reason for hiding this comment

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

I know this fixes the issue, but we probably ought to be setting the estimators shape as the active shells from the start.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is actually already being taken of, since the estimators only use the active shells, the simulation_state slice the active shells from the "raw" radiation_field_state.

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
26 changes: 26 additions & 0 deletions tardis/model/tests/test_csvy_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,32 @@ def test_compare_models(model_config_fnames, atomic_dataset):
)


def test_dimensionality_after_update_v_inner_boundary(
example_csvy_file_dir, atomic_dataset
):
"""Test that the dimensionality of SimulationState parameters after updating v_inner_boundary
in the context of csvy models specifically"""
csvy_config_file = example_csvy_file_dir / "radiative_csvy.yml"
config = Configuration.from_yaml(csvy_config_file)
csvy_model = SimulationState.from_csvy(config, atom_data=atomic_dataset)

new_config = config
new_config.model.v_inner_boundary = csvy_model.velocity[1]
new_csvy_model = SimulationState.from_csvy(
new_config, atom_data=atomic_dataset
)

assert new_csvy_model.no_of_raw_shells == csvy_model.no_of_raw_shells
assert new_csvy_model.no_of_shells == csvy_model.no_of_shells - 1
assert new_csvy_model.velocity.shape[0] == csvy_model.velocity.shape[0] - 1
assert new_csvy_model.density.shape[0] == csvy_model.density.shape[0] - 1
assert new_csvy_model.volume.shape[0] == csvy_model.volume.shape[0] - 1
assert (
new_csvy_model.t_radiative.shape[0]
== csvy_model.t_radiative.shape[0] - 1
)


@pytest.fixture(scope="module")
def csvy_model_test_abundances(example_csvy_file_dir, atomic_dataset):
"""Returns SimulationState to use to test abundances dataframes"""
Expand Down
2 changes: 1 addition & 1 deletion tardis/visualization/tools/sdec_plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ def from_simulation(cls, sim, packets_mode):
"line_id"
)
transport_state = sim.transport.transport_state
r_inner = sim.simulation_state.geometry.r_inner
r_inner = sim.simulation_state.geometry.r_inner_active
t_inner = sim.simulation_state.packet_source.temperature
time_of_simulation = (
transport_state.packet_collection.time_of_simulation * u.s
Expand Down
Loading