From 1bfafd34b13fccd5f10864eee966bd95051f28bd Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 11:47:27 -0500 Subject: [PATCH 01/10] Rename PhysicsGridCalculator to XsCalculator --- app/demo-interactor/HostKNDemoRunner.cc | 4 +- ...ysicsGridCalculator.hh => XsCalculator.hh} | 18 ++--- ...sGridCalculator.i.hh => XsCalculator.i.hh} | 36 ++++++---- test/CMakeLists.txt | 3 +- test/physics/em/LivermorePE.test.cc | 4 +- test/physics/grid/CalculatorTestBase.cc | 69 ++++++++++++++++++ test/physics/grid/CalculatorTestBase.hh | 56 +++++++++++++++ test/physics/grid/ValueGridBuilder.test.cc | 12 ++-- ...alculator.test.cc => XsCalculator.test.cc} | 71 ++++++------------- 9 files changed, 189 insertions(+), 84 deletions(-) rename src/physics/grid/{PhysicsGridCalculator.hh => XsCalculator.hh} (81%) rename src/physics/grid/{PhysicsGridCalculator.i.hh => XsCalculator.i.hh} (67%) create mode 100644 test/physics/grid/CalculatorTestBase.cc create mode 100644 test/physics/grid/CalculatorTestBase.hh rename test/physics/grid/{PhysicsGridCalculator.test.cc => XsCalculator.test.cc} (62%) diff --git a/app/demo-interactor/HostKNDemoRunner.cc b/app/demo-interactor/HostKNDemoRunner.cc index da960ccc06..85bfd2b9c4 100644 --- a/app/demo-interactor/HostKNDemoRunner.cc +++ b/app/demo-interactor/HostKNDemoRunner.cc @@ -19,7 +19,7 @@ #include "physics/base/Secondary.hh" #include "physics/base/SecondaryAllocatorView.hh" #include "physics/em/detail/KleinNishinaInteractor.hh" -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "DetectorView.hh" #include "HostStackAllocatorStore.hh" #include "HostDetectorStore.hh" @@ -75,7 +75,7 @@ auto HostKNDemoRunner::operator()(demo_interactor::KNDemoRunArgs args) // Physics calculator const auto& xs_host_ptrs = xsparams_->host_pointers(); - PhysicsGridCalculator calc_xs(xs_host_ptrs.xs, xs_host_ptrs.reals); + XsCalculator calc_xs(xs_host_ptrs.xs, xs_host_ptrs.reals); // Make secondary store HostStackAllocatorStore secondaries(args.max_steps); diff --git a/src/physics/grid/PhysicsGridCalculator.hh b/src/physics/grid/XsCalculator.hh similarity index 81% rename from src/physics/grid/PhysicsGridCalculator.hh rename to src/physics/grid/XsCalculator.hh index b50155e1fc..f625900e8f 100644 --- a/src/physics/grid/PhysicsGridCalculator.hh +++ b/src/physics/grid/XsCalculator.hh @@ -3,7 +3,7 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file PhysicsGridCalculator.hh +//! \file XsCalculator.hh //---------------------------------------------------------------------------// #pragma once @@ -14,7 +14,7 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Find and interpolate physics data based on a track's energy. + * Find and interpolate cross sections on a uniform log grid. * * \todo Currently this is hard-coded to use "cross section grid pointers" * which have energy coordinates uniform in log space. This should @@ -26,11 +26,11 @@ namespace celeritas * scaled by the energy. * * \code - PhysicsGridCalculator calc_xs(xs_grid, xs_params.reals); + XsCalculator calc_xs(xs_grid, xs_params.reals); real_type xs = calc_xs(particle); \endcode */ -class PhysicsGridCalculator +class XsCalculator { public: //!@{ @@ -42,17 +42,19 @@ class PhysicsGridCalculator public: // Construct from state-independent data inline CELER_FUNCTION - PhysicsGridCalculator(const XsGridData& grid, const Values& values); + XsCalculator(const XsGridData& grid, const Values& values); // Find and interpolate from the energy inline CELER_FUNCTION real_type operator()(Energy energy) const; private: - const XsGridData& data_; - Span values_; + const XsGridData& data_; + const Values& reals_; + + CELER_FORCEINLINE_FUNCTION real_type get(size_type index) const; }; //---------------------------------------------------------------------------// } // namespace celeritas -#include "PhysicsGridCalculator.i.hh" +#include "XsCalculator.i.hh" diff --git a/src/physics/grid/PhysicsGridCalculator.i.hh b/src/physics/grid/XsCalculator.i.hh similarity index 67% rename from src/physics/grid/PhysicsGridCalculator.i.hh rename to src/physics/grid/XsCalculator.i.hh index d9707bc873..cba0463ad0 100644 --- a/src/physics/grid/PhysicsGridCalculator.i.hh +++ b/src/physics/grid/XsCalculator.i.hh @@ -3,9 +3,8 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file PhysicsGridCalculator.i.hh +//! \file XsCalculator.i.hh //---------------------------------------------------------------------------// - #include #include "base/Interpolator.hh" #include "physics/grid/UniformGrid.hh" @@ -17,12 +16,11 @@ namespace celeritas * Construct from cross section data. */ CELER_FUNCTION -PhysicsGridCalculator::PhysicsGridCalculator(const XsGridData& grid, - const Values& values) - : data_(grid), values_(values[grid.value]) +XsCalculator::XsCalculator(const XsGridData& grid, const Values& values) + : data_(grid), reals_(values) { CELER_EXPECT(data_); - CELER_ASSERT(values_.size() == data_.log_energy.size); + CELER_ASSERT(grid.value.size() == data_.log_energy.size); } //---------------------------------------------------------------------------// @@ -32,10 +30,10 @@ PhysicsGridCalculator::PhysicsGridCalculator(const XsGridData& grid, * If needed, we can add a "log(energy/MeV)" accessor if we constantly reuse * that value and don't want to repeat the `std::log` operation. */ -CELER_FUNCTION real_type PhysicsGridCalculator::operator()(Energy energy) const +CELER_FUNCTION real_type XsCalculator::operator()(Energy energy) const { - UniformGrid loge_grid(data_.log_energy); - real_type loge = std::log(energy.value()); + const UniformGrid loge_grid(data_.log_energy); + const real_type loge = std::log(energy.value()); // Snap out-of-bounds values to closest grid points size_type lower_idx; @@ -43,12 +41,12 @@ CELER_FUNCTION real_type PhysicsGridCalculator::operator()(Energy energy) const if (loge <= loge_grid.front()) { lower_idx = 0; - result = values_[lower_idx]; + result = this->get(lower_idx); } else if (loge >= loge_grid.back()) { lower_idx = loge_grid.size() - 1; - result = values_[lower_idx]; + result = this->get(lower_idx); } else { @@ -56,8 +54,8 @@ CELER_FUNCTION real_type PhysicsGridCalculator::operator()(Energy energy) const lower_idx = loge_grid.find(loge); CELER_ASSERT(lower_idx + 1 < loge_grid.size()); - real_type upper_xs = values_[lower_idx + 1]; - real_type upper_energy = std::exp(loge_grid[lower_idx + 1]); + const real_type 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) { // Cross section data for the upper point has *already* been scaled @@ -67,7 +65,7 @@ CELER_FUNCTION real_type PhysicsGridCalculator::operator()(Energy energy) const // Interpolate *linearly* on energy using the lower_idx data. LinearInterpolator interpolate_xs( - {std::exp(loge_grid[lower_idx]), values_[lower_idx]}, + {std::exp(loge_grid[lower_idx]), this->get(lower_idx)}, {upper_energy, upper_xs}); result = interpolate_xs(energy.value()); } @@ -79,5 +77,15 @@ CELER_FUNCTION real_type PhysicsGridCalculator::operator()(Energy energy) const return result; } +//---------------------------------------------------------------------------// +/*! + * Get the raw cross section data at a particular index. + */ +CELER_FUNCTION real_type XsCalculator::get(size_type index) const +{ + CELER_EXPECT(index < data_.value.size()); + return reals_[data_.value[index]]; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fbfd077ade..cdda2a520b 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -167,6 +167,7 @@ add_library(CeleritasPhysicsTest physics/base/MockModel.cc physics/base/MockProcess.cc physics/base/PhysicsTestBase.cc + physics/grid/CalculatorTestBase.cc ) target_link_libraries(CeleritasPhysicsTest PRIVATE celeritas CeleritasTest) set(CELERITASTEST_LINK_LIBRARIES CeleritasPhysicsTest) @@ -178,10 +179,10 @@ celeritas_add_test(physics/base/PhysicsStepUtils.test.cc) celeritas_setup_tests(SERIAL PREFIX physics/grid) celeritas_add_test(physics/grid/GridIdFinder.test.cc) -celeritas_add_test(physics/grid/PhysicsGridCalculator.test.cc) celeritas_add_test(physics/grid/UniformGrid.test.cc) celeritas_add_test(physics/grid/ValueGridBuilder.test.cc) celeritas_add_test(physics/grid/ValueGridInserter.test.cc) +celeritas_add_test(physics/grid/XsCalculator.test.cc) celeritas_setup_tests(SERIAL PREFIX physics/material) celeritas_add_test(physics/material/ElementSelector.test.cc) diff --git a/test/physics/em/LivermorePE.test.cc b/test/physics/em/LivermorePE.test.cc index e8f42a5d24..3720eeadbc 100644 --- a/test/physics/em/LivermorePE.test.cc +++ b/test/physics/em/LivermorePE.test.cc @@ -23,7 +23,7 @@ #include "physics/em/LivermorePEParams.hh" #include "physics/em/PhotoelectricProcess.hh" #include "physics/em/LivermorePEMacroXsCalculator.hh" -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "physics/grid/ValueGridBuilder.hh" #include "physics/grid/ValueGridInserter.hh" #include "physics/material/MaterialTrackView.hh" @@ -504,7 +504,7 @@ TEST_F(LivermorePEInteractorTest, model) EXPECT_EQ(1, grid_storage.size()); // Test cross sections calculated from tables - celeritas::PhysicsGridCalculator calc_xs( + celeritas::XsCalculator calc_xs( grid_storage[ValueGridInserter::XsIndex{0}], real_storage); EXPECT_SOFT_EQ(0.1, calc_xs(MevEnergy{1e-3})); EXPECT_SOFT_EQ(1e-5, calc_xs(MevEnergy{1e2})); diff --git a/test/physics/grid/CalculatorTestBase.cc b/test/physics/grid/CalculatorTestBase.cc new file mode 100644 index 0000000000..147ebfb5bb --- /dev/null +++ b/test/physics/grid/CalculatorTestBase.cc @@ -0,0 +1,69 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file CalculatorTestBase.cc +//---------------------------------------------------------------------------// +#include "CalculatorTestBase.hh" + +#include +#include +#include "base/Interpolator.hh" +#include "base/PieBuilder.hh" +#include "base/Range.hh" +#include "base/SoftEqual.hh" + +using namespace celeritas; + +namespace celeritas_test +{ +//---------------------------------------------------------------------------// +void CalculatorTestBase::build(real_type emin, real_type emax, size_type count) +{ + CELER_EXPECT(count >= 2); + data_.log_energy + = UniformGridData::from_bounds(std::log(emin), std::log(emax), count); + + std::vector temp_xs(count); + + // Interpolate xs grid: linear in bin, log in energy + Interpolator calc_xs( + {0.0, emin}, {count - 1.0, emax}); + for (auto i : range(temp_xs.size())) + { + temp_xs[i] = calc_xs(i); + } + + value_storage_ = {}; + data_.value = make_pie_builder(&value_storage_) + .insert_back(temp_xs.begin(), temp_xs.end()); + value_ref_ = value_storage_; + + CELER_ENSURE(data_); + CELER_ENSURE(soft_equal(emax, value_ref_[data_.value].back())); +} + +//---------------------------------------------------------------------------// +/*! + * Set the index above which data is 1/E + */ +void CalculatorTestBase::set_prime_index(size_type i) +{ + CELER_EXPECT(data_); + CELER_EXPECT(i < data_.log_energy.size); + data_.prime_index = i; +} + +//---------------------------------------------------------------------------// +/*! + * Get cross sections that can be modified. + */ +auto CalculatorTestBase::mutable_values() -> SpanReal +{ + CELER_EXPECT(data_); + return value_storage_[data_.value]; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/physics/grid/CalculatorTestBase.hh b/test/physics/grid/CalculatorTestBase.hh new file mode 100644 index 0000000000..cac691c01d --- /dev/null +++ b/test/physics/grid/CalculatorTestBase.hh @@ -0,0 +1,56 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file CalculatorTestBase.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "gtest/Test.hh" +#include "base/Pie.hh" +#include "physics/grid/XsGridInterface.hh" + +namespace celeritas_test +{ +//---------------------------------------------------------------------------// +/*! + * Brief class description. + * + * Optional detailed class description, and possibly example usage: + * \code + CalculatorTestBase ...; + \endcode + */ +class CalculatorTestBase : public celeritas::Test +{ + public: + //!@{ + //! Type aliases + using real_type = celeritas::real_type; + using size_type = celeritas::size_type; + using XsGridData = celeritas::XsGridData; + using Pointers = celeritas::Pie; + using SpanReal = celeritas::Span; + //!@} + + public: + // Construct linear cross sections + void build(real_type emin, real_type emax, size_type count); + void set_prime_index(size_type i); + SpanReal mutable_values(); + + const XsGridData& data() const { return data_; } + const Pointers& values() const { return value_ref_; } + + private: + XsGridData data_; + celeritas::Pie + value_storage_; + Pointers value_ref_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/physics/grid/ValueGridBuilder.test.cc b/test/physics/grid/ValueGridBuilder.test.cc index 2e419f6e57..d2501fe515 100644 --- a/test/physics/grid/ValueGridBuilder.test.cc +++ b/test/physics/grid/ValueGridBuilder.test.cc @@ -9,7 +9,7 @@ #include #include -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "physics/grid/ValueGridInserter.hh" #include "celeritas_test.hh" @@ -26,7 +26,7 @@ class ValueGridBuilderTest : public celeritas::Test using SPConstBuilder = std::shared_ptr; using VecBuilder = std::vector; using VecReal = std::vector; - using Energy = PhysicsGridCalculator ::Energy; + using Energy = XsCalculator::Energy; using XsIndex = ValueGridInserter::XsIndex; protected: @@ -75,13 +75,13 @@ TEST_F(ValueGridBuilderTest, xs_grid) // Test results using the physics calculator ASSERT_EQ(2, grid_storage.size()); { - PhysicsGridCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); } { - PhysicsGridCalculator calc_xs(grid_storage[XsIndex{1}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{1}], real_storage); EXPECT_SOFT_EQ(10., calc_xs(Energy{1e-3})); EXPECT_SOFT_EQ(1., calc_xs(Energy{1e-2})); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e-1})); @@ -106,7 +106,7 @@ TEST_F(ValueGridBuilderTest, log_grid) // Test results using the physics calculator ASSERT_EQ(1, grid_storage.size()); { - PhysicsGridCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); @@ -133,7 +133,7 @@ TEST_F(ValueGridBuilderTest, DISABLED_generic_grid) // Test results using the physics calculator ASSERT_EQ(2, grid_storage.size()); { - PhysicsGridCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); diff --git a/test/physics/grid/PhysicsGridCalculator.test.cc b/test/physics/grid/XsCalculator.test.cc similarity index 62% rename from test/physics/grid/PhysicsGridCalculator.test.cc rename to test/physics/grid/XsCalculator.test.cc index 70d0e9ef63..6fa52ce051 100644 --- a/test/physics/grid/PhysicsGridCalculator.test.cc +++ b/test/physics/grid/XsCalculator.test.cc @@ -3,16 +3,13 @@ // See the top-level COPYRIGHT file for details. // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// -//! \file PhysicsGridCalculator.test.cc +//! \file XsCalculator.test.cc //---------------------------------------------------------------------------// -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include -#include -#include "base/Interpolator.hh" -#include "base/PieBuilder.hh" -#include "base/Range.hh" #include "celeritas_test.hh" +#include "CalculatorTestBase.hh" using namespace celeritas; @@ -20,51 +17,23 @@ using namespace celeritas; // TEST HARNESS //---------------------------------------------------------------------------// -class PhysicsGridCalculatorTest : public celeritas::Test +class XsCalculatorTest : public celeritas_test::CalculatorTestBase { protected: - using Energy = PhysicsGridCalculator::Energy; - using real_type = celeritas::real_type; - - void build(real_type emin, real_type emax, int count) - { - CELER_EXPECT(count >= 2); - data.log_energy = UniformGridData::from_bounds( - std::log(emin), std::log(emax), count); - - std::vector temp_xs(count); - - // Interpolate xs grid: linear in bin, log in energy - using celeritas::Interp; - celeritas::Interpolator calc_xs( - {0.0, emin}, {count - 1.0, emax}); - for (auto i : celeritas::range(temp_xs.size())) - { - temp_xs[i] = calc_xs(i); - } - - data.value = make_pie_builder(&stored_xs) - .insert_back(temp_xs.begin(), temp_xs.end()); - - CELER_ENSURE(data); - CELER_ENSURE(celeritas::soft_equal(emax, stored_xs[data.value].back())); - } - - XsGridData data; - Pie stored_xs; + using Energy = XsCalculator::Energy; }; //---------------------------------------------------------------------------// // TESTS //---------------------------------------------------------------------------// -TEST_F(PhysicsGridCalculatorTest, simple) +TEST_F(XsCalculatorTest, simple) { // Energy from 1 to 1e5 MeV with 5 grid points; XS should be the same // *No* magical 1/E scaling this->build(1.0, 1e5, 6); - PhysicsGridCalculator calc(this->data, this->stored_xs); + XsCalculator calc(this->data(), this->values()); // Test on grid points EXPECT_SOFT_EQ(1.0, calc(Energy{1})); @@ -80,14 +49,14 @@ TEST_F(PhysicsGridCalculatorTest, simple) EXPECT_SOFT_EQ(1e5, calc(Energy{1e7})); } -TEST_F(PhysicsGridCalculatorTest, scaled_lowest) +TEST_F(XsCalculatorTest, scaled_lowest) { // Energy from .1 to 1e4 MeV with 5 grid points; XS should be constant // since the constructor fills it with E this->build(0.1, 1e4, 6); - data.prime_index = 0; + this->set_prime_index(0); - PhysicsGridCalculator calc(this->data, this->stored_xs); + XsCalculator calc(this->data(), this->values()); // Test on grid points EXPECT_SOFT_EQ(1, calc(Energy{0.1})); @@ -105,13 +74,13 @@ TEST_F(PhysicsGridCalculatorTest, scaled_lowest) EXPECT_SOFT_EQ(0.1, calc(Energy{1e5})); } -TEST_F(PhysicsGridCalculatorTest, scaled_middle) +TEST_F(XsCalculatorTest, scaled_middle) { // Energy from .1 to 1e4 MeV with 5 grid points; XS should be constant // since the constructor fills it with E this->build(0.1, 1e4, 6); - data.prime_index = 3; - auto xs = this->stored_xs[data.value]; + this->set_prime_index(3); + auto xs = this->mutable_values(); std::fill(xs.begin(), xs.begin() + 3, 1.0); // Change constant to 3 just to shake things up @@ -120,7 +89,7 @@ TEST_F(PhysicsGridCalculatorTest, scaled_middle) xs *= 3; } - PhysicsGridCalculator calc(this->data, this->stored_xs); + XsCalculator calc(this->data(), this->values()); // Test on grid points EXPECT_SOFT_EQ(3, calc(Energy{0.1})); @@ -138,13 +107,13 @@ TEST_F(PhysicsGridCalculatorTest, scaled_middle) EXPECT_SOFT_EQ(0.3, calc(Energy{1e5})); } -TEST_F(PhysicsGridCalculatorTest, scaled_highest) +TEST_F(XsCalculatorTest, scaled_highest) { // values of 1, 10, 100 --> actual xs = {1, 10, 1} this->build(1, 100, 3); - data.prime_index = 2; + this->set_prime_index(2); - PhysicsGridCalculator calc(this->data, this->stored_xs); + XsCalculator calc(this->data(), this->values()); EXPECT_SOFT_EQ(1, calc(Energy{0.0001})); EXPECT_SOFT_EQ(1, calc(Energy{1})); EXPECT_SOFT_EQ(10, calc(Energy{10})); @@ -155,12 +124,12 @@ TEST_F(PhysicsGridCalculatorTest, scaled_highest) EXPECT_SOFT_EQ(.1, calc(Energy{1000})); } -TEST_F(PhysicsGridCalculatorTest, TEST_IF_CELERITAS_DEBUG(scaled_off_the_end)) +TEST_F(XsCalculatorTest, TEST_IF_CELERITAS_DEBUG(scaled_off_the_end)) { // values of 1, 10, 100 --> actual xs = {1, 10, 100} this->build(1, 100, 3); + XsGridData data(this->data()); data.prime_index = 3; // disallowed - EXPECT_THROW(PhysicsGridCalculator(this->data, this->stored_xs), - celeritas::DebugError); + EXPECT_THROW(XsCalculator(data, this->values()), celeritas::DebugError); } From be263a3f57a691221ce325fbb21288d048b2565f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 11:48:12 -0500 Subject: [PATCH 02/10] Add nonuniform grid --- src/physics/grid/NonuniformGrid.hh | 66 +++++++++++++++++++++++ src/physics/grid/NonuniformGrid.i.hh | 67 ++++++++++++++++++++++++ src/physics/grid/UniformGrid.hh | 10 ++-- src/physics/grid/UniformGrid.i.hh | 21 +------- test/CMakeLists.txt | 1 + test/physics/grid/NonuniformGrid.test.cc | 67 ++++++++++++++++++++++++ 6 files changed, 207 insertions(+), 25 deletions(-) create mode 100644 src/physics/grid/NonuniformGrid.hh create mode 100644 src/physics/grid/NonuniformGrid.i.hh create mode 100644 test/physics/grid/NonuniformGrid.test.cc diff --git a/src/physics/grid/NonuniformGrid.hh b/src/physics/grid/NonuniformGrid.hh new file mode 100644 index 0000000000..c81494c75f --- /dev/null +++ b/src/physics/grid/NonuniformGrid.hh @@ -0,0 +1,66 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file NonuniformGrid.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Macros.hh" +#include "base/Pie.hh" +#include "base/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Interact with a nonuniform grid of increasing values. + * + * This should have the same interface (aside from constructor) as + * UniformGrid. + */ +template +class NonuniformGrid +{ + public: + //!@{ + //! Type aliases + using size_type = ::celeritas::size_type; + using value_type = T; + using Values + = Pie; + //!@} + + public: + // Construct with data + explicit inline CELER_FUNCTION + NonuniformGrid(const PieSlice& values, const Values& data); + + //! Number of grid points + CELER_FORCEINLINE_FUNCTION size_type size() const { return data_.size(); } + + //! Minimum/first value + CELER_FORCEINLINE_FUNCTION value_type front() const + { + return data_.front(); + } + + //! Maximum/last value + CELER_FORCEINLINE_FUNCTION value_type back() const { return data_.back(); } + + // Calculate the value at the given grid point + inline CELER_FUNCTION value_type operator[](size_type i) const; + + // Find the index of the given value (*must* be in bounds) + inline CELER_FUNCTION size_type find(value_type value) const; + + private: + // TODO: change backend for effiency if needeed + Span data_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas + +#include "NonuniformGrid.i.hh" diff --git a/src/physics/grid/NonuniformGrid.i.hh b/src/physics/grid/NonuniformGrid.i.hh new file mode 100644 index 0000000000..c48314543b --- /dev/null +++ b/src/physics/grid/NonuniformGrid.i.hh @@ -0,0 +1,67 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file NonuniformGrid.i.hh +//---------------------------------------------------------------------------// +#include "base/Algorithms.hh" +#include "base/Assert.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct with data. + */ +template +CELER_FUNCTION +NonuniformGrid::NonuniformGrid(const PieSlice& values, + const Values& data) + : data_(data[values]) +{ + CELER_EXPECT(data_.size() >= 2); + CELER_EXPECT(data_.front() <= data_.back()); // Approximation for "sorted" +} + +//---------------------------------------------------------------------------// +/*! + * Get the value at the given grid point. + */ +template +CELER_FUNCTION auto NonuniformGrid::operator[](size_type i) const + -> value_type +{ + CELER_EXPECT(i < data_.size()); + return data_[i]; +} + +//---------------------------------------------------------------------------// +/*! + * Find the value bin such that data[result] <= value < data[result + 1]. + * + * The given value *must* be in range, because out-of-bounds values usually + * require different treatment (e.g. clipping to the boundary values rather + * than interpolating). It's easier to test the exceptional cases (final grid + * point) outside of the grid view. + */ +template +CELER_FUNCTION size_type NonuniformGrid::find(value_type value) const +{ + CELER_EXPECT(value >= this->front() && value < this->back()); + + auto iter = celeritas::lower_bound(data_.begin(), data_.end(), value); + CELER_ASSERT(iter != data_.end()); + + if (value != *iter) + { + // Exactly on end grid point, or not on a grid point at all: move to + // previous bin + --iter; + } + + return iter - data_.begin(); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/physics/grid/UniformGrid.hh b/src/physics/grid/UniformGrid.hh index d6f487af39..771dcd7ff6 100644 --- a/src/physics/grid/UniformGrid.hh +++ b/src/physics/grid/UniformGrid.hh @@ -33,14 +33,14 @@ class UniformGrid // Construct with data explicit inline CELER_FUNCTION UniformGrid(const UniformGridData& data); - // Number of grid points - inline CELER_FUNCTION size_type size() const; + //! Number of grid points + CELER_FORCEINLINE_FUNCTION size_type size() const { return data_.size; } //! Minimum/first value - CELER_FUNCTION value_type front() const { return data_.front; } + CELER_FORCEINLINE_FUNCTION value_type front() const { return data_.front; } - // Maximum/last value - inline CELER_FUNCTION value_type back() const; + //! Maximum/last value + CELER_FORCEINLINE_FUNCTION value_type back() const { return data_.back; } // Calculate the value at the given grid point inline CELER_FUNCTION value_type operator[](size_type i) const; diff --git a/src/physics/grid/UniformGrid.i.hh b/src/physics/grid/UniformGrid.i.hh index 2f25f25d8a..67ea3b5870 100644 --- a/src/physics/grid/UniformGrid.i.hh +++ b/src/physics/grid/UniformGrid.i.hh @@ -5,7 +5,6 @@ //---------------------------------------------------------------------------// //! \file UniformGrid.i.hh //---------------------------------------------------------------------------// -#include #include "base/Assert.hh" namespace celeritas @@ -22,25 +21,7 @@ UniformGrid::UniformGrid(const UniformGridData& data) : data_(data) //---------------------------------------------------------------------------// /*! - * Access the number of grid points. - */ -CELER_FUNCTION size_type UniformGrid::size() const -{ - return data_.size; -} - -//---------------------------------------------------------------------------// -/*! - * Access the highest value. - */ -CELER_FUNCTION auto UniformGrid::back() const -> value_type -{ - return data_.back; -} - -//---------------------------------------------------------------------------// -/*! - * Get the value value at the given grid point. + * Get the value at the given grid point. */ CELER_FUNCTION auto UniformGrid::operator[](size_type i) const -> value_type { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cdda2a520b..e6855ac1ad 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -179,6 +179,7 @@ celeritas_add_test(physics/base/PhysicsStepUtils.test.cc) celeritas_setup_tests(SERIAL PREFIX physics/grid) celeritas_add_test(physics/grid/GridIdFinder.test.cc) +celeritas_add_test(physics/grid/NonuniformGrid.test.cc) celeritas_add_test(physics/grid/UniformGrid.test.cc) celeritas_add_test(physics/grid/ValueGridBuilder.test.cc) celeritas_add_test(physics/grid/ValueGridInserter.test.cc) diff --git a/test/physics/grid/NonuniformGrid.test.cc b/test/physics/grid/NonuniformGrid.test.cc new file mode 100644 index 0000000000..2172b51af0 --- /dev/null +++ b/test/physics/grid/NonuniformGrid.test.cc @@ -0,0 +1,67 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file NonuniformGrid.test.cc +//---------------------------------------------------------------------------// +#include "physics/grid/NonuniformGrid.hh" + +#include "base/Pie.hh" +#include "base/PieBuilder.hh" +#include "celeritas_test.hh" + +using celeritas::NonuniformGrid; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class NonuniformGridTest : public celeritas::Test +{ + protected: + using GridT = NonuniformGrid; + + void SetUp() override + { + auto build = celeritas::make_pie_builder(&data); + slc = build.insert_back({0, 1, 3, 3, 7}); + } + + celeritas::PieSlice slc; + celeritas::Pie + data; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(NonuniformGridTest, accessors) +{ + GridT grid(slc, data); + + EXPECT_EQ(5, grid.size()); + EXPECT_EQ(0, grid.front()); + EXPECT_EQ(7, grid.back()); + EXPECT_EQ(0, grid[0]); + EXPECT_EQ(3, grid[2]); +} + +TEST_F(NonuniformGridTest, find) +{ + GridT grid(slc, data); + +#if CELERITAS_DEBUG + EXPECT_THROW(grid.find(-1), celeritas::DebugError); +#endif + EXPECT_EQ(0, grid.find(0)); + EXPECT_EQ(1, grid.find(1)); + EXPECT_EQ(1, grid.find(2)); + EXPECT_EQ(2, grid.find(3)); + EXPECT_EQ(3, grid.find(4)); +#if CELERITAS_DEBUG + EXPECT_THROW(grid.find(7), celeritas::DebugError); + EXPECT_THROW(grid.find(10), celeritas::DebugError); +#endif +} From 022a929a67da1f19e9eb27bed241e06f63d01891 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 12:03:46 -0500 Subject: [PATCH 03/10] Add range and inverse range calculators --- src/physics/grid/EnergyLossCalculator.hh | 19 +++++ src/physics/grid/InverseRangeCalculator.hh | 60 ++++++++++++++++ src/physics/grid/InverseRangeCalculator.i.hh | 67 +++++++++++++++++ src/physics/grid/RangeCalculator.hh | 56 +++++++++++++++ src/physics/grid/RangeCalculator.i.hh | 72 +++++++++++++++++++ test/CMakeLists.txt | 2 + .../grid/InverseRangeCalculator.test.cc | 68 ++++++++++++++++++ test/physics/grid/RangeCalculator.test.cc | 57 +++++++++++++++ 8 files changed, 401 insertions(+) create mode 100644 src/physics/grid/EnergyLossCalculator.hh create mode 100644 src/physics/grid/InverseRangeCalculator.hh create mode 100644 src/physics/grid/InverseRangeCalculator.i.hh create mode 100644 src/physics/grid/RangeCalculator.hh create mode 100644 src/physics/grid/RangeCalculator.i.hh create mode 100644 test/physics/grid/InverseRangeCalculator.test.cc create mode 100644 test/physics/grid/RangeCalculator.test.cc diff --git a/src/physics/grid/EnergyLossCalculator.hh b/src/physics/grid/EnergyLossCalculator.hh new file mode 100644 index 0000000000..07f9cff1ca --- /dev/null +++ b/src/physics/grid/EnergyLossCalculator.hh @@ -0,0 +1,19 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file EnergyLossCalculator.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "XsCalculator.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +//! For now, energy loss calculation has the same behavior as cross sections +using EnergyLossCalculator = XsCalculator; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/physics/grid/InverseRangeCalculator.hh b/src/physics/grid/InverseRangeCalculator.hh new file mode 100644 index 0000000000..f591ad2af4 --- /dev/null +++ b/src/physics/grid/InverseRangeCalculator.hh @@ -0,0 +1,60 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file InverseRangeCalculator.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Quantity.hh" +#include "NonuniformGrid.hh" +#include "XsGridInterface.hh" +#include "UniformGrid.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Calculate the energy that would limit a particle to a particular range. + * + * This should provide the inverse of the result of \c RangeCalculator. The + * given \c range is not allowed to be greater than the maximum range in the + * physics data. + * + * The range must be monotonically increasing in energy, since it's defined as + * the integral of the inverse of the stopping power (which is always + * positive). For ranges shorter than the minimum energy in the table, the + * resulting energy is scaled: + * \f[ + E = E_\textrm{min}} \left( \frac{r}{r_\textrm{min}} \right)^2 + * \f] + * This scaling is the inverse of the off-the-end energy scaling in the + * RangeCalculator. + */ +class InverseRangeCalculator +{ + public: + //!@{ + //! Type aliases + using Energy = Quantity; + using Values = Pie; + //!@} + + public: + // Construct from state-independent data + inline CELER_FUNCTION + InverseRangeCalculator(const XsGridData& grid, const Values& values); + + // Find and interpolate from the energy + inline CELER_FUNCTION Energy operator()(real_type range) const; + + private: + UniformGrid log_energy_; + NonuniformGrid range_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas + +#include "InverseRangeCalculator.i.hh" diff --git a/src/physics/grid/InverseRangeCalculator.i.hh b/src/physics/grid/InverseRangeCalculator.i.hh new file mode 100644 index 0000000000..ed21de634d --- /dev/null +++ b/src/physics/grid/InverseRangeCalculator.i.hh @@ -0,0 +1,67 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file InverseRangeCalculator.i.hh +//---------------------------------------------------------------------------// +#include +#include "base/Algorithms.hh" +#include "base/Assert.hh" +#include "base/Interpolator.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from range data. + * + * The range is expected to be monotonically increaing with energy. + * Lower-energy particles have shorter ranges. + */ +CELER_FUNCTION +InverseRangeCalculator::InverseRangeCalculator(const XsGridData& grid, + const Values& values) + : log_energy_(grid.log_energy), range_(grid.value, values) +{ + CELER_EXPECT(range_.size() == log_energy_.size()); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the energy of a particle that has the given range. + */ +CELER_FUNCTION auto InverseRangeCalculator::operator()(real_type range) const + -> Energy +{ + CELER_EXPECT(range >= 0 && range <= range_.back()); + + if (range < range_.front()) + { + // Very short range: this corresponds to "energy < emin" for range + // calculation: range = r[0] * sqrt(E / E[0]) + return Energy{std::exp(log_energy_.front()) + * ipow<2>(range / range_.front())}; + } + // Range should *never* exceed the longest range (highest energy) since + // that should have limited the step + if (CELER_UNLIKELY(range >= range_.back())) + { + CELER_ASSERT(range == range_.back()); + return Energy{std::exp(log_energy_.back())}; + } + + // Search for lower bin index + auto idx = range_.find(range); + CELER_ASSERT(idx + 1 < log_energy_.size()); + + // Interpolate: 'x' = range, y = log energy + LinearInterpolator interpolate_log_energy( + {range_[idx], std::exp(log_energy_[idx])}, + {range_[idx + 1], std::exp(log_energy_[idx + 1])}); + auto loge = interpolate_log_energy(range); + return Energy{loge}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/physics/grid/RangeCalculator.hh b/src/physics/grid/RangeCalculator.hh new file mode 100644 index 0000000000..3d866d5a63 --- /dev/null +++ b/src/physics/grid/RangeCalculator.hh @@ -0,0 +1,56 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file RangeCalculator.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "base/Quantity.hh" +#include "XsGridInterface.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Find and interpolate range on a uniform log grid. + * + * \code + RangeCalculator calc_range(xs_grid, xs_params.reals); + real_type range = calc_range(particle); + \endcode + * + * Below the minimum tabulated energy, the range is scaled: + * \f[ + r = r_\textrm{min}} \sqrt{ \frac{E}{E_\textrm{min}}} + * \f] + */ +class RangeCalculator +{ + public: + //!@{ + //! Type aliases + using Energy = Quantity; + using Values = Pie; + //!@} + + public: + // Construct from state-independent data + inline CELER_FUNCTION + RangeCalculator(const XsGridData& grid, const Values& values); + + // Find and interpolate from the energy + inline CELER_FUNCTION real_type operator()(Energy energy) const; + + private: + const XsGridData& data_; + const Values& reals_; + + CELER_FORCEINLINE_FUNCTION real_type get(size_type index) const; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas + +#include "RangeCalculator.i.hh" diff --git a/src/physics/grid/RangeCalculator.i.hh b/src/physics/grid/RangeCalculator.i.hh new file mode 100644 index 0000000000..ae680cba0a --- /dev/null +++ b/src/physics/grid/RangeCalculator.i.hh @@ -0,0 +1,72 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file RangeCalculator.i.hh +//---------------------------------------------------------------------------// +#include +#include "base/Interpolator.hh" +#include "physics/grid/UniformGrid.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from cross section data. + * + * Range tables should be uniform in energy, without extra scaling. + */ +CELER_FUNCTION +RangeCalculator::RangeCalculator(const XsGridData& grid, const Values& values) + : data_(grid), reals_(values) +{ + CELER_EXPECT(data_); + CELER_EXPECT(data_.prime_index == XsGridData::no_scaling()); +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the range. + */ +CELER_FUNCTION real_type RangeCalculator::operator()(Energy energy) const +{ + UniformGrid loge_grid(data_.log_energy); + const real_type loge = std::log(energy.value()); + + if (loge <= loge_grid.front()) + { + real_type result = this->get(0); + // Scale by sqrt(E/Emin) = exp(.5 (log E - log Emin)) + result *= std::exp(real_type(.5) * (loge - loge_grid.front())); + return result; + } + else if (loge >= loge_grid.back()) + { + // Clip to highest range value + return this->get(loge_grid.size() - 1); + } + + // Locate the energy bin + auto idx = loge_grid.find(loge); + CELER_ASSERT(idx + 1 < loge_grid.size()); + + // Interpolate *linearly* on energy + LinearInterpolator interpolate_xs( + {std::exp(loge_grid[idx]), this->get(idx)}, + {std::exp(loge_grid[idx + 1]), this->get(idx + 1)}); + return interpolate_xs(energy.value()); +} + +//---------------------------------------------------------------------------// +/*! + * Get the raw range data at a particular index. + */ +CELER_FUNCTION real_type RangeCalculator::get(size_type index) const +{ + CELER_EXPECT(index < reals_.size()); + return reals_[data_.value[index]]; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e6855ac1ad..f863c985e1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -179,7 +179,9 @@ celeritas_add_test(physics/base/PhysicsStepUtils.test.cc) celeritas_setup_tests(SERIAL PREFIX physics/grid) celeritas_add_test(physics/grid/GridIdFinder.test.cc) +celeritas_add_test(physics/grid/InverseRangeCalculator.test.cc) celeritas_add_test(physics/grid/NonuniformGrid.test.cc) +celeritas_add_test(physics/grid/RangeCalculator.test.cc) celeritas_add_test(physics/grid/UniformGrid.test.cc) celeritas_add_test(physics/grid/ValueGridBuilder.test.cc) celeritas_add_test(physics/grid/ValueGridInserter.test.cc) diff --git a/test/physics/grid/InverseRangeCalculator.test.cc b/test/physics/grid/InverseRangeCalculator.test.cc new file mode 100644 index 0000000000..15768d0313 --- /dev/null +++ b/test/physics/grid/InverseRangeCalculator.test.cc @@ -0,0 +1,68 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file InverseRangeCalculator.test.cc +//---------------------------------------------------------------------------// +#include "physics/grid/InverseRangeCalculator.hh" + +#include "base/SoftEqual.hh" +#include "CalculatorTestBase.hh" +#include "celeritas_test.hh" + +using celeritas::InverseRangeCalculator; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class InverseRangeCalculatorTest : public celeritas_test::CalculatorTestBase +{ + protected: + using Energy = InverseRangeCalculator::Energy; + + void SetUp() override + { + // Energy from 1e1 to 1e4 MeV with 3 bins (4 points) + this->build(10, 1e4, 4); + + // InverseRange is 1/20 of energy + auto value_span = this->mutable_values(); + for (real_type& xs : value_span) + { + xs *= .05; + } + + // Adjust final point for roundoff for exact top-of-range testing + CELER_ASSERT(celeritas::soft_equal(real_type(500), value_span.back())); + value_span.back() = 500; + } +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +// Note: these are all the same values as the RangeCalculator test. +TEST_F(InverseRangeCalculatorTest, all) +{ + InverseRangeCalculator calc_energy(this->data(), this->values()); + + // Values below should be scaled below emin + EXPECT_SOFT_EQ(1.0, calc_energy(.5 * std::sqrt(1. / 10.)).value()); + EXPECT_SOFT_EQ(2.0, calc_energy(.5 * std::sqrt(2. / 10.)).value()); + + // Values in range + EXPECT_SOFT_EQ(10.0, calc_energy(.5).value()); + EXPECT_SOFT_EQ(20.0, calc_energy(1).value()); + EXPECT_SOFT_EQ(100.0, calc_energy(5).value()); + + // Top of range + EXPECT_SOFT_EQ(1e4, calc_energy(500).value()); + +#if CELERITAS_DEBUG + // Above range + EXPECT_THROW(calc_energy(500.1), celeritas::DebugError); +#endif +} diff --git a/test/physics/grid/RangeCalculator.test.cc b/test/physics/grid/RangeCalculator.test.cc new file mode 100644 index 0000000000..4624add743 --- /dev/null +++ b/test/physics/grid/RangeCalculator.test.cc @@ -0,0 +1,57 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file RangeCalculator.test.cc +//---------------------------------------------------------------------------// +#include "physics/grid/RangeCalculator.hh" + +#include "celeritas_test.hh" +#include "CalculatorTestBase.hh" + +using celeritas::RangeCalculator; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class RangeCalculatorTest : public celeritas_test::CalculatorTestBase +{ + protected: + using Energy = RangeCalculator::Energy; + + void SetUp() override + { + // Energy from 1e1 to 1e4 MeV with 3 bins (4 points) + this->build(10, 1e4, 4); + + // Range is 1/20 of energy + for (real_type& xs : this->mutable_values()) + { + xs *= .05; + } + } +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(RangeCalculatorTest, all) +{ + RangeCalculator calc_range(this->data(), this->values()); + + // Values below should be scaled below emin + EXPECT_SOFT_EQ(.5 * std::sqrt(1. / 10.), calc_range(Energy{1})); + EXPECT_SOFT_EQ(.5 * std::sqrt(2. / 10.), calc_range(Energy{2})); + // Values in range + EXPECT_SOFT_EQ(0.5, calc_range(Energy{10})); + EXPECT_SOFT_EQ(1.0, calc_range(Energy{20})); + EXPECT_SOFT_EQ(5.0, calc_range(Energy{100})); + + // Top of range + EXPECT_SOFT_EQ(500, calc_range(Energy{1e4})); + // Above range + EXPECT_SOFT_EQ(500, calc_range(Energy{1.001e4})); +} From b1654c449f9e5c0c91963e1a887cfed95d7f8a77 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 21:13:49 -0500 Subject: [PATCH 04/10] WIP: step utils with wrong range Need to store *unmodified* range for each particle type for use in inverse range lookup. --- scripts/build/yuri.cmake | 8 +- scripts/build/yuri.sh | 7 +- src/physics/base/ParticleTrackView.hh | 5 +- src/physics/base/ParticleTrackView.i.hh | 4 +- src/physics/base/PhysicsInterface.hh | 15 +- src/physics/base/PhysicsParams.cc | 4 + src/physics/base/PhysicsParams.hh | 4 + src/physics/base/PhysicsStepUtils.hh | 16 +- src/physics/base/PhysicsStepUtils.i.hh | 128 +++++++++++---- src/physics/base/PhysicsTrackView.hh | 18 ++- src/physics/base/PhysicsTrackView.i.hh | 50 +++++- src/physics/grid/ValueGridBuilder.cc | 26 ++- src/physics/grid/ValueGridBuilder.hh | 28 +++- test/physics/base/MockProcess.cc | 27 ++-- test/physics/base/MockProcess.hh | 6 +- test/physics/base/Physics.test.cc | 22 +-- test/physics/base/Physics.test.hh | 6 +- test/physics/base/PhysicsStepUtils.test.cc | 178 +++++++++++++-------- test/physics/base/PhysicsTestBase.cc | 18 ++- test/physics/base/PhysicsTestBase.hh | 10 +- 20 files changed, 387 insertions(+), 193 deletions(-) diff --git a/scripts/build/yuri.cmake b/scripts/build/yuri.cmake index 2deef31d98..8172722043 100644 --- a/scripts/build/yuri.cmake +++ b/scripts/build/yuri.cmake @@ -2,15 +2,15 @@ macro(set_cache_var var type val) set(${var} "${val}" CACHE "${type}" "yuri.sh") endmacro() -set_cache_var(CELERITAS_BUILD_DOCS BOOL ON) +set_cache_var(CELERITAS_BUILD_DOCS BOOL OFF) # Dependency options set_cache_var(CELERITAS_USE_CUDA BOOL OFF) set_cache_var(CELERITAS_USE_Geant4 BOOL OFF) set_cache_var(CELERITAS_GIT_SUBMODULE BOOL OFF) -set_cache_var(CELERITAS_USE_MPI BOOL ON) -set_cache_var(CELERITAS_USE_ROOT BOOL ON) -set_cache_var(CELERITAS_USE_SWIG_Python BOOL ON) +set_cache_var(CELERITAS_USE_MPI BOOL OFF) +set_cache_var(CELERITAS_USE_ROOT BOOL OFF) +set_cache_var(CELERITAS_USE_SWIG_Python BOOL OFF) set_cache_var(CELERITAS_USE_VecGeom BOOL OFF) # Build options diff --git a/scripts/build/yuri.sh b/scripts/build/yuri.sh index 5750cd06fe..3774f0a7a0 100755 --- a/scripts/build/yuri.sh +++ b/scripts/build/yuri.sh @@ -4,8 +4,8 @@ SOURCE_DIR="$(cd "${BUILDSCRIPT_DIR}" && git rev-parse --show-toplevel)" printf "\e[2;37mBuilding from ${SOURCE_DIR}\e[0m\n" cd $SOURCE_DIR -mkdir build 2>/dev/null || true -cd build +mkdir build-xcode 2>/dev/null || true +cd build-xcode module load doxygen swig @@ -13,8 +13,7 @@ CELERITAS_ENV=${SPACK_ROOT}/var/spack/environments/celeritas/.spack-env/view export PATH=$CELERITAS_ENV/bin:${PATH} export CMAKE_PREFIX_PATH=$CELERITAS_ENV:${CMAKE_PREFIX_PATH} -cmake -C ${BUILDSCRIPT_DIR}/yuri.cmake -G Ninja \ - -DCMAKE_INSTALL_PREFIX:PATH=$SOURCE_DIR/install \ +cmake -C ${BUILDSCRIPT_DIR}/yuri.cmake -G Xcode \ .. ninja -v ctest --output-on-failure diff --git a/src/physics/base/ParticleTrackView.hh b/src/physics/base/ParticleTrackView.hh index 2d20595a0e..f1a515e36d 100644 --- a/src/physics/base/ParticleTrackView.hh +++ b/src/physics/base/ParticleTrackView.hh @@ -35,6 +35,7 @@ class ParticleTrackView = ParticleParamsData; using ParticleStateRef = ParticleStateData; + using Energy = units::MevEnergy; using Initializer_t = ParticleTrackState; //!@} @@ -49,7 +50,7 @@ class ParticleTrackView operator=(const Initializer_t& other); // Change the particle's energy [MeV] - inline CELER_FUNCTION void energy(units::MevEnergy); + inline CELER_FUNCTION void energy(Energy); //// DYNAMIC PROPERTIES (pure accessors, free) //// @@ -57,7 +58,7 @@ class ParticleTrackView CELER_FORCEINLINE_FUNCTION ParticleId particle_id() const; // Kinetic energy [MeV] - CELER_FORCEINLINE_FUNCTION units::MevEnergy energy() const; + CELER_FORCEINLINE_FUNCTION Energy energy() const; // Whether the particle is stopped (zero kinetic energy) CELER_FORCEINLINE_FUNCTION bool is_stopped() const; diff --git a/src/physics/base/ParticleTrackView.i.hh b/src/physics/base/ParticleTrackView.i.hh index c620004d54..821a4233a8 100644 --- a/src/physics/base/ParticleTrackView.i.hh +++ b/src/physics/base/ParticleTrackView.i.hh @@ -47,7 +47,7 @@ ParticleTrackView::operator=(const Initializer_t& other) * applications, the new energy should always be less than the starting energy. */ CELER_FUNCTION -void ParticleTrackView::energy(units::MevEnergy quantity) +void ParticleTrackView::energy(Energy quantity) { CELER_EXPECT(this->particle_id()); CELER_EXPECT(quantity >= zero_quantity()); @@ -69,7 +69,7 @@ CELER_FUNCTION ParticleId ParticleTrackView::particle_id() const /*! * Kinetic energy [MeV]. */ -CELER_FUNCTION units::MevEnergy ParticleTrackView::energy() const +CELER_FUNCTION auto ParticleTrackView::energy() const -> Energy { return states_.state[thread_].energy; } diff --git a/src/physics/base/PhysicsInterface.hh b/src/physics/base/PhysicsInterface.hh index 67fc834b90..cb2404df77 100644 --- a/src/physics/base/PhysicsInterface.hh +++ b/src/physics/base/PhysicsInterface.hh @@ -170,7 +170,7 @@ struct PhysicsParamsData //// USER-CONFIGURABLE CONSTANTS //// real_type scaling_min_range{}; //!< rho [cm] real_type scaling_fraction{}; //!< alpha [unitless] - // real_type max_eloss_fraction{}; //!< For scaled range calculation + real_type linear_loss_limit{}; //!< For scaled range calculation //// MEMBER FUNCTIONS //// @@ -178,7 +178,8 @@ struct PhysicsParamsData explicit CELER_FUNCTION operator bool() const { return !process_groups.empty() && max_particle_processes - && scaling_min_range > 0 && scaling_fraction > 0; + && scaling_min_range > 0 && scaling_fraction > 0 + && linear_loss_limit > 0; } //! Assign from another set of data @@ -201,6 +202,7 @@ struct PhysicsParamsData scaling_min_range = other.scaling_min_range; scaling_fraction = other.scaling_fraction; + linear_loss_limit = other.linear_loss_limit; return *this; } @@ -218,9 +220,12 @@ struct PhysicsParamsData */ struct PhysicsTrackState { - real_type interaction_mfp; //!< Remaining MFP to interaction - real_type step_length; //!< Maximum step length - real_type macro_xs; + real_type interaction_mfp; //!< Remaining MFP to interaction + + real_type step_length; //!< Overall physics step length + real_type range_limit; //!< Maximum step due to eloss + real_type macro_xs; //!< Total cross section + ModelId model_id; //!< Selected model if interacting ElementComponentId element_id; //!< Selected element during interaction }; diff --git a/src/physics/base/PhysicsParams.cc b/src/physics/base/PhysicsParams.cc index c4185d48a3..925f5e0dcb 100644 --- a/src/physics/base/PhysicsParams.cc +++ b/src/physics/base/PhysicsParams.cc @@ -126,8 +126,12 @@ void PhysicsParams::build_options(const Options& opts, HostValue* data) const "Non-positive max_step_over_range=" << opts.max_step_over_range); CELER_VALIDATE(opts.min_range > 0, "Non-positive min_range=" << opts.min_range); + CELER_VALIDATE( + opts.linear_loss_limit >= 0 && opts.linear_loss_limit <= 1, + "Non-fractional linear_loss_limit=" << opts.linear_loss_limit); data->scaling_min_range = opts.min_range; data->scaling_fraction = opts.max_step_over_range; + data->linear_loss_limit = opts.linear_loss_limit; } //---------------------------------------------------------------------------// diff --git a/src/physics/base/PhysicsParams.hh b/src/physics/base/PhysicsParams.hh index 066febcf35..b6becff256 100644 --- a/src/physics/base/PhysicsParams.hh +++ b/src/physics/base/PhysicsParams.hh @@ -40,6 +40,9 @@ class ParticleParams; * - \c max_step_over_range: at higher energy (longer range), gradually * decrease the maximum step length until it's this fraction of the tabulated * range. + * - \c linear_loss_limit: if the mean energy loss along a step is greater than + * this fractional value of the pre-step kinetic energy, recalculate the + * energy loss. */ class PhysicsParams { @@ -62,6 +65,7 @@ class PhysicsParams { real_type min_range = 1 * units::millimeter; //!< rho_R real_type max_step_over_range = 0.2; //!< alpha_r + real_type linear_loss_limit = 0.01; //!< xi }; //! Physics parameter construction arguments diff --git a/src/physics/base/PhysicsStepUtils.hh b/src/physics/base/PhysicsStepUtils.hh index 397bdb8403..1cfc08ec28 100644 --- a/src/physics/base/PhysicsStepUtils.hh +++ b/src/physics/base/PhysicsStepUtils.hh @@ -17,15 +17,15 @@ namespace celeritas //---------------------------------------------------------------------------// // INLINE HELPER FUNCTIONS //---------------------------------------------------------------------------// -inline CELER_FUNCTION real_type -calc_tabulated_physics_step(const MaterialTrackView& material, - const ParticleTrackView& particle, - PhysicsTrackView& physics); +inline CELER_FUNCTION void +update_physics_step(const MaterialTrackView& material, + const ParticleTrackView& particle, + PhysicsTrackView& physics); -inline CELER_FUNCTION real_type -calc_energy_loss(const ParticleTrackView& particle, - const PhysicsTrackView& physics, - real_type step_length); +inline CELER_FUNCTION ParticleTrackView::Energy + calc_energy_loss(const ParticleTrackView& particle, + const PhysicsTrackView& physics, + real_type step_length); template inline CELER_FUNCTION ModelId select_model(const ParticleTrackView& particle, diff --git a/src/physics/base/PhysicsStepUtils.i.hh b/src/physics/base/PhysicsStepUtils.i.hh index 182e5559b3..65f30fc969 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -10,21 +10,22 @@ #include "base/Algorithms.hh" #include "base/NumericLimits.hh" #include "base/Range.hh" +#include "physics/grid/EnergyLossCalculator.hh" +#include "physics/grid/InverseRangeCalculator.hh" +#include "physics/grid/RangeCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "Types.hh" namespace celeritas { //---------------------------------------------------------------------------// /*! - * Calculate physics steps based on cross sections and range limits. + * Calculate physics step limits based on cross sections and range limiters. */ -CELER_FUNCTION real_type -calc_tabulated_physics_step(const MaterialTrackView& material, - const ParticleTrackView& particle, - PhysicsTrackView& physics) +CELER_FUNCTION void update_physics_step(const MaterialTrackView& material, + const ParticleTrackView& particle, + PhysicsTrackView& physics) { - CELER_EXPECT(physics.has_interaction_mfp()); - constexpr real_type inf = numeric_limits::infinity(); using VGT = ValueGridType; @@ -32,10 +33,8 @@ calc_tabulated_physics_step(const MaterialTrackView& material, // type) and calculate cross section and particle range. real_type total_macro_xs = 0; real_type min_range = inf; - for (ParticleProcessId::value_type pp_idx : - range(physics.num_particle_processes())) + for (auto ppid : range(ParticleProcessId{physics.num_particle_processes()})) { - const ParticleProcessId ppid{pp_idx}; real_type process_xs = 0; if (auto model_id = physics.hardwired_model(ppid, particle.energy())) { @@ -51,7 +50,7 @@ calc_tabulated_physics_step(const MaterialTrackView& material, // Calculate macroscopic cross section for this process, then // accumulate it into the total cross section and save the cross // section for later. - auto calc_xs = physics.make_calculator(grid_id); + auto calc_xs = physics.make_calculator(grid_id); process_xs = calc_xs(particle.energy()); total_macro_xs += process_xs; } @@ -59,13 +58,12 @@ calc_tabulated_physics_step(const MaterialTrackView& material, if (auto grid_id = physics.value_grid(VGT::range, ppid)) { - auto calc_range = physics.make_calculator(grid_id); + auto calc_range = physics.make_calculator(grid_id); real_type process_range = calc_range(particle.energy()); - // TODO: scale range by sqrt(particle.energy() / minKE) - // if < minKE?? min_range = min(min_range, process_range); } } + physics.macro_xs(total_macro_xs); if (min_range != inf) { @@ -73,43 +71,107 @@ calc_tabulated_physics_step(const MaterialTrackView& material, // user options min_range = physics.range_to_step(min_range); } - - // Update step length with discrete interaction - return min(min_range, physics.interaction_mfp() / total_macro_xs); + physics.range_limit(min_range); } //---------------------------------------------------------------------------// /*! - * Calculate energy loss over the given "true" step length. + * Calculate mean energy loss over the given "true" step length. + * + * See section 7.2.4 Run Time Energy Loss Computation of the Geant4 physics + * manual. See also the longer discussions in section 8 of PHYS010 of the + * Geant3 manual (1993). + * + * In Geant4's calculation \c G4VEnergyLossProcess::AlongStepDoIt: + * - \c reduceFactor is 1 except for heavy ions + * - \c massRatio is 1 except for heavy ions + * - \c length = step + * - \c fRange = physics.range_limit() + * - \c linLossLimit = \c linear_loss_limit = \f$xi\f$ + * + * Energy loss rate is a function of differential cross section: the integral + * of low-energy secondaries (below \c T) produced as a function of energy: \f[ + * \frac{dE}{dx} = N_Z \int_0^T \frac{d \sigma_Z(E, T)}{dT} T dT + * \f] + * + * The stopping range \em R due to these low-energy processes is: + * \f[ + * R = \int_0 ^{E_0} - \frac{dx}{dE} dE . + * \f] + * + * Both Celeritas and Geant4 approximate the range limit as the minimum range + * over all processes, rather than the range as a result from integrating all + * energy loss processes over the allowed energy range. + * + * Geant4's stepping algorithm independently stores \c fRange for each process, + * then (looping externally over all processes) calculates energy loss, checks + * for the linear loss limit, and reduces the particle energy. Celeritas + * inverts this loop so the total energy loss from along-step processess (not + * including multiple scattering) is calculated first, then checked against + * being greater than the lineaer loss limit. + * + * \todo The geant3 manual makes the point that linear interpolation of energy + * loss rate results in a piecewise constant energy deposition curve, which is + * why they use spline interpolation. Investigate higher-order reconstruction + * of energy loss curve, e.g. through spline-based interpolation. + * + * \todo Range should be a integral of a function of energy loss rate summed + * over all processes, not independently considered. */ -CELER_FUNCTION real_type calc_energy_loss(const ParticleTrackView& particle, - const PhysicsTrackView& physics, - real_type step) +CELER_FUNCTION ParticleTrackView::Energy + calc_energy_loss(const ParticleTrackView& particle, + const PhysicsTrackView& physics, + real_type step) { CELER_EXPECT(step >= 0); + static_assert(ParticleTrackView::Energy::unit_type::value() + == EnergyLossCalculator::Energy::unit_type::value(), + "Incompatible energy types"); using VGT = ValueGridType; + const auto pre_step_energy = particle.energy(); - // Loop over all processes that apply to this track (based on particle - // type) and calculate cross section and particle range. + // Calculate the sum of energy loss rate over all processes. real_type total_eloss_rate = 0; - for (ParticleProcessId::value_type pp_idx : - range(physics.num_particle_processes())) + for (auto ppid : range(ParticleProcessId{physics.num_particle_processes()})) { - const ParticleProcessId ppid{pp_idx}; if (auto grid_id = physics.value_grid(VGT::energy_loss, ppid)) { - auto calc_eloss_rate = physics.make_calculator(grid_id); - total_eloss_rate += calc_eloss_rate(particle.energy()); + auto calc_eloss_rate + = physics.make_calculator(grid_id); + total_eloss_rate += calc_eloss_rate(pre_step_energy); } } - // TODO: reduce energy loss rate using range tables for individual - // processes?? aka "long step" with max_eloss_fraction. Unlike Geant4 where - // each process sequentially operates on the track, we can apply limits to - // all energy loss processes simultaneouly. + // Scale loss rate by step length + real_type eloss = total_eloss_rate * step; + + if (eloss > pre_step_energy.value() * physics.linear_loss_limit()) + { + // Enough energy is lost over this step that the dE/dx linear + // approximation is probably wrong. Reduce the energy loss rate + // using range tables for individual processes. + const real_type remaining_range = physics.range_limit() - step; + + // Find the highest energy such that the post-step range limiter is the + // same point as the pre-step estimate of the range. + real_type efinal = max(0, pre_step_energy.value() - eloss); + for (auto ppid : + range(ParticleProcessId{physics.num_particle_processes()})) + { + if (auto grid_id = physics.value_grid(VGT::range, ppid)) + { + auto calc_energy + = physics.make_calculator(grid_id); + efinal = max(efinal, calc_energy(remaining_range).value()); + } + } + + eloss = pre_step_energy.value() - efinal; + } - return total_eloss_rate * step; + CELER_ENSURE(eloss <= pre_step_energy.value()); + return ParticleTrackView::Energy{eloss}; } //---------------------------------------------------------------------------// diff --git a/src/physics/base/PhysicsTrackView.hh b/src/physics/base/PhysicsTrackView.hh index 35b9ec6bfd..ec843f794d 100644 --- a/src/physics/base/PhysicsTrackView.hh +++ b/src/physics/base/PhysicsTrackView.hh @@ -12,7 +12,6 @@ #include "base/Types.hh" #include "physics/base/Units.hh" #include "physics/grid/GridIdFinder.hh" -#include "physics/grid/PhysicsGridCalculator.hh" #include "physics/material/MaterialView.hh" #include "physics/material/Types.hh" #include "Types.hh" @@ -54,10 +53,13 @@ class PhysicsTrackView // Set the remaining MFP to interaction inline CELER_FUNCTION void interaction_mfp(real_type); - // Set the physics step length + // Set the overall physics step length inline CELER_FUNCTION void step_length(real_type); - // Set the total (process-integrated) macroscopic xs + // Set the maximum range from along-step energy loss [cm] + inline CELER_FUNCTION void range_limit(real_type); + + // Set the total (process-integrated) macroscopic xs [cm^-1] inline CELER_FUNCTION void macro_xs(real_type); // Select a model for the current interaction (or {} for no interaction) @@ -74,6 +76,9 @@ class PhysicsTrackView // Maximum step length [cm] CELER_FORCEINLINE_FUNCTION real_type step_length() const; + // Maximum range from along-step energy loss [cm] + CELER_FORCEINLINE_FUNCTION real_type range_limit() const; + // Total (process-integrated) macroscopic xs [cm^-1] CELER_FORCEINLINE_FUNCTION real_type macro_xs() const; @@ -106,14 +111,17 @@ class PhysicsTrackView // Calculate scaled step range inline CELER_FUNCTION real_type range_to_step(real_type range) const; + // Fractional energy loss allowed before post-step recalculation + inline CELER_FUNCTION real_type linear_loss_limit() const; + // Calculate macroscopic cross section on the fly for the given model inline CELER_FUNCTION real_type calc_xs_otf(ModelId model, MaterialView& material, MevEnergy energy) const; // Construct a grid calculator from a physics table - inline CELER_FUNCTION - PhysicsGridCalculator make_calculator(ValueGridId) const; + template + inline CELER_FUNCTION T make_calculator(ValueGridId) const; //// SCRATCH SPACE //// diff --git a/src/physics/base/PhysicsTrackView.i.hh b/src/physics/base/PhysicsTrackView.i.hh index c3b7bc8570..39d1200d7a 100644 --- a/src/physics/base/PhysicsTrackView.i.hh +++ b/src/physics/base/PhysicsTrackView.i.hh @@ -41,6 +41,7 @@ PhysicsTrackView::operator=(const Initializer_t&) { this->state().interaction_mfp = -1; this->state().step_length = -1; + this->state().range_limit = -1; this->state().macro_xs = -1; this->state().model_id = ModelId{}; return *this; @@ -78,6 +79,16 @@ CELER_FUNCTION void PhysicsTrackView::macro_xs(real_type inv_distance) this->state().macro_xs = inv_distance; } +//---------------------------------------------------------------------------// +/*! + * Set the maximum range from along-step energy loss. + */ +CELER_FUNCTION void PhysicsTrackView::range_limit(real_type distance) +{ + CELER_EXPECT(distance > 0); + this->state().range_limit = distance; +} + //---------------------------------------------------------------------------// /*! * Select a model ID for the current track. @@ -113,6 +124,10 @@ CELER_FUNCTION real_type PhysicsTrackView::interaction_mfp() const //---------------------------------------------------------------------------// /*! * Maximum step length. + * + * \todo Maybe calculate on the fly, or set at the same time as range + * limit/macro xs? Should be: + * min(this->range_limit(), this->interaction_mfp() / this->macro_xs()) */ CELER_FUNCTION real_type PhysicsTrackView::step_length() const { @@ -135,6 +150,17 @@ CELER_FUNCTION real_type PhysicsTrackView::macro_xs() const return xs; } +//---------------------------------------------------------------------------// +/*! + * Maximum range from along-step energy loss. + */ +CELER_FUNCTION real_type PhysicsTrackView::range_limit() const +{ + real_type limit = this->state().range_limit; + CELER_ENSURE(limit >= 0); + return limit; +} + //---------------------------------------------------------------------------// /*! * Access the model ID that has been selected for the current track. @@ -172,7 +198,7 @@ CELER_FUNCTION ProcessId PhysicsTrackView::process(ParticleProcessId ppid) const * Return value grid data for the given table type and process if available. * * If the result is not null, it can be used to instantiate a - * PhysicsGridCalculator. + * grid Calculator. * * If the result is null, it's likely because the process doesn't have the * associated value (e.g. if the table type is "energy_loss" and the process is @@ -262,6 +288,15 @@ CELER_FUNCTION real_type PhysicsTrackView::range_to_step(real_type range) const return range; } +//---------------------------------------------------------------------------// +/*! + * Fractional along-step energy loss allowed before recalculating from range. + */ +CELER_FUNCTION real_type PhysicsTrackView::linear_loss_limit() const +{ + return params_.linear_loss_limit; +} + //---------------------------------------------------------------------------// /*! * Calculate macroscopic cross section on the fly. @@ -290,20 +325,21 @@ CELER_FUNCTION real_type PhysicsTrackView::calc_xs_otf(ModelId model, //---------------------------------------------------------------------------// /*! - * Access data for constructing PhysicsGridCalculator. + * Construct a grid calculator of the given type. + * + * The calculator must take two arguments: a reference to XsGridData, and a + * reference to the Values data structure. */ -CELER_FUNCTION PhysicsGridCalculator -PhysicsTrackView::make_calculator(ValueGridId id) const +template +CELER_FUNCTION T PhysicsTrackView::make_calculator(ValueGridId id) const { CELER_EXPECT(id < params_.value_grids.size()); - return PhysicsGridCalculator(params_.value_grids[id], params_.reals); + return T{params_.value_grids[id], params_.reals}; } //---------------------------------------------------------------------------// /*! * Access scratch space for particle-process cross section calculations. - * - * \todo Try changing the ordering to coalesce memory for GPU access. */ CELER_FUNCTION real_type& PhysicsTrackView::per_process_xs(ParticleProcessId ppid) diff --git a/src/physics/grid/ValueGridBuilder.cc b/src/physics/grid/ValueGridBuilder.cc index 4ae7c3e9e0..ba09c43482 100644 --- a/src/physics/grid/ValueGridBuilder.cc +++ b/src/physics/grid/ValueGridBuilder.cc @@ -159,12 +159,14 @@ auto ValueGridXsBuilder::build(ValueGridInserter insert) const -> ValueGridId */ ValueGridLogBuilder::ValueGridLogBuilder(real_type emin, real_type emax, - VecReal xs) - : log_emin_(std::log(emin)), log_emax_(std::log(emax)), xs_(std::move(xs)) + VecReal value) + : log_emin_(std::log(emin)) + , log_emax_(std::log(emax)) + , value_(std::move(value)) { CELER_EXPECT(emin > 0); CELER_EXPECT(emax > emin); - CELER_EXPECT(xs_.size() >= 2); + CELER_EXPECT(value_.size() >= 2); } //---------------------------------------------------------------------------// @@ -174,8 +176,22 @@ ValueGridLogBuilder::ValueGridLogBuilder(real_type emin, auto ValueGridLogBuilder::build(ValueGridInserter insert) const -> ValueGridId { return insert( - UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()), - make_span(xs_)); + UniformGridData::from_bounds(log_emin_, log_emax_, value_.size()), + this->value()); +} + +//---------------------------------------------------------------------------// +// RANGE BUILDER +//---------------------------------------------------------------------------// +/*! + * Construct from raw data. + */ +ValueGridRangeBuilder::ValueGridRangeBuilder(real_type emin, + real_type emax, + VecReal xs) + : Base(emin, emax, std::move(xs)) +{ + CELER_EXPECT(is_monotonic_increasing(this->value())); } //---------------------------------------------------------------------------// diff --git a/src/physics/grid/ValueGridBuilder.hh b/src/physics/grid/ValueGridBuilder.hh index 6fa89bb200..f867a7f21a 100644 --- a/src/physics/grid/ValueGridBuilder.hh +++ b/src/physics/grid/ValueGridBuilder.hh @@ -84,13 +84,14 @@ class ValueGridXsBuilder final : public ValueGridBuilder * * This vector is still uniform in log(E). */ -class ValueGridLogBuilder final : public ValueGridBuilder +class ValueGridLogBuilder : public ValueGridBuilder { public: //!@{ //! Type aliases - using VecReal = std::vector; - using Id = PieId; + using VecReal = std::vector; + using SpanConstReal = Span; + using Id = PieId; //!@} public: @@ -100,10 +101,29 @@ class ValueGridLogBuilder final : public ValueGridBuilder // Construct in the given store ValueGridId build(ValueGridInserter) const final; + //! Access Values + SpanConstReal value() const { return make_span(value_); } + private: real_type log_emin_; real_type log_emax_; - VecReal xs_; + VecReal value_; +}; + +//---------------------------------------------------------------------------// +/*! + * Build a physics vector for range tables. + * + * Range tables are uniform in log(E), and range must monotonically increase + * with energy. + */ +class ValueGridRangeBuilder : public ValueGridLogBuilder +{ + using Base = ValueGridLogBuilder; + + public: + // Construct + ValueGridRangeBuilder(real_type emin, real_type emax, VecReal value); }; //---------------------------------------------------------------------------// diff --git a/test/physics/base/MockProcess.cc b/test/physics/base/MockProcess.cc index 35bc4a6243..078a8115e0 100644 --- a/test/physics/base/MockProcess.cc +++ b/test/physics/base/MockProcess.cc @@ -23,7 +23,6 @@ MockProcess::MockProcess(Input data) : data_(std::move(data)) CELER_EXPECT(data_.interact); CELER_EXPECT(data_.xs >= celeritas::zero_quantity()); CELER_EXPECT(data_.energy_loss >= 0); - CELER_EXPECT(data_.range >= 0); } //---------------------------------------------------------------------------// @@ -52,28 +51,26 @@ auto MockProcess::step_limits(Applicability range) const -> StepLimitBuilders StepLimitBuilders builders; if (data_.xs > zero_quantity()) { - real_type value = unit_cast(data_.xs) * numdens; + real_type xs = unit_cast(data_.xs) * numdens; builders[size_type(ValueGridType::macro_xs)] - = std::make_unique(range.lower.value(), - range.upper.value(), - VecReal{value, value}); + = std::make_unique( + range.lower.value(), range.upper.value(), VecReal{xs, xs}); } if (data_.energy_loss > 0) { - real_type value = data_.energy_loss * numdens; + real_type eloss_rate = data_.energy_loss * numdens; builders[size_type(ValueGridType::energy_loss)] - = std::make_unique(range.lower.value(), - range.upper.value(), - VecReal{value, value}); - } - if (data_.range > 0) - { - builders[size_type(ValueGridType::range)] = std::make_unique( range.lower.value(), range.upper.value(), - VecReal{data_.range * range.lower.value(), - data_.range * range.upper.value()}); + VecReal{eloss_rate, eloss_rate}); + + builders[size_type(ValueGridType::range)] + = std::make_unique( + range.lower.value(), + range.upper.value(), + VecReal{range.lower.value() / eloss_rate, + range.upper.value() / eloss_rate}); } return builders; diff --git a/test/physics/base/MockProcess.hh b/test/physics/base/MockProcess.hh index 11b9c73b5b..0d6b3c3674 100644 --- a/test/physics/base/MockProcess.hh +++ b/test/physics/base/MockProcess.hh @@ -27,7 +27,8 @@ namespace celeritas_test * constant with energy. * - Energy loss rate is also constant with energy and scales with the number * density. - * - Range scales linearly with energy. + * - Range is determined by the energy loss rate -- constant energy loss rate k + * means linear range of E/k. * * The given applicability vector has one element per model that it will * create. Each model can have a different particle type and/or energy range. @@ -53,8 +54,7 @@ class MockProcess : public celeritas::Process VecApplicability applic; //!< Applicablity per model ModelCallback interact; //!< MockModel::interact callback BarnMicroXs xs{}; //!< Constant per atom [bn] - real_type energy_loss{}; //!< Constant per atom [MeV/cm / (1/cm^3)] - real_type range{}; //!< Linearly increasing [cm/MeV] + real_type energy_loss{}; //!< Constant per atom [MeV/cm / cm^-3] }; public: diff --git a/test/physics/base/Physics.test.cc b/test/physics/base/Physics.test.cc index 918a7a5cad..d6f7d09024 100644 --- a/test/physics/base/Physics.test.cc +++ b/test/physics/base/Physics.test.cc @@ -12,6 +12,8 @@ #include "base/Range.hh" #include "base/PieStateStore.hh" #include "physics/base/ParticleParams.hh" +#include "physics/grid/RangeCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "PhysicsTestBase.hh" #include "Physics.test.hh" @@ -276,7 +278,7 @@ TEST_F(PhysicsTrackViewHostTest, calc_xs) auto scat_ppid = this->find_ppid(phys, "scattering"); auto id = phys.value_grid(ValueGridType::macro_xs, scat_ppid); ASSERT_TRUE(id); - auto calc_xs = phys.make_calculator(id); + auto calc_xs = phys.make_calculator(id); xs.push_back(calc_xs(MevEnergy{1.0})); } } @@ -302,7 +304,7 @@ TEST_F(PhysicsTrackViewHostTest, calc_range) ASSERT_TRUE(meow_ppid); auto id = phys.value_grid(ValueGridType::range, meow_ppid); ASSERT_TRUE(id); - auto calc_range = phys.make_calculator(id); + auto calc_range = phys.make_calculator(id); for (real_type energy : {0.01, 1.0, 1e2}) { range.push_back(calc_range(MevEnergy{energy})); @@ -310,8 +312,9 @@ TEST_F(PhysicsTrackViewHostTest, calc_range) } } - const double expected_range[] = {0.04, 4, 40, 0.04, 4, 40}; - const double expected_step[] = {0.04, 0.958, 8.1598, 0.04, 0.958, 8.1598}; + const double expected_range[] = {0.025, 2.5, 25, 0.025, 2.5, 25}; + const double expected_step[] + = {0.025, 0.6568, 5.15968, 0.025, 0.6568, 5.15968}; EXPECT_VEC_SOFT_EQ(expected_range, range); EXPECT_VEC_SOFT_EQ(expected_step, step); } @@ -345,17 +348,16 @@ TEST_F(PhysicsTrackViewHostTest, cuda_surrogate) } } - // PRINT_EXPECTED(step); const double expected_step[] = {166.6666666667, 166.6666666667, 166.6666666667, 166.6666666667, inf, - 0.003, - 0.003, - 0.7573333333333, - 8.1598, - 8.1598}; + 2.5e-05, + 0.00025, + 0.178, + 0.6568, + 0.6568}; EXPECT_VEC_SOFT_EQ(expected_step, step); } diff --git a/test/physics/base/Physics.test.hh b/test/physics/base/Physics.test.hh index a913628560..26048d0336 100644 --- a/test/physics/base/Physics.test.hh +++ b/test/physics/base/Physics.test.hh @@ -13,6 +13,8 @@ #include "physics/base/PhysicsInterface.hh" #include "physics/base/Units.hh" #include "physics/base/Types.hh" +#include "physics/grid/RangeCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "physics/material/Types.hh" // Kernel functions @@ -62,7 +64,7 @@ inline CELER_FUNCTION celeritas::real_type real_type process_xs = 0; if (auto id = phys.value_grid(ValueGridType::macro_xs, ppid)) { - auto calc_xs = phys.make_calculator(id); + auto calc_xs = phys.make_calculator(id); process_xs = calc_xs(energy); } @@ -85,7 +87,7 @@ inline CELER_FUNCTION celeritas::real_type { if (auto id = phys.value_grid(ValueGridType::range, ppid)) { - auto calc_range = phys.make_calculator(id); + auto calc_range = phys.make_calculator(id); step = min(step, calc_range(energy)); } } diff --git a/test/physics/base/PhysicsStepUtils.test.cc b/test/physics/base/PhysicsStepUtils.test.cc index 280fe05c4a..01207a608a 100644 --- a/test/physics/base/PhysicsStepUtils.test.cc +++ b/test/physics/base/PhysicsStepUtils.test.cc @@ -30,6 +30,18 @@ class PhysicsStepUtilsTest : public PhysicsTestBase using ParticleStateStore = PieStateStore; using PhysicsStateStore = PieStateStore; + PhysicsOptions build_physics_options() const override + { + PhysicsOptions opts; + if (false) + { + // Don't scale the range -- use exactly the analytic values our + // model has. + opts.min_range = inf; + } + return opts; + } + void SetUp() override { Base::SetUp(); @@ -40,6 +52,30 @@ class PhysicsStepUtilsTest : public PhysicsTestBase phys_state = PhysicsStateStore(*this->physics(), 1); } + PhysicsTrackView init_track(MaterialTrackView* mat, + MaterialId mid, + ParticleTrackView* par, + const char* name, + MevEnergy energy) + { + CELER_EXPECT(mat && par); + CELER_EXPECT(mid < this->materials()->size()); + *mat = MaterialTrackView::Initializer_t{mid}; + + ParticleTrackView::Initializer_t par_init; + par_init.particle_id = this->particles()->find(name); + CELER_EXPECT(par_init.particle_id); + par_init.energy = energy; + *par = par_init; + + PhysicsTrackView phys(this->physics()->host_pointers(), + phys_state.ref(), + par->particle_id(), + mat->material_id(), + ThreadId{0}); + return phys; + } + MaterialStateStore mat_state; ParticleStateStore par_state; PhysicsStateStore phys_state; @@ -49,88 +85,92 @@ class PhysicsStepUtilsTest : public PhysicsTestBase // TESTS //---------------------------------------------------------------------------// -TEST_F(PhysicsStepUtilsTest, calc_tabulated_physics_step) +TEST_F(PhysicsStepUtilsTest, update_physics_step) { MaterialTrackView material( this->materials()->host_pointers(), mat_state.ref(), ThreadId{0}); ParticleTrackView particle( this->particles()->host_pointers(), par_state.ref(), ThreadId{0}); - MaterialTrackView::Initializer_t mat_init; - ParticleTrackView::Initializer_t par_init; - // Test a variety of energy ranges and multiple material IDs + // Test a variety of energies and multiple material IDs { - mat_init.material_id = MaterialId{0}; - par_init.energy = MevEnergy{1}; - par_init.particle_id = this->particles()->find("gamma"); - material = mat_init; - particle = par_init; - PhysicsTrackView phys(this->physics()->host_pointers(), - phys_state.ref(), - par_init.particle_id, - mat_init.material_id, - ThreadId{0}); - phys.interaction_mfp(1); - real_type step - = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1. / 3.e-4, step); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); + celeritas::update_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(3.e-4, phys.macro_xs()); + EXPECT_SOFT_EQ(inf, phys.range_limit()); } { - mat_init.material_id = MaterialId{1}; - par_init.energy = MevEnergy{10}; - par_init.particle_id = this->particles()->find("celeriton"); - material = mat_init; - particle = par_init; - PhysicsTrackView phys(this->physics()->host_pointers(), - phys_state.ref(), - par_init.particle_id, - mat_init.material_id, - ThreadId{0}); - phys.interaction_mfp(1e-2); - real_type step - = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1.e-2 / 9.e-3, step); - - // Increase the distance to interaction so range limits the step length - phys.interaction_mfp(1); - step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(4.1595999999999984, step); - - // Decrease the particle energy - par_init.energy = MevEnergy{1e-2}; - particle = par_init; - step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(2.e-2, step); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{1}, &particle, "celeriton", MevEnergy{10}); + celeritas::update_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(9.e-3, phys.macro_xs()); + EXPECT_SOFT_EQ(0.6568, phys.range_limit()); } { - mat_init.material_id = MaterialId{2}; - par_init.energy = MevEnergy{1e-2}; - par_init.particle_id = this->particles()->find("anti-celeriton"); - material = mat_init; - particle = par_init; - PhysicsTrackView phys(this->physics()->host_pointers(), - phys_state.ref(), - par_init.particle_id, - mat_init.material_id, - ThreadId{0}); - phys.interaction_mfp(1e-2); - real_type step - = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1.e-2 / 9.e-1, step); - - // Increase the distance to interaction so range limits the step length - phys.interaction_mfp(1); - step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(0.03, step); - - // Increase the particle energy so interaction limits the step length - par_init.energy = MevEnergy{10}; - particle = par_init; - step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1. / 9.e-1, step); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{1}, &particle, "celeriton", MevEnergy{1e-2}); + celeritas::update_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(0.009, phys.macro_xs()); + EXPECT_SOFT_EQ(0.0025, phys.range_limit()); + } + { + PhysicsTrackView phys = this->init_track(&material, + MaterialId{2}, + &particle, + "anti-celeriton", + MevEnergy{1e-2}); + celeritas::update_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(9.e-1, phys.macro_xs()); + EXPECT_SOFT_EQ(2.5e-5, phys.range_limit()); + } + { + PhysicsTrackView phys = this->init_track(&material, + MaterialId{2}, + &particle, + "anti-celeriton", + MevEnergy{10}); + celeritas::update_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(9.e-1, phys.macro_xs()); + EXPECT_SOFT_EQ(0.025, phys.range_limit()); } } -TEST_F(PhysicsStepUtilsTest, calc_energy_loss) {} +TEST_F(PhysicsStepUtilsTest, calc_energy_loss) +{ + MaterialTrackView material( + this->materials()->host_pointers(), mat_state.ref(), ThreadId{0}); + ParticleTrackView particle( + this->particles()->host_pointers(), par_state.ref(), ThreadId{0}); + + // Test a variety of energies and multiple material IDs + { + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); + phys.macro_xs(3e-4); + phys.range_limit(inf); + + // Long step, but gamma means no energy loss + EXPECT_SOFT_EQ(0, celeritas::calc_energy_loss(particle, phys, 1e4)); + } + { + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "celeriton", MevEnergy{10}); + phys.macro_xs(9e-4); + phys.range_limit(25.0); // limit when min_range=inf + + // Tiny step: should still be linear loss (single process) + const real_type eloss_rate = 0.2 + 0.4; + EXPECT_SOFT_EQ(eloss_rate * 1e-6, + celeritas::calc_energy_loss(particle, phys, 1e-6)); + + // Long step (lose half energy) will call inverse lookup. The correct + // answer (if range table construction was done over energy loss) + // should be half since the slowing down rate is constant over all + const real_type step = 0.5 * particle.energy().value() / eloss_rate; + EXPECT_SOFT_EQ( + 5, celeritas::calc_energy_loss(particle, phys, step).value()); + } +} TEST_F(PhysicsStepUtilsTest, select_model) {} diff --git a/test/physics/base/PhysicsTestBase.cc b/test/physics/base/PhysicsTestBase.cc index e0a07e1a80..b7fa6456a1 100644 --- a/test/physics/base/PhysicsTestBase.cc +++ b/test/physics/base/PhysicsTestBase.cc @@ -76,6 +76,12 @@ auto PhysicsTestBase::build_particles() const -> SPConstParticles return std::make_shared(std::move(inp)); } +//---------------------------------------------------------------------------// +auto PhysicsTestBase::build_physics_options() const -> PhysicsOptions +{ + return {}; +} + //---------------------------------------------------------------------------// auto PhysicsTestBase::build_physics() const -> SPConstPhysics { @@ -83,6 +89,7 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics PhysicsParams::Input physics_inp; physics_inp.materials = this->materials(); physics_inp.particles = this->particles(); + physics_inp.options = this->build_physics_options(); // Add a few processes MockProcess::Input inp; @@ -94,7 +101,6 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics make_applicability("celeriton", 1, 100)}; inp.xs = Barn{1.0}; inp.energy_loss = {}; - inp.range = {}; physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -102,7 +108,6 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics inp.applic = {make_applicability("gamma", 1e-6, 100)}; inp.xs = Barn{2.0}; inp.energy_loss = {}; - inp.range = {}; physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -112,8 +117,7 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics make_applicability("celeriton", 1, 10), make_applicability("celeriton", 10, 100)}; inp.xs = Barn{3.0}; - inp.energy_loss = 0.2; - inp.range = 2; + inp.energy_loss = 0.2 * 1e-20; // 0.2 MeV/cm in celerogen physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -122,8 +126,7 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics inp.applic = {make_applicability("anti-celeriton", 1e-3, 1), make_applicability("anti-celeriton", 1, 100)}; inp.xs = Barn{4.0}; - inp.energy_loss = 0.3; - inp.range = 3; + inp.energy_loss = 0.3 * 1e-20; physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -131,8 +134,7 @@ auto PhysicsTestBase::build_physics() const -> SPConstPhysics inp.applic = {make_applicability("celeriton", 1e-3, 10), make_applicability("anti-celeriton", 1e-3, 10)}; inp.xs = Barn{5.0}; - inp.energy_loss = 0.4; - inp.range = 4; + inp.energy_loss = 0.4 * 1e-20; physics_inp.processes.push_back(std::make_shared(inp)); } return std::make_shared(std::move(physics_inp)); diff --git a/test/physics/base/PhysicsTestBase.hh b/test/physics/base/PhysicsTestBase.hh index 71aae1687c..df45f877a9 100644 --- a/test/physics/base/PhysicsTestBase.hh +++ b/test/physics/base/PhysicsTestBase.hh @@ -11,15 +11,9 @@ #include "physics/material/MaterialParams.hh" #include "physics/base/ParticleParams.hh" +#include "physics/base/PhysicsParams.hh" #include "MockProcess.hh" -namespace celeritas -{ -class MaterialParams; -class ParticleParams; -class PhysicsParams; -} // namespace celeritas - namespace celeritas_test { //---------------------------------------------------------------------------// @@ -45,6 +39,7 @@ class PhysicsTestBase : public celeritas::Test using SPConstMaterials = std::shared_ptr; using SPConstParticles = std::shared_ptr; using SPConstPhysics = std::shared_ptr; + using PhysicsOptions = celeritas::PhysicsParams::Options; using Applicability = celeritas::Applicability; using ModelId = celeritas::ModelId; using ModelCallback = std::function; @@ -57,6 +52,7 @@ class PhysicsTestBase : public celeritas::Test virtual SPConstMaterials build_materials() const; virtual SPConstParticles build_particles() const; + virtual PhysicsOptions build_physics_options() const; virtual SPConstPhysics build_physics() const; const SPConstMaterials& materials() const { return materials_; } From b6ef9e9991fd4d1f2f0eb3711597f72494e4fe4f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 21:42:31 -0500 Subject: [PATCH 05/10] Fix range calculation --- scripts/build/yuri.cmake | 8 +-- scripts/build/yuri.sh | 7 ++- src/physics/base/PhysicsInterface.hh | 2 - src/physics/base/PhysicsStepUtils.hh | 8 +-- src/physics/base/PhysicsStepUtils.i.hh | 47 +++++++++----- src/physics/base/PhysicsTrackView.hh | 6 -- src/physics/base/PhysicsTrackView.i.hh | 26 -------- test/physics/base/PhysicsStepUtils.test.cc | 72 +++++++++++++--------- 8 files changed, 86 insertions(+), 90 deletions(-) diff --git a/scripts/build/yuri.cmake b/scripts/build/yuri.cmake index 8172722043..2deef31d98 100644 --- a/scripts/build/yuri.cmake +++ b/scripts/build/yuri.cmake @@ -2,15 +2,15 @@ macro(set_cache_var var type val) set(${var} "${val}" CACHE "${type}" "yuri.sh") endmacro() -set_cache_var(CELERITAS_BUILD_DOCS BOOL OFF) +set_cache_var(CELERITAS_BUILD_DOCS BOOL ON) # Dependency options set_cache_var(CELERITAS_USE_CUDA BOOL OFF) set_cache_var(CELERITAS_USE_Geant4 BOOL OFF) set_cache_var(CELERITAS_GIT_SUBMODULE BOOL OFF) -set_cache_var(CELERITAS_USE_MPI BOOL OFF) -set_cache_var(CELERITAS_USE_ROOT BOOL OFF) -set_cache_var(CELERITAS_USE_SWIG_Python BOOL OFF) +set_cache_var(CELERITAS_USE_MPI BOOL ON) +set_cache_var(CELERITAS_USE_ROOT BOOL ON) +set_cache_var(CELERITAS_USE_SWIG_Python BOOL ON) set_cache_var(CELERITAS_USE_VecGeom BOOL OFF) # Build options diff --git a/scripts/build/yuri.sh b/scripts/build/yuri.sh index 3774f0a7a0..5750cd06fe 100755 --- a/scripts/build/yuri.sh +++ b/scripts/build/yuri.sh @@ -4,8 +4,8 @@ SOURCE_DIR="$(cd "${BUILDSCRIPT_DIR}" && git rev-parse --show-toplevel)" printf "\e[2;37mBuilding from ${SOURCE_DIR}\e[0m\n" cd $SOURCE_DIR -mkdir build-xcode 2>/dev/null || true -cd build-xcode +mkdir build 2>/dev/null || true +cd build module load doxygen swig @@ -13,7 +13,8 @@ CELERITAS_ENV=${SPACK_ROOT}/var/spack/environments/celeritas/.spack-env/view export PATH=$CELERITAS_ENV/bin:${PATH} export CMAKE_PREFIX_PATH=$CELERITAS_ENV:${CMAKE_PREFIX_PATH} -cmake -C ${BUILDSCRIPT_DIR}/yuri.cmake -G Xcode \ +cmake -C ${BUILDSCRIPT_DIR}/yuri.cmake -G Ninja \ + -DCMAKE_INSTALL_PREFIX:PATH=$SOURCE_DIR/install \ .. ninja -v ctest --output-on-failure diff --git a/src/physics/base/PhysicsInterface.hh b/src/physics/base/PhysicsInterface.hh index cb2404df77..3a204879a9 100644 --- a/src/physics/base/PhysicsInterface.hh +++ b/src/physics/base/PhysicsInterface.hh @@ -221,9 +221,7 @@ struct PhysicsParamsData struct PhysicsTrackState { real_type interaction_mfp; //!< Remaining MFP to interaction - real_type step_length; //!< Overall physics step length - real_type range_limit; //!< Maximum step due to eloss real_type macro_xs; //!< Total cross section ModelId model_id; //!< Selected model if interacting diff --git a/src/physics/base/PhysicsStepUtils.hh b/src/physics/base/PhysicsStepUtils.hh index 1cfc08ec28..bfc6f5c812 100644 --- a/src/physics/base/PhysicsStepUtils.hh +++ b/src/physics/base/PhysicsStepUtils.hh @@ -17,10 +17,10 @@ namespace celeritas //---------------------------------------------------------------------------// // INLINE HELPER FUNCTIONS //---------------------------------------------------------------------------// -inline CELER_FUNCTION void -update_physics_step(const MaterialTrackView& material, - const ParticleTrackView& particle, - PhysicsTrackView& physics); +inline CELER_FUNCTION real_type +calc_tabulated_physics_step(const MaterialTrackView& material, + const ParticleTrackView& particle, + PhysicsTrackView& physics); inline CELER_FUNCTION ParticleTrackView::Energy calc_energy_loss(const ParticleTrackView& particle, diff --git a/src/physics/base/PhysicsStepUtils.i.hh b/src/physics/base/PhysicsStepUtils.i.hh index 65f30fc969..2c1025c8c9 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -22,10 +22,13 @@ namespace celeritas /*! * Calculate physics step limits based on cross sections and range limiters. */ -CELER_FUNCTION void update_physics_step(const MaterialTrackView& material, - const ParticleTrackView& particle, - PhysicsTrackView& physics) +inline CELER_FUNCTION real_type +calc_tabulated_physics_step(const MaterialTrackView& material, + const ParticleTrackView& particle, + PhysicsTrackView& physics) { + CELER_EXPECT(physics.has_interaction_mfp()); + constexpr real_type inf = numeric_limits::infinity(); using VGT = ValueGridType; @@ -71,7 +74,9 @@ CELER_FUNCTION void update_physics_step(const MaterialTrackView& material, // user options min_range = physics.range_to_step(min_range); } - physics.range_limit(min_range); + + // Update step length with discrete interaction + return min(min_range, physics.interaction_mfp() / total_macro_xs); } //---------------------------------------------------------------------------// @@ -113,10 +118,11 @@ CELER_FUNCTION void update_physics_step(const MaterialTrackView& material, * \todo The geant3 manual makes the point that linear interpolation of energy * loss rate results in a piecewise constant energy deposition curve, which is * why they use spline interpolation. Investigate higher-order reconstruction - * of energy loss curve, e.g. through spline-based interpolation. + * of energy loss curve, e.g. through spline-based interpolation or log-log + * interpolation. * - * \todo Range should be a integral of a function of energy loss rate summed - * over all processes, not independently considered. + * \note The inverse range correction assumes range is always the integral of + * the stopping power/energy loss. */ CELER_FUNCTION ParticleTrackView::Energy calc_energy_loss(const ParticleTrackView& particle, @@ -149,25 +155,32 @@ CELER_FUNCTION ParticleTrackView::Energy if (eloss > pre_step_energy.value() * physics.linear_loss_limit()) { // Enough energy is lost over this step that the dE/dx linear - // approximation is probably wrong. Reduce the energy loss rate - // using range tables for individual processes. - const real_type remaining_range = physics.range_limit() - step; - - // Find the highest energy such that the post-step range limiter is the - // same point as the pre-step estimate of the range. - real_type efinal = max(0, pre_step_energy.value() - eloss); + // approximation is probably wrong. Use the definition of the range as + // the integral of 1/loss to back-calculate the actual energy loss + // along the curve given the actual step. + eloss = 0; for (auto ppid : range(ParticleProcessId{physics.num_particle_processes()})) { if (auto grid_id = physics.value_grid(VGT::range, ppid)) { + // Recalculate beginning-of-step range (instead of storing) + auto calc_range + = physics.make_calculator(grid_id); + real_type remaining_range = calc_range(pre_step_energy) - step; + CELER_ASSERT(remaining_range > 0); + + // Calculate energy along the range curve corresponding to the + // actual step taken: this gives the exact energy loss over the + // step due to this process. auto calc_energy = physics.make_calculator(grid_id); - efinal = max(efinal, calc_energy(remaining_range).value()); + eloss += (pre_step_energy.value() + - calc_energy(remaining_range).value()); } } - - eloss = pre_step_energy.value() - efinal; + CELER_ASSERT(eloss > 0); + CELER_ASSERT(eloss < pre_step_energy.value()); } CELER_ENSURE(eloss <= pre_step_energy.value()); diff --git a/src/physics/base/PhysicsTrackView.hh b/src/physics/base/PhysicsTrackView.hh index ec843f794d..911d4124b7 100644 --- a/src/physics/base/PhysicsTrackView.hh +++ b/src/physics/base/PhysicsTrackView.hh @@ -56,9 +56,6 @@ class PhysicsTrackView // Set the overall physics step length inline CELER_FUNCTION void step_length(real_type); - // Set the maximum range from along-step energy loss [cm] - inline CELER_FUNCTION void range_limit(real_type); - // Set the total (process-integrated) macroscopic xs [cm^-1] inline CELER_FUNCTION void macro_xs(real_type); @@ -76,9 +73,6 @@ class PhysicsTrackView // Maximum step length [cm] CELER_FORCEINLINE_FUNCTION real_type step_length() const; - // Maximum range from along-step energy loss [cm] - CELER_FORCEINLINE_FUNCTION real_type range_limit() const; - // Total (process-integrated) macroscopic xs [cm^-1] CELER_FORCEINLINE_FUNCTION real_type macro_xs() const; diff --git a/src/physics/base/PhysicsTrackView.i.hh b/src/physics/base/PhysicsTrackView.i.hh index 39d1200d7a..6f96622a38 100644 --- a/src/physics/base/PhysicsTrackView.i.hh +++ b/src/physics/base/PhysicsTrackView.i.hh @@ -41,7 +41,6 @@ PhysicsTrackView::operator=(const Initializer_t&) { this->state().interaction_mfp = -1; this->state().step_length = -1; - this->state().range_limit = -1; this->state().macro_xs = -1; this->state().model_id = ModelId{}; return *this; @@ -79,16 +78,6 @@ CELER_FUNCTION void PhysicsTrackView::macro_xs(real_type inv_distance) this->state().macro_xs = inv_distance; } -//---------------------------------------------------------------------------// -/*! - * Set the maximum range from along-step energy loss. - */ -CELER_FUNCTION void PhysicsTrackView::range_limit(real_type distance) -{ - CELER_EXPECT(distance > 0); - this->state().range_limit = distance; -} - //---------------------------------------------------------------------------// /*! * Select a model ID for the current track. @@ -124,10 +113,6 @@ CELER_FUNCTION real_type PhysicsTrackView::interaction_mfp() const //---------------------------------------------------------------------------// /*! * Maximum step length. - * - * \todo Maybe calculate on the fly, or set at the same time as range - * limit/macro xs? Should be: - * min(this->range_limit(), this->interaction_mfp() / this->macro_xs()) */ CELER_FUNCTION real_type PhysicsTrackView::step_length() const { @@ -150,17 +135,6 @@ CELER_FUNCTION real_type PhysicsTrackView::macro_xs() const return xs; } -//---------------------------------------------------------------------------// -/*! - * Maximum range from along-step energy loss. - */ -CELER_FUNCTION real_type PhysicsTrackView::range_limit() const -{ - real_type limit = this->state().range_limit; - CELER_ENSURE(limit >= 0); - return limit; -} - //---------------------------------------------------------------------------// /*! * Access the model ID that has been selected for the current track. diff --git a/test/physics/base/PhysicsStepUtils.test.cc b/test/physics/base/PhysicsStepUtils.test.cc index 01207a608a..19135048fe 100644 --- a/test/physics/base/PhysicsStepUtils.test.cc +++ b/test/physics/base/PhysicsStepUtils.test.cc @@ -85,7 +85,7 @@ class PhysicsStepUtilsTest : public PhysicsTestBase // TESTS //---------------------------------------------------------------------------// -TEST_F(PhysicsStepUtilsTest, update_physics_step) +TEST_F(PhysicsStepUtilsTest, calc_tabulated_physics_step) { MaterialTrackView material( this->materials()->host_pointers(), mat_state.ref(), ThreadId{0}); @@ -96,23 +96,30 @@ TEST_F(PhysicsStepUtilsTest, update_physics_step) { PhysicsTrackView phys = this->init_track( &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); - celeritas::update_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(3.e-4, phys.macro_xs()); - EXPECT_SOFT_EQ(inf, phys.range_limit()); + phys.interaction_mfp(1); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(1. / 3.e-4, step); } { PhysicsTrackView phys = this->init_track( &material, MaterialId{1}, &particle, "celeriton", MevEnergy{10}); - celeritas::update_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(9.e-3, phys.macro_xs()); - EXPECT_SOFT_EQ(0.6568, phys.range_limit()); + phys.interaction_mfp(1e-4); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(1.e-4 / 9.e-3, step); + + // Increase the distance to interaction so range limits the step length + phys.interaction_mfp(1); + step = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(0.6568, step); } { PhysicsTrackView phys = this->init_track( &material, MaterialId{1}, &particle, "celeriton", MevEnergy{1e-2}); - celeritas::update_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(0.009, phys.macro_xs()); - EXPECT_SOFT_EQ(0.0025, phys.range_limit()); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(2.5e-3, step); } { PhysicsTrackView phys = this->init_track(&material, @@ -120,9 +127,15 @@ TEST_F(PhysicsStepUtilsTest, update_physics_step) &particle, "anti-celeriton", MevEnergy{1e-2}); - celeritas::update_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(9.e-1, phys.macro_xs()); - EXPECT_SOFT_EQ(2.5e-5, phys.range_limit()); + phys.interaction_mfp(1e-6); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(1.e-6 / 9.e-1, step); + + // Increase the distance to interaction so range limits the step length + phys.interaction_mfp(1); + step = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(2.5e-5, step); } { PhysicsTrackView phys = this->init_track(&material, @@ -130,9 +143,9 @@ TEST_F(PhysicsStepUtilsTest, update_physics_step) &particle, "anti-celeriton", MevEnergy{10}); - celeritas::update_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(9.e-1, phys.macro_xs()); - EXPECT_SOFT_EQ(0.025, phys.range_limit()); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(2.5e-2, step); } } @@ -143,33 +156,36 @@ TEST_F(PhysicsStepUtilsTest, calc_energy_loss) ParticleTrackView particle( this->particles()->host_pointers(), par_state.ref(), ThreadId{0}); - // Test a variety of energies and multiple material IDs { + // Long step, but gamma means no energy loss PhysicsTrackView phys = this->init_track( &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); - phys.macro_xs(3e-4); - phys.range_limit(inf); - - // Long step, but gamma means no energy loss - EXPECT_SOFT_EQ(0, celeritas::calc_energy_loss(particle, phys, 1e4)); + EXPECT_SOFT_EQ( + 0, celeritas::calc_energy_loss(particle, phys, 1e4).value()); } { PhysicsTrackView phys = this->init_track( &material, MaterialId{0}, &particle, "celeriton", MevEnergy{10}); - phys.macro_xs(9e-4); - phys.range_limit(25.0); // limit when min_range=inf + const real_type eloss_rate = 0.2 + 0.4; // Tiny step: should still be linear loss (single process) - const real_type eloss_rate = 0.2 + 0.4; - EXPECT_SOFT_EQ(eloss_rate * 1e-6, - celeritas::calc_energy_loss(particle, phys, 1e-6)); + EXPECT_SOFT_EQ( + eloss_rate * 1e-6, + celeritas::calc_energy_loss(particle, phys, 1e-6).value()); // Long step (lose half energy) will call inverse lookup. The correct // answer (if range table construction was done over energy loss) // should be half since the slowing down rate is constant over all - const real_type step = 0.5 * particle.energy().value() / eloss_rate; + real_type step = 0.5 * particle.energy().value() / eloss_rate; EXPECT_SOFT_EQ( 5, celeritas::calc_energy_loss(particle, phys, step).value()); + + // Long step (lose half energy) will call inverse lookup. The correct + // answer (if range table construction was done over energy loss) + // should be half since the slowing down rate is constant over all + step = 0.999 * particle.energy().value() / eloss_rate; + EXPECT_SOFT_EQ( + 9.99, celeritas::calc_energy_loss(particle, phys, step).value()); } } From 7983de0cfcd7e6156019841ad393aaa7b0d5a0cb Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Mon, 1 Mar 2021 21:59:18 -0500 Subject: [PATCH 06/10] Update documentation --- src/physics/base/PhysicsStepUtils.i.hh | 31 +++++++++++++------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/physics/base/PhysicsStepUtils.i.hh b/src/physics/base/PhysicsStepUtils.i.hh index 2c1025c8c9..366489cacb 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -87,15 +87,10 @@ calc_tabulated_physics_step(const MaterialTrackView& material, * manual. See also the longer discussions in section 8 of PHYS010 of the * Geant3 manual (1993). * - * In Geant4's calculation \c G4VEnergyLossProcess::AlongStepDoIt: - * - \c reduceFactor is 1 except for heavy ions - * - \c massRatio is 1 except for heavy ions - * - \c length = step - * - \c fRange = physics.range_limit() - * - \c linLossLimit = \c linear_loss_limit = \f$xi\f$ - * - * Energy loss rate is a function of differential cross section: the integral - * of low-energy secondaries (below \c T) produced as a function of energy: \f[ + * Energy loss rate (stopping power) is a function of differential cross + * section: the integral of low-energy secondaries (below \c T) produced as a + * function of energy: + * \f[ * \frac{dE}{dx} = N_Z \int_0^T \frac{d \sigma_Z(E, T)}{dT} T dT * \f] * @@ -106,23 +101,29 @@ calc_tabulated_physics_step(const MaterialTrackView& material, * * Both Celeritas and Geant4 approximate the range limit as the minimum range * over all processes, rather than the range as a result from integrating all - * energy loss processes over the allowed energy range. + * energy loss processes over the allowed energy range. This is usuallly not + * a problem in practice because the range will get automatically decreased by + * \c range_to_step . * - * Geant4's stepping algorithm independently stores \c fRange for each process, + * Geant4's stepping algorithm independently stores the range for each process, * then (looping externally over all processes) calculates energy loss, checks * for the linear loss limit, and reduces the particle energy. Celeritas * inverts this loop so the total energy loss from along-step processess (not * including multiple scattering) is calculated first, then checked against - * being greater than the lineaer loss limit. + * being greater than the linear loss limit. + * + * If energy loss is greater than the loss limit, we loop over all + * processes with range tables and recalculate the pre-step range and solve for + * the exact post-step energy loss. + * + * \note The inverse range correction assumes range is always the integral of + * the stopping power/energy loss. * * \todo The geant3 manual makes the point that linear interpolation of energy * loss rate results in a piecewise constant energy deposition curve, which is * why they use spline interpolation. Investigate higher-order reconstruction * of energy loss curve, e.g. through spline-based interpolation or log-log * interpolation. - * - * \note The inverse range correction assumes range is always the integral of - * the stopping power/energy loss. */ CELER_FUNCTION ParticleTrackView::Energy calc_energy_loss(const ParticleTrackView& particle, From ceed24b26c07aad386d2557b8aff0722075736ad Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 2 Mar 2021 08:21:46 -0500 Subject: [PATCH 07/10] Fix missing header and CUDA names --- app/demo-interactor/KNDemoKernel.cu | 4 ++-- app/demo-interactor/KernelUtils.hh | 14 +++++++------- src/physics/grid/ValueGridBuilder.cc | 1 + 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/app/demo-interactor/KNDemoKernel.cu b/app/demo-interactor/KNDemoKernel.cu index bc451de52c..c8aa39cfd8 100644 --- a/app/demo-interactor/KNDemoKernel.cu +++ b/app/demo-interactor/KNDemoKernel.cu @@ -16,7 +16,7 @@ #include "physics/base/SecondaryAllocatorView.hh" #include "physics/em/detail/KleinNishinaInteractor.hh" #include "random/cuda/RngEngine.hh" -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "DetectorView.hh" #include "KernelUtils.hh" @@ -75,7 +75,7 @@ move_kernel(ParamsDeviceRef const params, StateDeviceRef const states) RngEngine rng(states.rng, ThreadId(tid)); // Move to collision - PhysicsGridCalculator calc_xs(params.tables.xs, params.tables.reals); + XsCalculator calc_xs(params.tables.xs, params.tables.reals); demo_interactor::move_to_collision(particle, calc_xs, states.direction[tid], diff --git a/app/demo-interactor/KernelUtils.hh b/app/demo-interactor/KernelUtils.hh index 499ba6ffa3..9d2c28f29c 100644 --- a/app/demo-interactor/KernelUtils.hh +++ b/app/demo-interactor/KernelUtils.hh @@ -10,7 +10,7 @@ #include "base/ArrayUtils.hh" #include "base/Macros.hh" #include "physics/base/ParticleTrackView.hh" -#include "physics/grid/PhysicsGridCalculator.hh" +#include "physics/grid/XsCalculator.hh" #include "random/distributions/ExponentialDistribution.hh" namespace demo_interactor @@ -19,12 +19,12 @@ namespace demo_interactor template inline CELER_FUNCTION void -move_to_collision(const celeritas::ParticleTrackView& particle, - const celeritas::PhysicsGridCalculator& calc_xs, - const celeritas::Real3& direction, - celeritas::Real3* position, - celeritas::real_type* time, - Rng& rng) +move_to_collision(const celeritas::ParticleTrackView& particle, + const celeritas::XsCalculator& calc_xs, + const celeritas::Real3& direction, + celeritas::Real3* position, + celeritas::real_type* time, + Rng& rng) { CELER_EXPECT(position && time); using celeritas::real_type; diff --git a/src/physics/grid/ValueGridBuilder.cc b/src/physics/grid/ValueGridBuilder.cc index ba09c43482..8a60294f78 100644 --- a/src/physics/grid/ValueGridBuilder.cc +++ b/src/physics/grid/ValueGridBuilder.cc @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #include "ValueGridBuilder.hh" +#include #include #include "base/SoftEqual.hh" #include "physics/grid/UniformGrid.hh" From a551677a072852b32700ea2e7bd10aeea73d4e51 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 2 Mar 2021 08:39:43 -0500 Subject: [PATCH 08/10] Update test results for changes to range/stopping power --- test/physics/base/Physics.test.cc | 40 +++++++------------------------ 1 file changed, 9 insertions(+), 31 deletions(-) diff --git a/test/physics/base/Physics.test.cc b/test/physics/base/Physics.test.cc index d6f7d09024..f6e0956df5 100644 --- a/test/physics/base/Physics.test.cc +++ b/test/physics/base/Physics.test.cc @@ -424,37 +424,15 @@ TEST_F(PHYS_DEVICE_TEST, all) phys_cuda_test(inp); std::vector step_host(step.size()); step.copy_to_host(make_span(step_host)); - // clang-format on - const double expected_step_host[] = {1666.666666667, - 0.002, - 0.003, - 1666.666666667, - 0.002, - 0.003, - 1666.666666667, - 0.556, - 0.7573333333333, - 1666.666666667, - 8.1598, - 8.1598, - inf, - 8.1598, - 8.1598, - 1.666666666667, - 0.002, - 0.003, - 1.666666666667, - 0.002, - 0.003, - 1.666666666667, - 0.5555555555556, - 0.5555555555556, - 1.666666666667, - 1.25, - 1.25, - inf, - 8.1598, - 8.1598}; // clang-format off + const double expected_step_host[] = { + 1666.666666667, 0.00025, 0.00025, 1666.666666667, 0.0025, + 0.0025, 1666.666666667, 0.6568, 0.6568, 1666.666666667, + 5.15968, 5.15968, inf, 5.15968, 5.15968, + 1.666666666667, 2.5e-07, 2.5e-07, 1.666666666667, 2.5e-06, + 2.5e-06, 1.666666666667, 0.0025, 0.0025, 1.666666666667, + 0.025, 0.025, inf, 0.025, 0.025}; + // clang-format on + PRINT_EXPECTED(step_host); EXPECT_VEC_SOFT_EQ(expected_step_host, step_host); } From 6c1e091672d38237056db641edea98a4a1c3a482 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Mar 2021 13:01:10 -0500 Subject: [PATCH 09/10] Make conversion assigment explicit for Collection. This should prevent accidental construction-from-temporaries like the one causing test failures in ValueGridBuilder. --- scripts/docker/ci/launch-testing.sh | 9 +++++---- src/base/Collection.hh | 16 ++++++++-------- src/base/Collection.i.hh | 12 ++++++------ test/physics/em/LivermorePE.test.cc | 4 +++- test/physics/grid/NonuniformGrid.test.cc | 9 +++++++-- test/physics/grid/ValueGridBuilder.test.cc | 10 ++++++---- 6 files changed, 35 insertions(+), 25 deletions(-) diff --git a/scripts/docker/ci/launch-testing.sh b/scripts/docker/ci/launch-testing.sh index 9d550210d9..d22a4675c0 100755 --- a/scripts/docker/ci/launch-testing.sh +++ b/scripts/docker/ci/launch-testing.sh @@ -1,11 +1,11 @@ -#!/bin/sh -ex +#!/bin/sh -e if [ -z "$1" ]; then echo "Usage: $0 pull-request-id" 1>&2 exit 2; fi -CONTAINER=$(docker run -t -d celeritas/ci-cuda11:latest) +CONTAINER=$(docker run -t -d celeritas/ci-cuda11:2021-01-27) echo "Launched container: ${CONTAINER}" docker exec -i $CONTAINER bash -l < - inline Collection(const Collection& other); + explicit inline Collection(const Collection& other); // Construct from another collection (mutable) template - inline Collection(Collection& other); + explicit inline Collection(Collection& other); //!@{ //! Default assignment @@ -133,13 +133,13 @@ class Collection Collection& operator=(Collection&& other) = default; //!@} - // Assign from another collection in the same memory space - template - inline Collection& operator=(const Collection& other); + // Assign from another collectio + template + inline Collection& operator=(const Collection& other); - // Assign (mutable!) from another collection in the same memory space - template - inline Collection& operator=(Collection& other); + // Assign (mutable!) from another collection + template + inline Collection& operator=(Collection& other); //// ACCESS //// diff --git a/src/base/Collection.i.hh b/src/base/Collection.i.hh index 9b44366b23..730c7a8097 100644 --- a/src/base/Collection.i.hh +++ b/src/base/Collection.i.hh @@ -37,12 +37,12 @@ Collection::Collection(Collection& other) //---------------------------------------------------------------------------// /*! - * Assign from another collection in the same memory space. + * Assign from another collection. */ template -template +template Collection& -Collection::operator=(const Collection& other) +Collection::operator=(const Collection& other) { storage_ = detail::CollectionAssigner()(other.storage_); detail::CollectionStorageValidator()(this->size(), @@ -52,12 +52,12 @@ Collection::operator=(const Collection& other) //---------------------------------------------------------------------------// /*! - * Assign (mutable!) from another collection in the same memory space. + * Assign (mutable!) from another collection. */ template -template +template Collection& -Collection::operator=(Collection& other) +Collection::operator=(Collection& other) { storage_ = detail::CollectionAssigner()(other.storage_); detail::CollectionStorageValidator()(this->size(), diff --git a/test/physics/em/LivermorePE.test.cc b/test/physics/em/LivermorePE.test.cc index f327fa4bd3..3e8453ed39 100644 --- a/test/physics/em/LivermorePE.test.cc +++ b/test/physics/em/LivermorePE.test.cc @@ -505,8 +505,10 @@ TEST_F(LivermorePEInteractorTest, model) EXPECT_EQ(1, grid_storage.size()); // Test cross sections calculated from tables + Collection real_ref{ + real_storage}; celeritas::XsCalculator calc_xs( - grid_storage[ValueGridInserter::XsIndex{0}], real_storage); + grid_storage[ValueGridInserter::XsIndex{0}], real_ref); EXPECT_SOFT_EQ(0.1, calc_xs(MevEnergy{1e-3})); EXPECT_SOFT_EQ(1e-5, calc_xs(MevEnergy{1e2})); EXPECT_SOFT_EQ(1e-9, calc_xs(MevEnergy{1e6})); diff --git a/test/physics/grid/NonuniformGrid.test.cc b/test/physics/grid/NonuniformGrid.test.cc index a73c8a352d..c0abe000a5 100644 --- a/test/physics/grid/NonuniformGrid.test.cc +++ b/test/physics/grid/NonuniformGrid.test.cc @@ -26,11 +26,16 @@ class NonuniformGridTest : public celeritas::Test { auto build = celeritas::make_builder(&data); irange = build.insert_back({0, 1, 3, 3, 7}); + ref = data; } celeritas::ItemRange irange; celeritas::Collection data; + celeritas::Collection + ref; }; //---------------------------------------------------------------------------// @@ -39,7 +44,7 @@ class NonuniformGridTest : public celeritas::Test TEST_F(NonuniformGridTest, accessors) { - GridT grid(irange, data); + GridT grid(irange, ref); EXPECT_EQ(5, grid.size()); EXPECT_EQ(0, grid.front()); @@ -50,7 +55,7 @@ TEST_F(NonuniformGridTest, accessors) TEST_F(NonuniformGridTest, find) { - GridT grid(irange, data); + GridT grid(irange, ref); #if CELERITAS_DEBUG EXPECT_THROW(grid.find(-1), celeritas::DebugError); diff --git a/test/physics/grid/ValueGridBuilder.test.cc b/test/physics/grid/ValueGridBuilder.test.cc index d5b152b39e..f96992ae5e 100644 --- a/test/physics/grid/ValueGridBuilder.test.cc +++ b/test/physics/grid/ValueGridBuilder.test.cc @@ -40,9 +40,11 @@ class ValueGridBuilderTest : public celeritas::Test { b->build(insert); } + real_ref = real_storage; } Collection real_storage; + Collection real_ref; Collection grid_storage; }; @@ -75,13 +77,13 @@ TEST_F(ValueGridBuilderTest, xs_grid) // Test results using the physics calculator ASSERT_EQ(2, grid_storage.size()); { - XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_ref); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); } { - XsCalculator calc_xs(grid_storage[XsIndex{1}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{1}], real_ref); EXPECT_SOFT_EQ(10., calc_xs(Energy{1e-3})); EXPECT_SOFT_EQ(1., calc_xs(Energy{1e-2})); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e-1})); @@ -106,7 +108,7 @@ TEST_F(ValueGridBuilderTest, log_grid) // Test results using the physics calculator ASSERT_EQ(1, grid_storage.size()); { - XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_ref); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); @@ -133,7 +135,7 @@ TEST_F(ValueGridBuilderTest, DISABLED_generic_grid) // Test results using the physics calculator ASSERT_EQ(2, grid_storage.size()); { - XsCalculator calc_xs(grid_storage[XsIndex{0}], real_storage); + XsCalculator calc_xs(grid_storage[XsIndex{0}], real_ref); EXPECT_SOFT_EQ(0.1, calc_xs(Energy{1e1})); EXPECT_SOFT_EQ(0.2, calc_xs(Energy{1e2})); EXPECT_SOFT_EQ(0.3, calc_xs(Energy{1e3})); From a3eb048b62ba34cdcbc64dea48a76dcdaac8b62a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 6 Mar 2021 20:13:31 -0500 Subject: [PATCH 10/10] Address reviewer feedback --- src/base/Collection.i.hh | 21 ++++++++------------- src/physics/base/PhysicsStepUtils.i.hh | 22 ++++++++++++++-------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/src/base/Collection.i.hh b/src/base/Collection.i.hh index 730c7a8097..876486ac35 100644 --- a/src/base/Collection.i.hh +++ b/src/base/Collection.i.hh @@ -10,8 +10,14 @@ namespace celeritas { //---------------------------------------------------------------------------// +//!@{ /*! - * Construct from another collection. + * Construct or assign from another collection. + * + * These are generally used to create "references" to "values" (same memory + * space) but can also be used to copy from device to host. The \c + * detail::CollectionAssigner class statically checks for allowable + * transformations and memory moves. */ template template @@ -22,10 +28,6 @@ Collection::Collection(const Collection& other) other.storage().size()); } -//---------------------------------------------------------------------------// -/*! - * Construct from another collection (mutable). - */ template template Collection::Collection(Collection& other) @@ -35,10 +37,6 @@ Collection::Collection(Collection& other) other.storage().size()); } -//---------------------------------------------------------------------------// -/*! - * Assign from another collection. - */ template template Collection& @@ -50,10 +48,6 @@ Collection::operator=(const Collection& other) return *this; } -//---------------------------------------------------------------------------// -/*! - * Assign (mutable!) from another collection. - */ template template Collection& @@ -64,6 +58,7 @@ Collection::operator=(Collection& other) other.storage().size()); return *this; } +//!@} //---------------------------------------------------------------------------// /*! diff --git a/src/physics/base/PhysicsStepUtils.i.hh b/src/physics/base/PhysicsStepUtils.i.hh index 366489cacb..52acb5ffdf 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -85,16 +85,21 @@ calc_tabulated_physics_step(const MaterialTrackView& material, * * See section 7.2.4 Run Time Energy Loss Computation of the Geant4 physics * manual. See also the longer discussions in section 8 of PHYS010 of the - * Geant3 manual (1993). - * - * Energy loss rate (stopping power) is a function of differential cross - * section: the integral of low-energy secondaries (below \c T) produced as a - * function of energy: + * Geant3 manual (1993). Stopping power is an integral over low-exiting-energy + * secondaries. Above some threshold energy \em T_c we treat exiting + * secondaries discretely; below it, we lump them into this continuous loss + * term that varies based on the energy, the atomic number density, and the + * element number: * \f[ - * \frac{dE}{dx} = N_Z \int_0^T \frac{d \sigma_Z(E, T)}{dT} T dT + * \frac{dE}{dx} = N_Z \int_0^{T_c} \frac{d \sigma_Z(E, T)}{dT} T dT * \f] + * Here, the cross section is a function of the primary's energy \em E and the + * exiting secondary energy \em T. * - * The stopping range \em R due to these low-energy processes is: + * The stopping range \em R due to these low-energy processes is an integral + * over the inverse of the stopping power: basically the distance that will + * take a particle (without discrete processes at play) from its current energy + * to an energy of zero. * \f[ * R = \int_0 ^{E_0} - \frac{dx}{dE} dE . * \f] @@ -103,7 +108,8 @@ calc_tabulated_physics_step(const MaterialTrackView& material, * over all processes, rather than the range as a result from integrating all * energy loss processes over the allowed energy range. This is usuallly not * a problem in practice because the range will get automatically decreased by - * \c range_to_step . + * \c range_to_step , and the above range calculation neglects energy loss by + * discrete processes. * * Geant4's stepping algorithm independently stores the range for each process, * then (looping externally over all processes) calculates energy loss, checks