diff --git a/app/demo-interactor/HostKNDemoRunner.cc b/app/demo-interactor/HostKNDemoRunner.cc index 9f21a019ff..e1a5184d25 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/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/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..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,14 +37,10 @@ Collection::Collection(Collection& other) other.storage().size()); } -//---------------------------------------------------------------------------// -/*! - * Assign from another collection in the same memory space. - */ template -template +template Collection& -Collection::operator=(const Collection& other) +Collection::operator=(const Collection& other) { storage_ = detail::CollectionAssigner()(other.storage_); detail::CollectionStorageValidator()(this->size(), @@ -50,20 +48,17 @@ Collection::operator=(const Collection& other) return *this; } -//---------------------------------------------------------------------------// -/*! - * Assign (mutable!) from another collection in the same memory space. - */ template -template +template Collection& -Collection::operator=(Collection& other) +Collection::operator=(Collection& other) { storage_ = detail::CollectionAssigner()(other.storage_); detail::CollectionStorageValidator()(this->size(), other.storage().size()); return *this; } +//!@} //---------------------------------------------------------------------------// /*! 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 7feb2b7aff..617d9a6917 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,10 @@ 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 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 0b71ea2230..d8be66a5e8 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 006ac5d562..7dd3d26029 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..bfc6f5c812 100644 --- a/src/physics/base/PhysicsStepUtils.hh +++ b/src/physics/base/PhysicsStepUtils.hh @@ -22,10 +22,10 @@ calc_tabulated_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 2368917a83..52acb5ffdf 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -10,15 +10,19 @@ #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 +inline CELER_FUNCTION real_type calc_tabulated_physics_step(const MaterialTrackView& material, const ParticleTrackView& particle, PhysicsTrackView& physics) @@ -32,10 +36,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::size_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 +53,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 +61,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) { @@ -80,36 +81,117 @@ calc_tabulated_physics_step(const MaterialTrackView& material, //---------------------------------------------------------------------------// /*! - * 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). 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_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 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] + * + * 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. This is usuallly not + * a problem in practice because the range will get automatically decreased by + * \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 + * 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 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. */ -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::size_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. 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); + eloss += (pre_step_energy.value() + - calc_energy(remaining_range).value()); + } + } + CELER_ASSERT(eloss > 0); + CELER_ASSERT(eloss < pre_step_energy.value()); + } - 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 e183130508..e0d3dbf864 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,10 @@ 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 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) @@ -106,14 +105,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 3db20d0704..a0db161ac5 100644 --- a/src/physics/base/PhysicsTrackView.i.hh +++ b/src/physics/base/PhysicsTrackView.i.hh @@ -172,7 +172,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 +262,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 +299,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/em/AtomicRelaxationParams.cc b/src/physics/em/AtomicRelaxationParams.cc index dfad08b2db..52e8a3d44b 100644 --- a/src/physics/em/AtomicRelaxationParams.cc +++ b/src/physics/em/AtomicRelaxationParams.cc @@ -108,9 +108,9 @@ AtomicRelaxationParams::AtomicRelaxationParams(const Input& inp) AtomicRelaxParamsPointers AtomicRelaxationParams::host_pointers() const { AtomicRelaxParamsPointers result; - result.elements = make_span(host_elements_); - result.electron_id = electron_id_; - result.gamma_id = gamma_id_; + result.elements = make_span(host_elements_); + result.electron_id = electron_id_; + result.gamma_id = gamma_id_; result.electron_cut = electron_cut_; result.gamma_cut = gamma_cut_; @@ -127,9 +127,9 @@ AtomicRelaxParamsPointers AtomicRelaxationParams::device_pointers() const CELER_EXPECT(celeritas::device()); AtomicRelaxParamsPointers result; - result.elements = device_elements_.device_pointers(); - result.electron_id = electron_id_; - result.gamma_id = gamma_id_; + result.elements = device_elements_.device_pointers(); + result.electron_id = electron_id_; + result.gamma_id = gamma_id_; result.electron_cut = electron_cut_; result.gamma_cut = gamma_cut_; 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..14acc6733c --- /dev/null +++ b/src/physics/grid/InverseRangeCalculator.hh @@ -0,0 +1,62 @@ +//----------------------------------*-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/Collection.hh" +#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 + = Collection; + //!@} + + 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/NonuniformGrid.hh b/src/physics/grid/NonuniformGrid.hh new file mode 100644 index 0000000000..d114e9afd5 --- /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/Collection.hh" +#include "base/Macros.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 + = Collection; + //!@} + + public: + // Construct with data + explicit inline CELER_FUNCTION + NonuniformGrid(const ItemRange& 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..b16d501404 --- /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 ItemRange& 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/RangeCalculator.hh b/src/physics/grid/RangeCalculator.hh new file mode 100644 index 0000000000..068ddf4e2f --- /dev/null +++ b/src/physics/grid/RangeCalculator.hh @@ -0,0 +1,58 @@ +//----------------------------------*-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/Collection.hh" +#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 + = Collection; + //!@} + + 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/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/src/physics/grid/ValueGridBuilder.cc b/src/physics/grid/ValueGridBuilder.cc index 4ae7c3e9e0..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" @@ -159,12 +160,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 +177,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 3d30542826..4145efc606 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 = ItemId; + using VecReal = std::vector; + using SpanConstReal = Span; + using Id = ItemId; //!@} 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/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 950b07195d..ef95c7d24e 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: //!@{ @@ -43,17 +43,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 8d3a41a53d..617aadf281 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,13 @@ 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/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) +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/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 3c557ec06f..b37f3bb1df 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/CollectionStateStore.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); } @@ -423,37 +425,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); } diff --git a/test/physics/base/Physics.test.hh b/test/physics/base/Physics.test.hh index 6c5d63f9e0..a033463a2b 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 @@ -64,7 +66,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); } @@ -87,7 +89,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 5aafe37bc5..36eaba2e6f 100644 --- a/test/physics/base/PhysicsStepUtils.test.cc +++ b/test/physics/base/PhysicsStepUtils.test.cc @@ -33,6 +33,18 @@ class PhysicsStepUtilsTest : public PhysicsTestBase using PhysicsStateStore = CollectionStateStore; + 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(); @@ -43,6 +55,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; @@ -58,82 +94,102 @@ TEST_F(PhysicsStepUtilsTest, calc_tabulated_physics_step) 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}); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); phys.interaction_mfp(1); real_type step = celeritas::calc_tabulated_physics_step(material, particle, phys); EXPECT_SOFT_EQ(1. / 3.e-4, step); } { - 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); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{1}, &particle, "celeriton", MevEnergy{10}); + phys.interaction_mfp(1e-4); real_type step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1.e-2 / 9.e-3, step); + 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(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); + EXPECT_SOFT_EQ(0.6568, step); } { - 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); + PhysicsTrackView phys = this->init_track( + &material, MaterialId{1}, &particle, "celeriton", MevEnergy{1e-2}); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(2.5e-3, step); + } + { + PhysicsTrackView phys = this->init_track(&material, + MaterialId{2}, + &particle, + "anti-celeriton", + MevEnergy{1e-2}); + phys.interaction_mfp(1e-6); real_type step = celeritas::calc_tabulated_physics_step(material, particle, phys); - EXPECT_SOFT_EQ(1.e-2 / 9.e-1, step); + 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(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); + EXPECT_SOFT_EQ(2.5e-5, step); + } + { + PhysicsTrackView phys = this->init_track(&material, + MaterialId{2}, + &particle, + "anti-celeriton", + MevEnergy{10}); + real_type step + = celeritas::calc_tabulated_physics_step(material, particle, phys); + EXPECT_SOFT_EQ(2.5e-2, step); } } -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}); + + { + // Long step, but gamma means no energy loss + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "gamma", MevEnergy{1}); + EXPECT_SOFT_EQ( + 0, celeritas::calc_energy_loss(particle, phys, 1e4).value()); + } + { + PhysicsTrackView phys = this->init_track( + &material, MaterialId{0}, &particle, "celeriton", MevEnergy{10}); + const real_type eloss_rate = 0.2 + 0.4; + + // Tiny step: should still be linear loss (single process) + 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 + 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()); + } +} 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_; } diff --git a/test/physics/em/LivermorePE.test.cc b/test/physics/em/LivermorePE.test.cc index 4c0824a419..3e8453ed39 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" @@ -157,10 +157,10 @@ class LivermorePEInteractorTest : public celeritas_test::InteractorHostTestBase } protected: - AtomicRelaxationParams::Input relax_inp_; - std::shared_ptr relax_params_; - std::shared_ptr livermore_params_; - celeritas::detail::LivermorePEPointers pointers_; + AtomicRelaxationParams::Input relax_inp_; + std::shared_ptr relax_params_; + std::shared_ptr livermore_params_; + celeritas::detail::LivermorePEPointers pointers_; celeritas_test::HostStackAllocatorStore vacancies_; }; @@ -505,8 +505,10 @@ TEST_F(LivermorePEInteractorTest, model) EXPECT_EQ(1, grid_storage.size()); // Test cross sections calculated from tables - celeritas::PhysicsGridCalculator calc_xs( - grid_storage[ValueGridInserter::XsIndex{0}], real_storage); + Collection real_ref{ + real_storage}; + celeritas::XsCalculator calc_xs( + 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/CalculatorTestBase.cc b/test/physics/grid/CalculatorTestBase.cc new file mode 100644 index 0000000000..09d64c4e5d --- /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/CollectionBuilder.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_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..2d9f9c219f --- /dev/null +++ b/test/physics/grid/CalculatorTestBase.hh @@ -0,0 +1,59 @@ +//----------------------------------*-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/Collection.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::Collection; + 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::Collection + value_storage_; + Pointers value_ref_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas_test 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/NonuniformGrid.test.cc b/test/physics/grid/NonuniformGrid.test.cc new file mode 100644 index 0000000000..c0abe000a5 --- /dev/null +++ b/test/physics/grid/NonuniformGrid.test.cc @@ -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 NonuniformGrid.test.cc +//---------------------------------------------------------------------------// +#include "physics/grid/NonuniformGrid.hh" + +#include "base/Collection.hh" +#include "base/CollectionBuilder.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_builder(&data); + irange = build.insert_back({0, 1, 3, 3, 7}); + ref = data; + } + + celeritas::ItemRange irange; + celeritas::Collection + data; + celeritas::Collection + ref; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(NonuniformGridTest, accessors) +{ + GridT grid(irange, ref); + + 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(irange, ref); + +#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 +} 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})); +} diff --git a/test/physics/grid/ValueGridBuilder.test.cc b/test/physics/grid/ValueGridBuilder.test.cc index 76828bff7d..f96992ae5e 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: @@ -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()); { - PhysicsGridCalculator 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})); } { - PhysicsGridCalculator 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()); { - PhysicsGridCalculator 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()); { - PhysicsGridCalculator 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})); diff --git a/test/physics/grid/PhysicsGridCalculator.test.cc b/test/physics/grid/XsCalculator.test.cc similarity index 64% rename from test/physics/grid/PhysicsGridCalculator.test.cc rename to test/physics/grid/XsCalculator.test.cc index fd054147c7..3d3dffa795 100644 --- a/test/physics/grid/PhysicsGridCalculator.test.cc +++ b/test/physics/grid/XsCalculator.test.cc @@ -3,16 +3,17 @@ // 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/CollectionBuilder.hh" +#include "base/Interpolator.hh" #include "base/Range.hh" #include "celeritas_test.hh" +#include "CalculatorTestBase.hh" using namespace celeritas; @@ -20,51 +21,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_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; - Collection 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 +53,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 +78,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 +93,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 +111,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 +128,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); }