Fix FPEs occurring when density is zero in MC

sxs-collaboration · Nov 21, 2024 · 8d28e53 · 8d28e53
1 parent f72ab25
commit 8d28e53
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 5 deletions.
diff --git a/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp b/src/Evolution/Particles/MonteCarlo/ImplicitMonteCarloCorrections.tpp
@@ -155,11 +155,13 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::
     // We have all we need to calculate beta here. We take the largest of the
     // two variations
     get(beta)[i] = fabs(get(neutrino_energy_1)[i] - get(neutrino_energy_0)[i]) /
-                   fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]);
+                   (fabs(get(fluid_energy_1)[i] - get(fluid_energy_0)[i]) +
+                    1.e-16);
     const double beta_lepton =
         fabs(get(neutrino_lepton_number_1)[i] -
              get(neutrino_lepton_number_0)[i]) /
-        fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]);
+        (fabs(get(fluid_lepton_number_1)[i] - get(fluid_lepton_number_0)[i]) +
+         1.e-16);
     get(beta)[i] =
         std::max(get(beta)[i], beta_lepton);
   }

diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp
@@ -309,9 +309,12 @@ void NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>::
   double ye = 0.0;
   double log_rho = 0.0;
   double log_temp = 0.0;
+  const double minimum_density = exp(table_log_density[0]);
   for (size_t p = 0; p < get(electron_fraction).size(); p++) {
     ye = get(electron_fraction)[p];
-    log_rho = log(get(rest_mass_density)[p]);
+    log_rho = get(rest_mass_density)[p] > minimum_density
+                  ? log(get(rest_mass_density)[p])
+                  : log(minimum_density);
     log_temp = get(temperature)[p] > minimum_temperature
                    ? log(get(temperature)[p])
                    : log(minimum_temperature);
@@ -351,6 +354,12 @@ void NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>::
             square(square(temperature_correction_factor));
         gsl::at(gsl::at(*scattering_opacity, ns), ng)[p] *=
             square(square(temperature_correction_factor));
+        gsl::at(gsl::at(*absorption_opacity, ns), ng)[p] =
+            std::clamp(gsl::at(gsl::at(*absorption_opacity, ns), ng)[p],
+                       min_kappa, max_kappa);
+        gsl::at(gsl::at(*scattering_opacity, ns), ng)[p] =
+            std::clamp(gsl::at(gsl::at(*scattering_opacity, ns), ng)[p],
+                       min_kappa, max_kappa);
       }
     }
   }

diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp
@@ -22,7 +22,6 @@ namespace{
 
 void test_flat_space_time_step(const bool skip) {
   CHECK(true);
-  // This test FPEs because the ghost zone handling is broken.
   if (skip) {
     return;
   }
@@ -270,5 +269,5 @@ void test_flat_space_time_step(const bool skip) {
 
 SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloTakeTimeStep",
                   "[Unit][Evolution]") {
-  test_flat_space_time_step(true);
+  test_flat_space_time_step(false);
 }