From ac6cd34850471e5138f4b1edadc1ba865633fcb1 Mon Sep 17 00:00:00 2001 From: Alexander Holas Date: Wed, 1 Nov 2023 15:38:26 -0400 Subject: [PATCH 1/4] assertion for first guess --- tardis/plasma/properties/nlte_rate_equation_solver.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tardis/plasma/properties/nlte_rate_equation_solver.py b/tardis/plasma/properties/nlte_rate_equation_solver.py index a2e8b5f028f..e3be83c383c 100644 --- a/tardis/plasma/properties/nlte_rate_equation_solver.py +++ b/tardis/plasma/properties/nlte_rate_equation_solver.py @@ -110,6 +110,10 @@ def calculate( number_density[shell], initial_electron_densities[shell], ) + # All first guess values have to be positive + assert ( + np.greater_equal(first_guess, 0.0).all() + ).all(), "First guess for NLTE solver has negative values, something went wrong." solution = root( self.population_objective_function, first_guess, From 2a4c3d8e950b3eccc5e4905bd2ff524dc11a4974 Mon Sep 17 00:00:00 2001 From: Alexander Holas Date: Wed, 1 Nov 2023 16:10:21 -0400 Subject: [PATCH 2/4] Add checks to decay --- tardis/io/decay.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/tardis/io/decay.py b/tardis/io/decay.py index a61a54a3cb8..c838c9fb4f5 100644 --- a/tardis/io/decay.py +++ b/tardis/io/decay.py @@ -2,10 +2,12 @@ from radioactivedecay import Nuclide, Inventory from radioactivedecay.utils import Z_to_elem from astropy import units as u +import logging +logger = logging.getLogger(__name__) -class IsotopeAbundances(pd.DataFrame): +class IsotopeAbundances(pd.DataFrame): _metadata = ["time_0"] def __init__(self, *args, **kwargs): @@ -92,9 +94,18 @@ def decay(self, t): t_second = ( u.Quantity(t, u.day).to(u.s).value - self.time_0.to(u.s).value ) + logger.info(f"Decaying abundances for {t_second} seconds") + if t_second < 0: + logger.warning( + f"Decay time {t_second} is negative. This could indicate a miss-specified input model." + f" A negative decay time can potentially lead to negative abundances." + ) decayed_inventories = [item.decay(t_second) for item in inventories] df = IsotopeAbundances.from_inventories(decayed_inventories) df.sort_index(inplace=True) + assert ( + df.gt(0).all().all() + ), "Negative abundances detected. Please make sure your input abundances are correct." return df def as_atoms(self): From 7a738622facdede6b57bb98070045b972a1060da Mon Sep 17 00:00:00 2001 From: Alexander Holas Date: Wed, 1 Nov 2023 16:47:13 -0400 Subject: [PATCH 3/4] Fix greater equal --- tardis/io/decay.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/io/decay.py b/tardis/io/decay.py index c838c9fb4f5..7005fa50e6e 100644 --- a/tardis/io/decay.py +++ b/tardis/io/decay.py @@ -104,7 +104,7 @@ def decay(self, t): df = IsotopeAbundances.from_inventories(decayed_inventories) df.sort_index(inplace=True) assert ( - df.gt(0).all().all() + df.ge(0.0).all().all() ), "Negative abundances detected. Please make sure your input abundances are correct." return df From b98d02bce2a8a98cb8d72855df56f1af7ae48702 Mon Sep 17 00:00:00 2001 From: Alexander Holas Date: Wed, 1 Nov 2023 16:57:25 -0400 Subject: [PATCH 4/4] Fix deprecated series index --- tardis/plasma/properties/nlte_rate_equation_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/plasma/properties/nlte_rate_equation_solver.py b/tardis/plasma/properties/nlte_rate_equation_solver.py index e3be83c383c..487b93ad5a3 100644 --- a/tardis/plasma/properties/nlte_rate_equation_solver.py +++ b/tardis/plasma/properties/nlte_rate_equation_solver.py @@ -595,7 +595,7 @@ def prepare_first_guess( """ first_guess = pd.Series(0.0, index=rate_matrix_index) for atomic_number in atomic_numbers: - first_guess.loc[(atomic_number, 1)][0] = number_density.loc[ + first_guess.loc[(atomic_number, 1)].iloc[0] = number_density.loc[ atomic_number ] # TODO: After the first iteration, the new guess can be the old solution.