Add min/max energy methods to XsCalculator

src/celeritas/em/UrbanMscParams.cc

@@ -22,6 +22,7 @@
 #include "corecel/sys/ScopedMem.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/grid/PolyEvaluator.hh"
+#include "celeritas/grid/XsCalculator.hh"
 #include "celeritas/io/ImportData.hh"
 #include "celeritas/io/ImportProcess.hh"
 #include "celeritas/mat/MaterialParams.hh"
@@ -79,11 +80,6 @@ UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
     // Save electron mass
     host_data.electron_mass = particles.get(host_data.ids.electron).mass();
 
-    // Save high/low energy limits
-    auto const& grid = host_data.xs[ItemId<XsGridData>(0)].log_energy;
-    host_data.params.low_energy_limit = MevEnergy(std::exp(grid.front));
-    host_data.params.high_energy_limit = MevEnergy(std::exp(grid.back));
-
     // Coefficients for scaled Z
     static Array<double, 2> const a_coeff{{0.87, 0.70}};
     static Array<double, 2> const b_coeff{{2.0 / 3, 1.0 / 2}};
@@ -134,6 +130,13 @@ UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
                 == host_data.par_mat_data.size());
         }
     }
+
+    // Save high/low energy limits
+    XsCalculator calc_xs(host_data.xs[ItemId<XsGridData>(0)],
+                         make_const_ref(host_data).reals);
+    host_data.params.low_energy_limit = calc_xs.energy_min();
+    host_data.params.high_energy_limit = calc_xs.energy_max();
+
     CELER_ASSERT(host_data);
 
     // Move to mirrored data, copying to device

src/celeritas/em/WentzelVIMscParams.cc

@@ -13,6 +13,7 @@
 #include "corecel/io/Logger.hh"
 #include "corecel/sys/ScopedMem.hh"
 #include "celeritas/em/data/WentzelVIMscData.hh"
+#include "celeritas/grid/XsCalculator.hh"
 #include "celeritas/io/ImportData.hh"
 #include "celeritas/mat/MaterialParams.hh"
 #include "celeritas/phys/ParticleParams.hh"
@@ -65,9 +66,10 @@ WentzelVIMscParams::WentzelVIMscParams(ParticleParams const& particles,
     host_data.electron_mass = particles.get(host_data.ids.electron).mass();
 
     // Save high/low energy limits
-    auto const& grid = host_data.xs[ItemId<XsGridData>(0)].log_energy;
-    host_data.params.low_energy_limit = MevEnergy(std::exp(grid.front));
-    host_data.params.high_energy_limit = MevEnergy(std::exp(grid.back));
+    XsCalculator calc_xs(host_data.xs[ItemId<XsGridData>(0)],
+                         make_const_ref(host_data).reals);
+    host_data.params.low_energy_limit = calc_xs.energy_min();
+    host_data.params.high_energy_limit = calc_xs.energy_max();
 
     CELER_ASSERT(host_data);
 

src/celeritas/grid/XsCalculator.hh

@@ -56,9 +56,22 @@ class XsCalculator
     // Get the cross section at the given index
     inline CELER_FUNCTION real_type operator[](size_type index) const;
 
+    // Get the minimum energy
+    CELER_FUNCTION Energy energy_min() const
+    {
+        return Energy(std::exp(loge_grid_.front()));
+    }
+
+    // Get the maximum energy
+    CELER_FUNCTION Energy energy_max() const
+    {
+        return Energy(std::exp(loge_grid_.back()));
+    }
+
   private:
     XsGridData const& data_;
     Values const& reals_;
+    UniformGrid loge_grid_;
 
     CELER_FORCEINLINE_FUNCTION real_type get(size_type index) const;
 };
@@ -71,7 +84,7 @@ class XsCalculator
  */
 CELER_FUNCTION
 XsCalculator::XsCalculator(XsGridData const& grid, Values const& values)
-    : data_(grid), reals_(values)
+    : data_(grid), reals_(values), loge_grid_(data_.log_energy)
 {
     CELER_EXPECT(data_);
     CELER_ASSERT(grid.value.size() == data_.log_energy.size);
@@ -86,29 +99,28 @@ XsCalculator::XsCalculator(XsGridData const& grid, Values const& values)
  */
 CELER_FUNCTION real_type XsCalculator::operator()(Energy energy) const
 {
-    UniformGrid const loge_grid(data_.log_energy);
     real_type const loge = std::log(energy.value());
 
     // Snap out-of-bounds values to closest grid points
     size_type lower_idx;
     real_type result;
-    if (loge <= loge_grid.front())
+    if (loge <= loge_grid_.front())
     {
         lower_idx = 0;
         result = this->get(lower_idx);
     }
-    else if (loge >= loge_grid.back())
+    else if (loge >= loge_grid_.back())
     {
-        lower_idx = loge_grid.size() - 1;
+        lower_idx = loge_grid_.size() - 1;
         result = this->get(lower_idx);
     }
     else
     {
         // Locate the energy bin
-        lower_idx = loge_grid.find(loge);
-        CELER_ASSERT(lower_idx + 1 < loge_grid.size());
+        lower_idx = loge_grid_.find(loge);
+        CELER_ASSERT(lower_idx + 1 < loge_grid_.size());
 
-        real_type const upper_energy = std::exp(loge_grid[lower_idx + 1]);
+        real_type const upper_energy = std::exp(loge_grid_[lower_idx + 1]);
         real_type upper_xs = this->get(lower_idx + 1);
         if (lower_idx + 1 == data_.prime_index)
         {
@@ -119,7 +131,7 @@ CELER_FUNCTION real_type XsCalculator::operator()(Energy energy) const
 
         // Interpolate *linearly* on energy using the lower_idx data.
         LinearInterpolator<real_type> interpolate_xs(
-            {std::exp(loge_grid[lower_idx]), this->get(lower_idx)},
+            {std::exp(loge_grid_[lower_idx]), this->get(lower_idx)},
             {upper_energy, upper_xs});
         result = interpolate_xs(energy.value());
     }
@@ -137,8 +149,7 @@ CELER_FUNCTION real_type XsCalculator::operator()(Energy energy) const
  */
 CELER_FUNCTION real_type XsCalculator::operator[](size_type index) const
 {
-    UniformGrid const loge_grid(data_.log_energy);
-    real_type energy = std::exp(loge_grid[index]);
+    real_type energy = std::exp(loge_grid_[index]);
     real_type result = this->get(index);
 
     if (index >= data_.prime_index)

test/celeritas/grid/XsCalculator.test.cc

@@ -59,6 +59,10 @@ TEST_F(XsCalculatorTest, simple)
     // Test out-of-bounds
     EXPECT_SOFT_EQ(1.0, calc(Energy{0.0001}));
     EXPECT_SOFT_EQ(1e5, calc(Energy{1e7}));
+
+    // Test energy grid bounds
+    EXPECT_SOFT_EQ(1.0, value_as<Energy>(calc.energy_min()));
+    EXPECT_SOFT_EQ(1e5, value_as<Energy>(calc.energy_max()));
 }
 
 TEST_F(XsCalculatorTest, scaled_lowest)
@@ -89,6 +93,10 @@ TEST_F(XsCalculatorTest, scaled_lowest)
     // this might not be the best behavior for the lower energy value)
     EXPECT_SOFT_EQ(1000, calc(Energy{0.0001}));
     EXPECT_SOFT_EQ(0.1, calc(Energy{1e5}));
+
+    // Test energy grid bounds
+    EXPECT_SOFT_EQ(0.1, value_as<Energy>(calc.energy_min()));
+    EXPECT_SOFT_EQ(1e4, value_as<Energy>(calc.energy_max()));
 }
 
 TEST_F(XsCalculatorTest, scaled_middle)
@@ -127,6 +135,10 @@ TEST_F(XsCalculatorTest, scaled_middle)
     // this might not be the right behavior for
     EXPECT_SOFT_EQ(3, calc(Energy{0.0001}));
     EXPECT_SOFT_EQ(0.3, calc(Energy{1e5}));
+
+    // Test energy grid bounds
+    EXPECT_SOFT_EQ(0.1, value_as<Energy>(calc.energy_min()));
+    EXPECT_SOFT_EQ(1e4, value_as<Energy>(calc.energy_max()));
 }
 
 TEST_F(XsCalculatorTest, scaled_highest)
@@ -149,6 +161,10 @@ TEST_F(XsCalculatorTest, scaled_highest)
     // Final point and higher are scaled by 1/E
     EXPECT_SOFT_EQ(1, calc(Energy{100}));
     EXPECT_SOFT_EQ(.1, calc(Energy{1000}));
+
+    // Test energy grid bounds
+    EXPECT_SOFT_EQ(1, value_as<Energy>(calc.energy_min()));
+    EXPECT_SOFT_EQ(100, value_as<Energy>(calc.energy_max()));
 }
 
 TEST_F(XsCalculatorTest, TEST_IF_CELERITAS_DEBUG(scaled_off_the_end))