Use ValueGridInserter and support linear interpolation

src/celeritas/grid/RangeGridCalculator.cc

@@ -14,11 +14,15 @@
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Default constructor.
+ */
+RangeGridCalculator::RangeGridCalculator() : RangeGridCalculator(BC::size_) {}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct with boundary conditions for spline interpolation.
- *
- * \todo support polynomial interpolation as well?
  */
 RangeGridCalculator::RangeGridCalculator(BC bc) : bc_(bc) {}
 
@@ -40,25 +44,29 @@ auto RangeGridCalculator::operator()(XsGridData const& orig_data,
     XsGridData data;
 
     auto calc_dedx = [&] {
-        if (!orig_data.derivative.empty() || bc_ == BC::size_
-            || orig_data.value.size() < 5)
+        if (orig_data.value.size() < 5 || bc_ == BC::size_)
         {
-            // Either the derivatives have already been calculated or linear
-            // interpolation should be used
-            return EnergyLossCalculator(orig_data, orig_reals);
+            // Use linear interpolation
+            data = orig_data;
+            data.derivative = {};
+            return EnergyLossCalculator(data, orig_reals);
         }
+        else if (orig_data.derivative.empty())
+        {
+            // Calculate the second derivatives for cubic spline interpolation
+            auto deriv = SplineDerivCalculator(bc_)(orig_data, orig_reals);
 
-        // Calculate the second derivatives for cubic spline interpolation
-        auto deriv = SplineDerivCalculator(bc_)(orig_data, orig_reals);
-
-        // Create a copy of the grid data with the derivatives
-        CollectionBuilder build(&host_reals);
-        data.log_energy = orig_data.log_energy;
-        data.value = build.insert_back(orig_reals[orig_data.value].begin(),
-                                       orig_reals[orig_data.value].end());
-        data.derivative = build.insert_back(deriv.begin(), deriv.end());
-        reals = host_reals;
-        return EnergyLossCalculator(data, reals);
+            // Create a copy of the grid data with the derivatives
+            CollectionBuilder build(&host_reals);
+            data.log_energy = orig_data.log_energy;
+            data.value = build.insert_back(orig_reals[orig_data.value].begin(),
+                                           orig_reals[orig_data.value].end());
+            data.derivative = build.insert_back(deriv.begin(), deriv.end());
+            reals = host_reals;
+            return EnergyLossCalculator(data, reals);
+        }
+        // The derivatives have already been calculated
+        return EnergyLossCalculator(orig_data, orig_reals);
     }();
 
     UniformGrid loge_grid(orig_data.log_energy);
@@ -71,6 +79,7 @@ auto RangeGridCalculator::operator()(XsGridData const& orig_data,
     real_type cum_range = 2 * std::exp(loge_grid[0]) / calc_dedx[0];
     result[0] = cum_range;
 
+    // Integrate the range from the energy loss
     for (size_type i = 1; i < loge_grid.size(); ++i)
     {
         real_type energy_lower = std::exp(loge_grid[i - 1]);
@@ -81,14 +90,17 @@ auto RangeGridCalculator::operator()(XsGridData const& orig_data,
         {
             energy -= delta_energy;
             real_type dedx = calc_dedx(units::MevEnergy(energy));
+
+            // Spline interpolation can exhibit oscillations that greatly
+            // affect the accuracy when the number of grid points is small and
+            // the scale of the x grid is large
             CELER_VALIDATE(dedx > 0,
                            << "negative value in range calculation: the "
                               "interpolation method may be unstable");
             cum_range += delta_energy / dedx;
         }
         result[i] = cum_range;
     }
-
     return result;
 }
 

src/celeritas/grid/RangeGridCalculator.hh

@@ -35,6 +35,12 @@ namespace celeritas
  * possible with what we've been importing from Geant4, this performs the same
  * calculation as in Geant4's \c G4LossTableBuilder::BuildRangeTable, which
  * uses the midpoint rule with 100 substeps for improved accuracy.
+ *
+ * The calculator is constructed with the boundary conditions for cubic spline
+ * interpolation. If the default constructor is used, or if the number of grid
+ * points is less than 5, linear interpolation will be used instead.
+ *
+ * \todo support polynomial interpolation as well?
  */
 class RangeGridCalculator
 {
@@ -48,6 +54,9 @@ class RangeGridCalculator
     //!@}
 
   public:
+    // Default constructor
+    RangeGridCalculator();
+
     // Construct with boundary conditions for spline interpolation
     explicit RangeGridCalculator(BC);
 

src/celeritas/grid/ValueGridInserter.cc

@@ -31,16 +31,18 @@ auto ValueGridInserter::operator()(UniformGridData const& log_grid,
                                    size_type prime_index,
                                    SpanConstDbl values) -> XsIndex
 {
-    CELER_EXPECT(log_grid);
-    CELER_EXPECT(log_grid.size == values.size());
-    CELER_EXPECT(prime_index <= log_grid.size
-                 || prime_index == XsGridData::no_scaling());
+    return this->insert_xs(log_grid, prime_index, values);
+}
 
-    XsGridData grid;
-    grid.log_energy = log_grid;
-    grid.prime_index = prime_index;
-    grid.value = values_.insert_back(values.begin(), values.end());
-    return xs_grids_.push_back(grid);
+//---------------------------------------------------------------------------//
+/*!
+ * Add a grid of physics xs data.
+ */
+auto ValueGridInserter::operator()(UniformGridData const& log_grid,
+                                   size_type prime_index,
+                                   SpanConstFlt values) -> XsIndex
+{
+    return this->insert_xs(log_grid, prime_index, values);
 }
 
 //---------------------------------------------------------------------------//
@@ -53,5 +55,36 @@ auto ValueGridInserter::operator()(UniformGridData const& log_grid,
     return (*this)(log_grid, XsGridData::no_scaling(), values);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Add a grid of log-spaced data without 1/E scaling.
+ */
+auto ValueGridInserter::operator()(UniformGridData const& log_grid,
+                                   SpanConstFlt values) -> XsIndex
+{
+    return (*this)(log_grid, XsGridData::no_scaling(), values);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Add a grid of physics xs data.
+ */
+template<class T>
+auto ValueGridInserter::insert_xs(UniformGridData const& log_grid,
+                                  size_type prime_index,
+                                  Span<T const> values) -> XsIndex
+{
+    CELER_EXPECT(log_grid);
+    CELER_EXPECT(log_grid.size == values.size());
+    CELER_EXPECT(prime_index <= log_grid.size
+                 || prime_index == XsGridData::no_scaling());
+
+    XsGridData grid;
+    grid.log_energy = log_grid;
+    grid.prime_index = prime_index;
+    grid.value = values_.insert_back(values.begin(), values.end());
+    return xs_grids_.push_back(grid);
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/grid/ValueGridInserter.hh

@@ -41,11 +41,12 @@ class ValueGridInserter
   public:
     //!@{
     //! \name Type aliases
+    using SpanConstFlt = Span<float const>;
+    using SpanConstDbl = Span<double const>;
     using RealCollection
         = Collection<real_type, Ownership::value, MemSpace::host>;
     using XsGridCollection
         = Collection<XsGridData, Ownership::value, MemSpace::host>;
-    using SpanConstDbl = Span<double const>;
     using XsIndex = ItemId<XsGridData>;
     //!@}
 
@@ -57,13 +58,20 @@ class ValueGridInserter
     XsIndex operator()(UniformGridData const& log_grid,
                        size_type prime_index,
                        SpanConstDbl values);
+    XsIndex operator()(UniformGridData const& log_grid,
+                       size_type prime_index,
+                       SpanConstFlt values);
 
     // Add a grid of uniform log-grid data
     XsIndex operator()(UniformGridData const& log_grid, SpanConstDbl values);
+    XsIndex operator()(UniformGridData const& log_grid, SpanConstFlt values);
 
   private:
     CollectionBuilder<real_type, MemSpace::host, ItemId<real_type>> values_;
     CollectionBuilder<XsGridData, MemSpace::host, ItemId<XsGridData>> xs_grids_;
+
+    template<class T>
+    XsIndex insert_xs(UniformGridData const&, size_type, Span<T const>);
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/phys/PhysicsParams.cc

@@ -564,17 +564,13 @@ void PhysicsParams::build_xs(Options const& opts,
                 if (auto grid_id = eloss_grid_ids[mat_idx])
                 {
                     //! \todo make the interpolation method configurable?
+
+                    // Build the range grid from the energy loss
                     auto const& grid = data->value_grids[grid_id];
-                    auto range = RangeGridCalculator(BC::geant)(
+                    auto const range = RangeGridCalculator(BC::geant)(
                         grid, make_const_ref(*data).reals);
-
-                    XsGridData grid_data;
-                    grid_data.log_energy = grid.log_energy;
-                    grid_data.value
-                        = make_builder(&data->reals)
-                              .insert_back(range.begin(), range.end());
                     range_grid_ids[mat_idx]
-                        = make_builder(&data->value_grids).push_back(grid_data);
+                        = insert_grid(grid.log_energy, make_span(range));
                 }
 
                 if (use_integral_xs)

test/celeritas/grid/RangeGridCalculator.test.cc

@@ -166,7 +166,7 @@ TEST_F(RangeGridCalculatorTest, sparse)
     }
     {
         // Linear interpolation
-        auto range = RangeGridCalculator(BC::size_)(dedx_grid, reals_ref);
+        auto range = RangeGridCalculator()(dedx_grid, reals_ref);
         static double const expected_range[] = {
             2.3818927937551e-07,
             1.7075905920016e-06,

test/celeritas/grid/ValueGridInserter.test.cc

@@ -7,6 +7,7 @@
 #include "celeritas/grid/ValueGridInserter.hh"
 
 #include <algorithm>
+#include <vector>
 
 #include "corecel/cont/Range.hh"
 
@@ -36,7 +37,7 @@ TEST_F(ValueGridInserterTest, all)
     ValueGridInserter insert(&real_storage, &grid_storage);
 
     {
-        double const values[] = {10, 20, 3};
+        std::vector<double> const values = {10, 20, 3};
 
         auto idx = insert(
             UniformGridData::from_bounds(0.0, 1.0, 3), 1, make_span(values));
@@ -48,7 +49,7 @@ TEST_F(ValueGridInserterTest, all)
         EXPECT_VEC_SOFT_EQ(values, real_storage[inserted.value]);
     }
     {
-        double const values[] = {1, 2, 4, 6, 8};
+        std::vector<double> const values = {1, 2, 4, 6, 8};
 
         auto idx = insert(UniformGridData::from_bounds(0.0, 10.0, 5),
                           make_span(values));