Add stopping power calculation

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<Secondary> secondaries(args.max_steps);

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],

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<class Rng>
 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;

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 <<EOF || echo "*BUILD FAILED*"
 set -e
@@ -16,5 +16,6 @@ git checkout mr/$1
 SOURCE_DIR=. BUILD_DIR=build entrypoint-shell ./scripts/build/docker.sh
 EOF
 docker stop --time=0 $CONTAINER
-echo "To resume:  docker exec -i -e 'TERM=xterm-256color' $CONTAINER bash -l"
-echo "To delete:  docker rm -f $CONTAINER"
+echo "To resume: docker start $CONTAINER \\"
+echo "           && docker exec -it -e 'TERM=xterm-256color' $CONTAINER bash -l"
+echo "To delete: docker rm -f $CONTAINER"

src/base/Collection.hh

@@ -121,25 +121,25 @@ class Collection
 
     // Construct from another collection
     template<Ownership W2, MemSpace M2>
-    inline Collection(const Collection<T, W2, M2, I>& other);
+    explicit inline Collection(const Collection<T, W2, M2, I>& other);
 
     // Construct from another collection (mutable)
     template<Ownership W2, MemSpace M2>
-    inline Collection(Collection<T, W2, M2, I>& other);
+    explicit inline Collection(Collection<T, W2, M2, I>& other);
 
     //!@{
     //! Default assignment
     Collection& operator=(const Collection& other) = default;
     Collection& operator=(Collection&& other) = default;
     //!@}
 
-    // Assign from another collection in the same memory space
-    template<Ownership W2>
-    inline Collection& operator=(const Collection<T, W2, M, I>& other);
+    // Assign from another collectio
+    template<Ownership W2, MemSpace M2>
+    inline Collection& operator=(const Collection<T, W2, M2, I>& other);
 
-    // Assign (mutable!) from another collection in the same memory space
-    template<Ownership W2>
-    inline Collection& operator=(Collection<T, W2, M, I>& other);
+    // Assign (mutable!) from another collection
+    template<Ownership W2, MemSpace M2>
+    inline Collection& operator=(Collection<T, W2, M2, I>& other);
 
     //// ACCESS ////
 

src/base/Collection.i.hh

@@ -37,12 +37,12 @@ Collection<T, W, M, I>::Collection(Collection<T, W2, M2, I>& other)
 
 //---------------------------------------------------------------------------//
 /*!
- * Assign from another collection in the same memory space.
+ * Assign from another collection.
  */
 template<class T, Ownership W, MemSpace M, class I>
-template<Ownership W2>
+template<Ownership W2, MemSpace M2>
 Collection<T, W, M, I>&
-Collection<T, W, M, I>::operator=(const Collection<T, W2, M, I>& other)
+Collection<T, W, M, I>::operator=(const Collection<T, W2, M2, I>& other)
 {
     storage_ = detail::CollectionAssigner<W, M>()(other.storage_);
     detail::CollectionStorageValidator<W2>()(this->size(),
@@ -52,12 +52,12 @@ Collection<T, W, M, I>::operator=(const Collection<T, W2, M, I>& other)
 
 //---------------------------------------------------------------------------//
 /*!
- * Assign (mutable!) from another collection in the same memory space.
+ * Assign (mutable!) from another collection.
  */
 template<class T, Ownership W, MemSpace M, class I>
-template<Ownership W2>
+template<Ownership W2, MemSpace M2>
 Collection<T, W, M, I>&
-Collection<T, W, M, I>::operator=(Collection<T, W2, M, I>& other)
+Collection<T, W, M, I>::operator=(Collection<T, W2, M2, I>& other)
 {
     storage_ = detail::CollectionAssigner<W, M>()(other.storage_);
     detail::CollectionStorageValidator<W2>()(this->size(),

src/physics/base/ParticleTrackView.hh

@@ -35,6 +35,7 @@ class ParticleTrackView
         = ParticleParamsData<Ownership::const_reference, MemSpace::native>;
     using ParticleStateRef
         = ParticleStateData<Ownership::reference, MemSpace::native>;
+    using Energy        = units::MevEnergy;
     using Initializer_t = ParticleTrackState;
     //!@}
 
@@ -49,15 +50,15 @@ 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) ////
 
     // Unique particle type identifier
     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;

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;
 }

src/physics/base/PhysicsInterface.hh

@@ -170,15 +170,16 @@ 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 ////
 
     //! True if assigned
     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
 };

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;
 }
 
 //---------------------------------------------------------------------------//

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

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<class Engine>
 inline CELER_FUNCTION ModelId select_model(const ParticleTrackView& particle,