From a7b9386830ebc5f58274d5860e81cb185f339159 Mon Sep 17 00:00:00 2001 From: Amanda Lund Date: Sat, 8 Feb 2025 20:26:01 +0000 Subject: [PATCH] Use ValueGridInserter and support linear interpolation --- src/celeritas/grid/RangeGridCalculator.cc | 50 +++++++++++------- src/celeritas/grid/RangeGridCalculator.hh | 9 ++++ src/celeritas/grid/ValueGridInserter.cc | 51 +++++++++++++++---- src/celeritas/grid/ValueGridInserter.hh | 10 +++- src/celeritas/phys/PhysicsParams.cc | 12 ++--- .../grid/RangeGridCalculator.test.cc | 2 +- test/celeritas/grid/ValueGridInserter.test.cc | 5 +- 7 files changed, 99 insertions(+), 40 deletions(-) diff --git a/src/celeritas/grid/RangeGridCalculator.cc b/src/celeritas/grid/RangeGridCalculator.cc index a1ead090c7..df3970ec9b 100644 --- a/src/celeritas/grid/RangeGridCalculator.cc +++ b/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,6 +90,10 @@ 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"); @@ -88,7 +101,6 @@ auto RangeGridCalculator::operator()(XsGridData const& orig_data, } result[i] = cum_range; } - return result; } diff --git a/src/celeritas/grid/RangeGridCalculator.hh b/src/celeritas/grid/RangeGridCalculator.hh index 783a560e0c..2dc3e15496 100644 --- a/src/celeritas/grid/RangeGridCalculator.hh +++ b/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); diff --git a/src/celeritas/grid/ValueGridInserter.cc b/src/celeritas/grid/ValueGridInserter.cc index 8c6fb75dff..402732d8b3 100644 --- a/src/celeritas/grid/ValueGridInserter.cc +++ b/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 +auto ValueGridInserter::insert_xs(UniformGridData const& log_grid, + size_type prime_index, + Span 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 diff --git a/src/celeritas/grid/ValueGridInserter.hh b/src/celeritas/grid/ValueGridInserter.hh index 82aa5c0b1a..ac910f77c0 100644 --- a/src/celeritas/grid/ValueGridInserter.hh +++ b/src/celeritas/grid/ValueGridInserter.hh @@ -41,11 +41,12 @@ class ValueGridInserter public: //!@{ //! \name Type aliases + using SpanConstFlt = Span; + using SpanConstDbl = Span; using RealCollection = Collection; using XsGridCollection = Collection; - using SpanConstDbl = Span; using XsIndex = ItemId; //!@} @@ -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> values_; CollectionBuilder> xs_grids_; + + template + XsIndex insert_xs(UniformGridData const&, size_type, Span); }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/phys/PhysicsParams.cc b/src/celeritas/phys/PhysicsParams.cc index b760dea83e..712606f5e0 100644 --- a/src/celeritas/phys/PhysicsParams.cc +++ b/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) diff --git a/test/celeritas/grid/RangeGridCalculator.test.cc b/test/celeritas/grid/RangeGridCalculator.test.cc index 021f67e117..f9a0182906 100644 --- a/test/celeritas/grid/RangeGridCalculator.test.cc +++ b/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, diff --git a/test/celeritas/grid/ValueGridInserter.test.cc b/test/celeritas/grid/ValueGridInserter.test.cc index 78e892f67a..ac120e1061 100644 --- a/test/celeritas/grid/ValueGridInserter.test.cc +++ b/test/celeritas/grid/ValueGridInserter.test.cc @@ -7,6 +7,7 @@ #include "celeritas/grid/ValueGridInserter.hh" #include +#include #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 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 const values = {1, 2, 4, 6, 8}; auto idx = insert(UniformGridData::from_bounds(0.0, 10.0, 5), make_span(values));