Add a common macroscopic cross section calculator (celeritas-project#…

…1145)

* Store neutron-nucleus elastic cross sections in units::barn

* Use a common method for calculating macroscopic cross sections with a templated micro scross section calculator and associated updates

* Add MacroXsCalculator

* Add a new quantity BarnXs

* Return units::BarnXs (quantity) from microscopic cross section calculators

* Use value_as in some places

src/celeritas/Quantities.hh

@@ -33,6 +33,7 @@ using AmuMass = Quantity<Amu>;
 //---------------------------------------------------------------------------//
 //!@{
 //! \name Units for manual input and/or test harnesses
+using BarnXs = Quantity<Barn>;
 using CmLength = Quantity<Centimeter>;
 using InvCmXs = Quantity<UnitInverse<Centimeter>>;
 using InvCcDensity = Quantity<InvCentimeterCubed>;

src/celeritas/em/interactor/LivermorePEInteractor.hh

@@ -216,7 +216,7 @@ CELER_FUNCTION Interaction LivermorePEInteractor::operator()(Engine& rng)
     }
     result.secondaries = secondaries;
 
-    CELER_ENSURE(result.energy_deposition.value() >= 0);
+    CELER_ENSURE(result.energy_deposition >= zero_quantity());
     return result;
 }
 
@@ -231,7 +231,9 @@ CELER_FUNCTION SubshellId LivermorePEInteractor::sample_subshell(Engine& rng) co
     auto const& shells = shared_.xs.shells[el.shells];
     size_type shell_id = 0;
 
-    real_type const cutoff = generate_canonical(rng) * calc_micro_xs_(el_id_);
+    using Xs = Quantity<LivermoreSubshell::XsUnits>;
+    real_type const cutoff = generate_canonical(rng)
+                             * value_as<Xs>(calc_micro_xs_(el_id_));
     if (Energy{inc_energy_} < el.thresh_lo)
     {
         // Accumulate discrete PDF for tabulated shell cross sections

src/celeritas/em/xs/LivermorePEMicroXsCalculator.hh

@@ -28,17 +28,18 @@ class LivermorePEMicroXsCalculator
   public:
     //!@{
     //! \name Type aliases
-    using XsUnits = LivermoreSubshell::XsUnits;
+    using ParamsRef = LivermorePERef;
     using Energy = Quantity<LivermoreSubshell::EnergyUnits>;
+    using BarnXs = units::BarnXs;
     //!@}
 
   public:
     // Construct with shared and state data
     inline CELER_FUNCTION
-    LivermorePEMicroXsCalculator(LivermorePERef const& shared, Energy energy);
+    LivermorePEMicroXsCalculator(ParamsRef const& shared, Energy energy);
 
     // Compute cross section
-    inline CELER_FUNCTION real_type operator()(ElementId el_id) const;
+    inline CELER_FUNCTION BarnXs operator()(ElementId el_id) const;
 
   private:
     // Shared constant physics properties
@@ -54,7 +55,7 @@ class LivermorePEMicroXsCalculator
  * Construct with shared and state data.
  */
 CELER_FUNCTION LivermorePEMicroXsCalculator::LivermorePEMicroXsCalculator(
-    LivermorePERef const& shared, Energy energy)
+    ParamsRef const& shared, Energy energy)
     : shared_(shared), inc_energy_(energy.value())
 {
 }
@@ -64,7 +65,7 @@ CELER_FUNCTION LivermorePEMicroXsCalculator::LivermorePEMicroXsCalculator(
  * Compute cross section
  */
 CELER_FUNCTION
-real_type LivermorePEMicroXsCalculator::operator()(ElementId el_id) const
+auto LivermorePEMicroXsCalculator::operator()(ElementId el_id) const -> BarnXs
 {
     CELER_EXPECT(el_id);
     LivermoreElement const& el = shared_.xs.elements[el_id];
@@ -101,7 +102,7 @@ real_type LivermorePEMicroXsCalculator::operator()(ElementId el_id) const
         GenericCalculator calc_xs(el.xs_lo, shared_.xs.reals);
         result = ipow<3>(inv_energy) * calc_xs(energy.value());
     }
-    return result;
+    return BarnXs{result};
 }
 
 //---------------------------------------------------------------------------//

src/celeritas/io/NeutronXsReader.cc

@@ -15,6 +15,7 @@
 #include "corecel/io/Logger.hh"
 #include "corecel/math/Algorithms.hh"
 #include "corecel/sys/Environment.hh"
+#include "celeritas/Quantities.hh"
 #include "celeritas/Units.hh"
 #include "celeritas/io/ImportPhysicsVector.hh"
 
@@ -86,10 +87,12 @@ NeutronXsReader::operator()(AtomicNumber atomic_number) const
         for (size_type i = 0; i < size; ++i)
         {
             CELER_ASSERT(infile);
-            // Convert to the celeritas native length 1/[len^2] from
-            // clhep::mm^2 as stored in G4PARTICLEXS/neutron/el data
+            // Convert to the celeritas units::barn (units::BarnXs.value())
+            // from clhep::mm^2 as stored in G4PARTICLEXS/neutron/el data
             infile >> result.x[i] >> input_xs.value();
-            result.y[i] = native_value_from(input_xs);
+            result.y[i]
+                = native_value_to<units::BarnXs>(native_value_from(input_xs))
+                      .value();
         }
     }
 

src/celeritas/mat/ElementSelector.hh

@@ -104,7 +104,8 @@ CELER_FUNCTION ElementSelector::ElementSelector(MaterialView const& material,
     CELER_EXPECT(storage.size() >= material.num_elements());
     for (auto i : range<size_type>(elements_.size()))
     {
-        real_type const micro_xs = calc_micro_xs(elements_[i].element);
+        real_type const micro_xs
+            = (calc_micro_xs(elements_[i].element)).value();
         CELER_ASSERT(micro_xs >= 0);
         real_type const frac = elements_[i].fraction;
 

src/celeritas/neutron/xs/NeutronElasticMicroXsCalculator.hh

@@ -27,24 +27,24 @@ class NeutronElasticMicroXsCalculator
   public:
     //!@{
     //! \name Type aliases
+    using ParamsRef = NeutronElasticRef;
     using Energy = units::MevEnergy;
-    using XsUnits = units::Native;  // [len^2]
+    using BarnXs = units::BarnXs;
     //!@}
 
   public:
     // Construct with shared and state data
     inline CELER_FUNCTION
-    NeutronElasticMicroXsCalculator(NeutronElasticRef const& shared,
-                                    Energy energy);
+    NeutronElasticMicroXsCalculator(ParamsRef const& shared, Energy energy);
 
     // Compute cross section
-    inline CELER_FUNCTION real_type operator()(ElementId el_id) const;
+    inline CELER_FUNCTION BarnXs operator()(ElementId el_id) const;
 
   private:
     // Shared constant physics properties
     NeutronElasticRef const& shared_;
     // Incident neutron energy
-    Energy const inc_energy_;
+    real_type const inc_energy_;
 };
 
 //---------------------------------------------------------------------------//
@@ -54,7 +54,7 @@ class NeutronElasticMicroXsCalculator
  * Construct with shared and state data.
  */
 CELER_FUNCTION NeutronElasticMicroXsCalculator::NeutronElasticMicroXsCalculator(
-    NeutronElasticRef const& shared, Energy energy)
+    ParamsRef const& shared, Energy energy)
     : shared_(shared), inc_energy_(energy.value())
 {
 }
@@ -64,7 +64,8 @@ CELER_FUNCTION NeutronElasticMicroXsCalculator::NeutronElasticMicroXsCalculator(
  * Compute microscopic (element) cross section
  */
 CELER_FUNCTION
-real_type NeutronElasticMicroXsCalculator::operator()(ElementId el_id) const
+auto NeutronElasticMicroXsCalculator::operator()(ElementId el_id) const
+    -> BarnXs
 {
     CELER_EXPECT(el_id < shared_.micro_xs.size());
 
@@ -73,9 +74,9 @@ real_type NeutronElasticMicroXsCalculator::operator()(ElementId el_id) const
 
     // Calculate micro cross section at the given energy
     GenericCalculator calc_xs(grid, shared_.reals);
-    real_type result = calc_xs(inc_energy_.value());
+    real_type result = calc_xs(inc_energy_);
 
-    return result;
+    return BarnXs{result};
 }
 
 //---------------------------------------------------------------------------//

src/celeritas/neutron/xs/NeutronElasticMacroXsCalculator.hh → src/celeritas/phys/MacroXsCalculator.hh

@@ -3,45 +3,48 @@
 // See the top-level COPYRIGHT file for details.
 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
 //---------------------------------------------------------------------------//
-//! \file celeritas/neutron/xs/NeutronElasticMacroXsCalculator.hh
+//! \file celeritas/phys/MacroXsCalculator.hh
 //---------------------------------------------------------------------------//
 #pragma once
 
 #include "corecel/Assert.hh"
 #include "corecel/Macros.hh"
 #include "corecel/Types.hh"
 #include "corecel/cont/Span.hh"
+#include "celeritas/Quantities.hh"
 #include "celeritas/UnitTypes.hh"
 #include "celeritas/mat/MaterialView.hh"
 
-#include "NeutronElasticMicroXsCalculator.hh"
-
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
 /*!
  * Calculates the macroscopic cross section.
+ *
+ * \tparam MicroXsT microscopic (element) cross section calculator
  */
-class NeutronElasticMacroXsCalculator
+template<class MicroXsT>
+class MacroXsCalculator
 {
   public:
     //!@{
     //! \name Type aliases
-    using Energy = units::MevEnergy;
-    using XsUnits = units::Native;  // [1/len]
+    using ParamsRef = typename MicroXsT::ParamsRef;
+    using Energy = typename MicroXsT::Energy;
+    using MicroXs = units::BarnXs;
+    using MacroXsUnits = units::Native;  // [1/len]
     //!@}
 
   public:
-    // Construct with shared data and material
+    // Construct with shared and material
     inline CELER_FUNCTION
-    NeutronElasticMacroXsCalculator(NeutronElasticRef const& shared,
-                                    MaterialView const& material);
+    MacroXsCalculator(ParamsRef const& shared, MaterialView const& material);
 
-    // Compute cross section on the fly at the given energy
+    // Compute the macroscopic cross section on the fly at the given energy
     inline CELER_FUNCTION real_type operator()(Energy energy) const;
 
   private:
-    NeutronElasticRef const& shared_;
+    ParamsRef const& shared_;
     Span<MatElementComponent const> elements_;
     real_type number_density_;
 };
@@ -50,10 +53,12 @@ class NeutronElasticMacroXsCalculator
 // INLINE DEFINITIONS
 //---------------------------------------------------------------------------//
 /*!
- * Construct with shared model and material data.
+ * Construct with shared and material data.
  */
-CELER_FUNCTION NeutronElasticMacroXsCalculator::NeutronElasticMacroXsCalculator(
-    NeutronElasticRef const& shared, MaterialView const& material)
+template<class MicroXsT>
+CELER_FUNCTION
+MacroXsCalculator<MicroXsT>::MacroXsCalculator(ParamsRef const& shared,
+                                               MaterialView const& material)
     : shared_(shared)
     , elements_(material.elements())
     , number_density_(material.number_density())
@@ -63,21 +68,21 @@ CELER_FUNCTION NeutronElasticMacroXsCalculator::NeutronElasticMacroXsCalculator(
 
 //---------------------------------------------------------------------------//
 /*!
- * Compute macroscopic cross section for the neutron elastic on the fly at
- * the given energy.
+ * Compute the macroscopic cross section on the fly at the given energy.
  */
+template<class MicroXsT>
 CELER_FUNCTION real_type
-NeutronElasticMacroXsCalculator::operator()(Energy energy) const
+MacroXsCalculator<MicroXsT>::operator()(Energy energy) const
 {
     real_type result = 0.;
-    NeutronElasticMicroXsCalculator calc_micro_xs(shared_, energy);
+    MicroXsT calc_micro_xs(shared_, energy);
     for (auto const& el_comp : elements_)
     {
-        real_type const micro_xs = calc_micro_xs(el_comp.element);
+        real_type micro_xs = value_as<MicroXs>(calc_micro_xs(el_comp.element));
         CELER_ASSERT(micro_xs >= 0);
         result += micro_xs * el_comp.fraction;
     }
-    result *= XsUnits::value() * number_density_;
+    result = native_value_from(MicroXs{result}) * number_density_;
     CELER_ENSURE(result >= 0);
     return result;
 }

src/celeritas/phys/PhysicsTrackView.hh

@@ -14,12 +14,13 @@
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/em/xs/EPlusGGMacroXsCalculator.hh"
-#include "celeritas/em/xs/LivermorePEMacroXsCalculator.hh"
+#include "celeritas/em/xs/LivermorePEMicroXsCalculator.hh"
 #include "celeritas/grid/GridIdFinder.hh"
 #include "celeritas/grid/XsCalculator.hh"
 #include "celeritas/mat/MaterialView.hh"
 #include "celeritas/mat/TabulatedElementSelector.hh"
-#include "celeritas/neutron/xs/NeutronElasticMacroXsCalculator.hh"
+#include "celeritas/neutron/xs/NeutronElasticMicroXsCalculator.hh"
+#include "celeritas/phys/MacroXsCalculator.hh"
 
 #include "PhysicsData.hh"
 
@@ -407,7 +408,7 @@ CELER_FUNCTION real_type PhysicsTrackView::calc_xs(ParticleProcessId ppid,
         // hardwired processes.
         if (model_id == params_.hardwired.livermore_pe)
         {
-            auto calc_xs = LivermorePEMacroXsCalculator(
+            auto calc_xs = MacroXsCalculator<LivermorePEMicroXsCalculator>(
                 params_.hardwired.livermore_pe_data, material);
             result = calc_xs(energy);
         }
@@ -419,7 +420,7 @@ CELER_FUNCTION real_type PhysicsTrackView::calc_xs(ParticleProcessId ppid,
         }
         else if (model_id == params_.hardwired.chips)
         {
-            auto calc_xs = NeutronElasticMacroXsCalculator(
+            auto calc_xs = MacroXsCalculator<NeutronElasticMicroXsCalculator>(
                 params_.hardwired.chips_data, material);
             result = calc_xs(energy);
         }