Support native CLHEP unit system (#1085)

.github/workflows/build-full.yml

@@ -13,6 +13,7 @@ jobs:
   gpu:
     name: gpu
     strategy:
+      fail-fast: false
       matrix:
         special: [null]
         geometry: ['orange', 'vecgeom']
@@ -46,6 +47,10 @@ jobs:
           - geometry: 'vecgeom'
             buildtype: 'reldeb'
             image: 'ubuntu-cuda'
+          - special: 'clhep'
+            geometry: 'vecgeom'
+            buildtype: 'debug'
+            image: 'ubuntu-cuda'
     env:
       ASAN_OPTIONS: "detect_leaks=0"
       CELER_TEST_STRICT: 1

CMakeLists.txt

@@ -136,6 +136,15 @@ celeritas_setup_option(CELERITAS_REAL_TYPE float)
 celeritas_define_options(CELERITAS_REAL_TYPE
   "Global runtime precision for real numbers")
 
+if((CELERITAS_CORE_GEO STREQUAL "ORANGE")
+    AND (NOT CELERITAS_UNITS STREQUAL "CGS"))
+  celeritas_error_incompatible_option(
+    "ORANGE currently requires CGS units"
+    CELERITAS_UNITS
+    CGS
+  )
+endif()
+
 #----------------------------------------------------------------------------#
 # CMAKE VERSION CHECKS
 #----------------------------------------------------------------------------#

app/celer-g4/DetectorConstruction.cc

@@ -205,7 +205,7 @@ auto DetectorConstruction::construct_field() const -> FieldData
         if (norm(field_val) > 0)
         {
             CELER_LOG_LOCAL(info)
-                << "Using a uniform field " << field_val << " [tesla]";
+                << "Using a uniform field " << field_val << " [T]";
             g4field = std::make_shared<G4UniformMagField>(
                 convert_to_geant(field_val, CLHEP::tesla));
         }
@@ -239,9 +239,9 @@ void DetectorConstruction::ConstructSDandField()
         auto const& field_options = GlobalSetup::Instance()->GetFieldOptions();
         auto chord_finder = std::make_unique<G4ChordFinder>(
             mag_field_.get(),
-            convert_to_geant(field_options.minimum_step, CLHEP::cm));
+            convert_to_geant(field_options.minimum_step, clhep_length));
         chord_finder->SetDeltaChord(
-            convert_to_geant(field_options.delta_chord, CLHEP::cm));
+            convert_to_geant(field_options.delta_chord, clhep_length));
 
         // Construct the magnetic field
         G4FieldManager* field_manager
@@ -251,7 +251,7 @@ void DetectorConstruction::ConstructSDandField()
         field_manager->SetChordFinder(chord_finder.release());
         field_manager->SetMinimumEpsilonStep(field_options.epsilon_step);
         field_manager->SetDeltaIntersection(
-            convert_to_geant(field_options.delta_intersection, CLHEP::cm));
+            convert_to_geant(field_options.delta_intersection, clhep_length));
     }
 
     auto sd_type = GlobalSetup::Instance()->input().sd_type;

app/celer-g4/PGPrimaryGeneratorAction.cc

@@ -74,7 +74,7 @@ void PGPrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
     {
         gun_.SetParticleDefinition(particle_def_[i % particle_def_.size()]);
         gun_.SetParticlePosition(
-            convert_to_geant(sample_pos_(rng_), CLHEP::cm));
+            convert_to_geant(sample_pos_(rng_), clhep_length));
         gun_.SetParticleMomentumDirection(
             convert_to_geant(sample_dir_(rng_), 1));
         gun_.SetParticleEnergy(

app/celer-g4/test-harness.py

@@ -132,7 +132,7 @@ def strtobool(text):
         except:
             pass
         else:
-            outfilename = f'{run_name}.out.failed.json'
+            outfilename = f'{problem_name}.out.failed.json'
             with open(outfilename, 'w') as f:
                 json.dump(j, f, indent=1)
             print("Failure written to", outfilename, file=stderr)

interface/celeritas.i

@@ -69,6 +69,11 @@
 #include "celeritas/Constants.hh"
 %}
 
+namespace celeritas
+{
+enum class UnitSystem;
+}
+
 %include "celeritas/Units.hh"
 %include "celeritas/Constants.hh"
 
@@ -111,6 +116,8 @@ namespace celeritas
 %include "celeritas/io/ImportPhysicsVector.hh"
 %template(VecImportPhysicsVector) std::vector<celeritas::ImportPhysicsVector>;
 
+%include "celeritas/io/ImportUnits.hh"
+
 %include "celeritas/io/ImportPhysicsTable.hh"
 %template(VecImportPhysicsTable) std::vector<celeritas::ImportPhysicsTable>;
 

scripts/cmake-presets/ci-ubuntu-cuda.json

@@ -73,6 +73,15 @@
         "CELERITAS_USE_OpenMP": {"type": "BOOL", "value": "OFF"}
       }
     },
+    {
+      "name": "debug-vecgeom-clhep",
+      "description": "Build with VecGeom and CLHEP units",
+      "inherits": [".vecgeom", ".debug", "base"],
+      "cacheVariables": {
+        "CELERITAS_USE_CUDA": {"type": "BOOL", "value": "OFF"},
+        "CELERITAS_UNITS": "CLHEP"
+      }
+    },
     {
       "name": "reldeb-vecgeom",
       "description": "Build with RelWithDebInfo, assertions, and VecGeom",

scripts/cmake-presets/goldfinger.json

@@ -58,6 +58,15 @@
         "CELERITAS_REAL_TYPE": "float"
       }
     },
+    {
+      "name": "clhep",
+      "displayName": "With CLHEP units",
+      "inherits": [".base", ".debug", "default"],
+      "cacheVariables": {
+        "CELERITAS_UNITS": "CLHEP",
+        "CELERITAS_USE_VecGeom": {"type": "BOOL",   "value": "ON"}
+      }
+    },
     {
       "name": "vecgeom",
       "displayName": "With vecgeom",

scripts/cmake-presets/wildstyle.json

@@ -72,6 +72,15 @@
         "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"}
       }
     },
+    {
+      "name": "clhep",
+      "displayName": "Debug with VecGeom and clhep units",
+      "inherits": [".debug", ".base"],
+      "cacheVariables": {
+        "CELERITAS_UNITS": "CLHEP",
+        "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"}
+      }
+    },
     {
       "name": "reldeb",
       "displayName": "Everything but vecgeom in release mode with debug symbols and assertions",

src/accel/AlongStepFactory.cc

@@ -11,6 +11,7 @@
 
 #include "corecel/io/Logger.hh"
 #include "corecel/math/ArrayUtils.hh"
+#include "corecel/math/QuantityIO.hh"
 #include "celeritas/em/UrbanMscParams.hh"
 #include "celeritas/ext/Convert.geant.hh"
 #include "celeritas/field/RZMapFieldInput.hh"
@@ -47,13 +48,14 @@ auto UniformAlongStepFactory::operator()(AlongStepFactoryInput const& input) con
 {
     // Get the field strength in tesla (or zero if accessor is undefined)
     auto field_params = get_field_ ? get_field_() : UniformFieldParams{};
-    real_type magnitude_tesla = norm(field_params.field) / units::tesla;
+    auto magnitude
+        = native_value_to<units::FieldTesla>(norm(field_params.field));
 
-    if (magnitude_tesla > 0)
+    if (magnitude > zero_quantity())
     {
         // Create a uniform field
         CELER_LOG(info) << "Creating along-step action with field strength "
-                        << magnitude_tesla << "T";
+                        << magnitude;
         return celeritas::AlongStepUniformMscAction::from_params(
             input.action_id,
             *input.material,

src/accel/ExceptionConverter.cc

@@ -102,7 +102,7 @@ void log_state(Logger::Message& msg,
         msg << "\n- Particle type ID: " << kce.particle();
     }
     msg << "\n- Energy: " << kce.energy() << "\n- Position: " << kce.pos()
-        << " (cm)"
+        << " [" << units::NativeTraits::Length::label() << "]"
         << "\n- Direction: " << kce.dir();
 
     if (core_params && kce.volume())

src/accel/LocalTransporter.cc

@@ -156,9 +156,9 @@ void LocalTransporter::Push(G4Track const& g4track)
                    << g4track.GetDefinition()->GetParticleName()
                    << "' particles");
 
-    track.position = convert_from_geant(g4track.GetPosition(), CLHEP::cm);
+    track.position = convert_from_geant(g4track.GetPosition(), clhep_length);
     track.direction = convert_from_geant(g4track.GetMomentumDirection(), 1);
-    track.time = convert_from_geant(g4track.GetGlobalTime(), CLHEP::s);
+    track.time = convert_from_geant(g4track.GetGlobalTime(), clhep_time);
 
     // TODO: Celeritas track IDs are independent from Geant4 track IDs, since
     // they must be sequential from zero for a given event. We may need to save

src/accel/RZMapMagneticField.hh

@@ -62,8 +62,7 @@ RZMapMagneticField::RZMapMagneticField(SPConstFieldParams params)
 void RZMapMagneticField::GetFieldValue(double const pos[3], double* field) const
 {
     // Calculate the magnetic field value in the native Celeritas unit system
-    Real3 result
-        = calc_field_(convert_from_geant(pos, CLHEP::cm) * units::centimeter);
+    Real3 result = calc_field_(convert_from_geant(pos, clhep_length));
     for (auto i = 0; i < 3; ++i)
     {
         // Return values of the field vector in CLHEP::tesla for Geant4

src/accel/detail/HitProcessor.cc

@@ -196,8 +196,8 @@ void HitProcessor::operator()(DetectorStepOutput const& out) const
             {
                 continue;
             }
-            HP_SET(points[sp]->SetGlobalTime, out.points[sp].time, CLHEP::s);
-            HP_SET(points[sp]->SetPosition, out.points[sp].pos, CLHEP::cm);
+            HP_SET(points[sp]->SetGlobalTime, out.points[sp].time, clhep_time);
+            HP_SET(points[sp]->SetPosition, out.points[sp].pos, clhep_length);
             HP_SET(points[sp]->SetKineticEnergy,
                    out.points[sp].energy,
                    CLHEP::MeV);

src/accel/detail/TouchableUpdater.cc

@@ -57,7 +57,7 @@ bool TouchableUpdater::operator()(Real3 const& pos,
                                   Real3 const& dir,
                                   G4LogicalVolume const* lv)
 {
-    auto g4pos = convert_to_geant(pos, CLHEP::cm);
+    auto g4pos = convert_to_geant(pos, clhep_length);
     auto g4dir = convert_to_geant(dir, 1);
 
     // Locate pre-step point
@@ -74,9 +74,9 @@ bool TouchableUpdater::operator()(Real3 const& pos,
         return true;
     }
 
-    constexpr double g4max_step = convert_to_geant(max_step(), CLHEP::cm);
+    constexpr double g4max_step = convert_to_geant(max_step(), clhep_length);
     constexpr double g4max_quiet_step
-        = convert_to_geant(max_quiet_step(), CLHEP::cm);
+        = convert_to_geant(max_quiet_step(), clhep_length);
     double g4safety{-1};
     double g4step{-1};
 
@@ -157,7 +157,7 @@ bool TouchableUpdater::operator()(Real3 const& pos,
     }
 
     // Reset the position and flip the direction
-    g4pos = convert_to_geant(pos, CLHEP::cm);
+    g4pos = convert_to_geant(pos, clhep_length);
     g4dir *= -1;
     find_next_step();
     if (try_cross_boundary())

src/celeritas/CMakeLists.txt

@@ -46,13 +46,16 @@ list(APPEND SOURCES
   grid/ValueGridInserter.cc
   grid/VectorUtils.cc
   io/AtomicRelaxationReader.cc
+  io/ImportData.cc
   io/ImportMaterial.cc
   io/ImportModel.cc
   io/ImportPhysicsTable.cc
   io/ImportPhysicsVector.cc
   io/ImportProcess.cc
+  io/ImportUnits.cc
   io/LivermorePEReader.cc
   io/SeltzerBergerReader.cc
+  io/detail/ImportDataConverter.cc
   mat/MaterialParams.cc
   mat/MaterialParamsOutput.cc
   mat/detail/Utils.cc
@@ -194,6 +197,7 @@ if(CELERITAS_USE_ROOT)
     "celeritas/io/ImportPhysicsTable.hh"
     "celeritas/io/ImportPhysicsVector.hh"
     "celeritas/io/ImportProcess.hh"
+    "celeritas/io/ImportUnits.hh"
     "celeritas/io/ImportVolume.hh"
     "celeritas/io/EventData.hh"
     NOINSTALL

src/celeritas/Quantities.hh

@@ -37,6 +37,7 @@ using CmLength = Quantity<Centimeter>;
 using InvCmXs = Quantity<UnitInverse<Centimeter>>;
 using InvCcDensity = Quantity<InvCentimeterCubed>;
 using MolCcDensity = Quantity<MolPerCentimeterCubed>;
+using GramCcDensity = Quantity<GramPerCentimeterCubed>;
 using FieldTesla = Quantity<Tesla>;
 //!@}
 //---------------------------------------------------------------------------//

src/celeritas/UnitTypes.hh

@@ -151,6 +151,12 @@ struct MolPerCentimeterCubed : UnitProduct<Mol, InvCentimeterCubed>
     static char const* label() { return "mol/cm^3"; }
 };
 
+//! Mass density
+struct GramPerCentimeterCubed : UnitProduct<Gram, InvCentimeterCubed>
+{
+    static char const* label() { return "g/cm^3"; }
+};
+
 //!@}
 //---------------------------------------------------------------------------//
 //!@{

src/celeritas/Units.hh

@@ -43,6 +43,12 @@ namespace units
  * Additionally:
  * - radians are used for measures of angle (unitless)
  * - steradians are used for measures of solid angle (unitless)
+ *
+ * TODO: if we're serious about supporting single-precision arithmetic, we
+ * should define a helper class that stores the constant as full precision but
+ * when multiplied by a single/double is truncated to that precision.
+ * Otherwise, if \c real_type is single-precision, then we lose accuracy in
+ * places like the GeantImporter where everything is double precision.
  */
 
 #define CELER_ICRT inline constexpr real_type

src/celeritas/em/data/RayleighData.hh

@@ -23,9 +23,9 @@ namespace celeritas
  * \f[
  *  FF(E,cos)^2 = \Sigma_{j} \frac{a_j}{[1 + b_j x]^{n}}
  * \f]
- * where \f$ x= E^{2}(1-cos\theta) \f$ and \em n is the high energy slope of
- * the form factor and \em a and \em b are free parameters to obtain the best
- * fit to the form factor. The unit for the energy (\em E) is in MeV.
+ * where \f$ x = E^{2}(1 - \cos\theta) \f$ and \em n is the high energy slope
+ * of the form factor and \em a and \em b are free parameters to obtain the
+ * best fit to the form factor. The unit for the energy (\em E) is in MeV.
  */
 struct RayleighParameters
 {
@@ -79,8 +79,6 @@ struct RayleighData
     }
 };
 
-using RayleighDeviceRef = DeviceCRef<RayleighData>;
-using RayleighHostRef = HostCRef<RayleighData>;
 using RayleighRef = NativeCRef<RayleighData>;
-
+//---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/em/interactor/RayleighInteractor.hh

@@ -58,13 +58,6 @@ class RayleighInteractor
 
     //// CONSTANTS ////
 
-    //! cm/hc in the MeV energy unit
-    static CELER_CONSTEXPR_FUNCTION real_type hc_factor()
-    {
-        return units::centimeter * native_value_from(units::MevEnergy{1.0})
-               / (constants::c_light * constants::h_planck);
-    }
-
     //! A point where the functional form of the form factor fit changes
     static CELER_CONSTEXPR_FUNCTION real_type fit_slice() { return 0.02; }
 

src/celeritas/em/model/RayleighModel.cc

@@ -153,6 +153,9 @@ void RayleighModel::build_data(HostValue* data, MaterialParams const& materials)
  * Parameters for Z = 0 are dropped as they are zeros and not used.
  * Reshaped as [el][param][3] with params.T.reshape((100, 3, 3)),
  * then updated 'n' with params[:,2,:] -= 1
+ *
+ * The fit data embeds centimeters as a unit system: this is accounted for in
+ * the interactor.
  */
 auto RayleighModel::get_el_parameters(AtomicNumber z) -> ElScatParams const&
 {

src/celeritas/em/msc/detail/MscStepToGeo.hh

@@ -77,8 +77,8 @@ class MscStepToGeo
 
     struct result_type
     {
-        real_type step{};  //!< Geometrical step length
-        real_type alpha{0};  //!< Scaled MFP slope
+        real_type step{};  //!< Geometrical step length [len]
+        real_type alpha{0};  //!< Scaled MFP slope [1/len]
     };
 
   public:

src/celeritas/em/msc/detail/UrbanMscHelper.hh

@@ -30,6 +30,8 @@ namespace detail
 /*!
  * This is a helper class for the UrbanMscStepLimit and UrbanMscScatter.
  *
+ * NOTE: units are "native" units, listed here as CGS.
+ *
  * \todo Refactor to UrbanMscTrackView .
  */
 class UrbanMscHelper

src/celeritas/em/xs/LPMCalculator.hh

@@ -62,7 +62,7 @@ class LPMCalculator
 
     // Current element
     ElementView const& element_;
-    // Electron density of the current material [1/cm^3]
+    // Electron density of the current material [1/len^3]
     real_type const electron_density_;
     // Characteristic energy for the LPM effect for this material [MeV]
     real_type const lpm_energy_;