Define unit traits and support arbitrary RZ field map units (celerita…

…s-project#1088)

* Add UnitSystem enum and traits visitor
* Fix RZ field map input
* Fix constexpr warnings
* Apply fixes

src/celeritas/Types.cc

@@ -8,9 +8,40 @@
 #include "Types.hh"
 
 #include "corecel/io/EnumStringMapper.hh"
+#include "corecel/io/StringEnumMapper.hh"
+
+#include "UnitTypes.hh"
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Get a string corresponding to a unit system.
+ */
+char const* to_cstring(UnitSystem value)
+{
+    static_assert(static_cast<int>(UnitSystem::cgs) == CELERITAS_UNITS_CGS);
+    static_assert(static_cast<int>(UnitSystem::si) == CELERITAS_UNITS_SI);
+    static_assert(static_cast<int>(UnitSystem::clhep) == CELERITAS_UNITS_CLHEP);
+    static_assert(static_cast<int>(UnitSystem::native) == CELERITAS_UNITS);
+
+    static EnumStringMapper<UnitSystem> const to_cstring_impl{
+        "none", "cgs", "si", "clhep"};
+    return to_cstring_impl(value);
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Get a unit system corresponding to a string value.
+ */
+UnitSystem to_unit_system(std::string const& s)
+{
+    static auto const from_string
+        = StringEnumMapper<UnitSystem>::from_cstring_func(to_cstring,
+                                                          "unit system");
+    return from_string(s);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Get a string corresponding to a state of matter.

src/celeritas/Types.hh

@@ -10,6 +10,7 @@
 
 #include <cstdint>
 
+#include "celeritas_config.h"
 // IWYU pragma: begin_exports
 #include "corecel/Assert.hh"
 #include "corecel/Macros.hh"
@@ -77,6 +78,18 @@ using SubshellId = OpaqueId<struct Subshell_>;
 
 //---------------------------------------------------------------------------//
 // ENUMERATIONS
+//---------------------------------------------------------------------------//
+//! Unit system used by Celeritas
+enum class UnitSystem
+{
+    none,  //!< Invalid unit system
+    cgs,  //!< Gaussian CGS
+    si,  //!< International System
+    clhep,  //!< Geant4 native
+    size_,
+    native = CELERITAS_UNITS,  //!< Compile time selected system
+};
+
 //---------------------------------------------------------------------------//
 //! Interpolation type
 enum class Interp
@@ -169,6 +182,12 @@ struct StepLimit
 // HELPER FUNCTIONS (HOST)
 //---------------------------------------------------------------------------//
 
+// Get a string corresponding to a unit system
+char const* to_cstring(UnitSystem);
+
+// Get a unit system corresponding to a string
+UnitSystem to_unit_system(std::string const& s);
+
 // Get a string corresponding to a material state
 char const* to_cstring(MatterState);
 

src/celeritas/UnitTypes.hh

@@ -8,10 +8,13 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include <utility>
+
 #include "celeritas_config.h"
 #include "corecel/math/UnitUtils.hh"
 
 #include "Constants.hh"
+#include "Types.hh"
 #include "Units.hh"
 
 namespace celeritas
@@ -251,43 +254,75 @@ struct LogMev
 //---------------------------------------------------------------------------//
 
 //! Traits class for units
-template<int>
-struct UnitTraits;
+template<UnitSystem>
+struct UnitSystemTraits;
 
 //! CGS unit traits
 template<>
-struct UnitTraits<CELERITAS_UNITS_CGS>
+struct UnitSystemTraits<UnitSystem::cgs>
 {
     using Length = Centimeter;
     using Mass = Gram;
     using Time = Second;
     using BField = Gauss;  //!< Magnetic flux density
+
+    static char const* label() { return "cgs"; }
 };
 
 //! SI unit traits
 template<>
-struct UnitTraits<CELERITAS_UNITS_SI>
+struct UnitSystemTraits<UnitSystem::si>
 {
     using Length = Meter;
     using Mass = Kilogram;
     using Time = Second;
     using BField = Tesla;
+
+    static char const* label() { return "si"; }
 };
 
 //! CLHEP unit traits
 template<>
-struct UnitTraits<CELERITAS_UNITS_CLHEP>
+struct UnitSystemTraits<UnitSystem::clhep>
 {
     using Length = Millimeter;
     using Mass = ClhepUnitMass;
     using Time = Nanosecond;
     using BField = ClhepUnitBField;
+
+    static char const* label() { return "clhep"; }
 };
 
-using CgsTraits = UnitTraits<CELERITAS_UNITS_CGS>;
-using SiTraits = UnitTraits<CELERITAS_UNITS_SI>;
-using ClhepTraits = UnitTraits<CELERITAS_UNITS_CLHEP>;
-using NativeTraits = UnitTraits<CELERITAS_UNITS>;
+using CgsTraits = UnitSystemTraits<UnitSystem::cgs>;
+using SiTraits = UnitSystemTraits<UnitSystem::si>;
+using ClhepTraits = UnitSystemTraits<UnitSystem::clhep>;
+using NativeTraits = UnitSystemTraits<UnitSystem::native>;
+
+//---------------------------------------------------------------------------//
+/*!
+ * Expand a macro to a switch statement over all possible unit system types.
+ *
+ * This is *not* a \c CELER_FUNCTION because unit conversion should be done
+ * only during preprocessing on the CPU.
+ */
+template<class F>
+constexpr decltype(auto) visit_unit_system(F&& func, UnitSystem sys)
+{
+#define CELER_US_VISIT_CASE(TYPE) \
+    case UnitSystem::TYPE:        \
+        return std::forward<F>(func)(UnitSystemTraits<UnitSystem::TYPE>{});
+
+    switch (sys)
+    {
+        CELER_US_VISIT_CASE(cgs);
+        CELER_US_VISIT_CASE(si);
+        CELER_US_VISIT_CASE(clhep);
+        default:
+            CELER_ASSERT_UNREACHABLE();
+    }
+
+#undef CELER_US_VISIT_CASE
+}
 
 //---------------------------------------------------------------------------//
 }  // namespace units

src/celeritas/field/RZMapFieldInput.hh

@@ -25,8 +25,9 @@ namespace celeritas
  * The magnetic field is discretized at nodes on an R-Z grid, and each point
  * the field vector is approximated by a 2-D vector in R-Z. The input units of
  * this field are in *NATIVE UNITS* (cm/gauss when CGS). An optional \c _units
- * field in the input can specify whether the input is in tesla or native
- * units, with allowable values of "T", "tesla", "gauss", or "native".
+ * field in the input can specify whether the input is in SI or CGS
+ * units, with allowable values of "si", "cgs", or "clhep". The native CLHEP
+ * unit strength is 1000*tesla.
  *
  * The field values are all indexed with R having stride 1: [Z][R]
  */
@@ -35,8 +36,8 @@ struct RZMapFieldInput
     unsigned int num_grid_z{};
     unsigned int num_grid_r{};
     double min_z{};  //!< Lower z coordinate [len]
-    double min_r{};  //!< Lower r coordinate [len]
     double max_z{};  //!< Last z coordinate [len]
+    double min_r{};  //!< Lower r coordinate [len]
     double max_r{};  //!< Last r coordinate [len]
     std::vector<double> field_z;  //!< Flattened Z field component [bfield]
     std::vector<double> field_r;  //!< Flattened R field component [bfield]

src/celeritas/field/RZMapFieldInputIO.json.cc

@@ -13,6 +13,7 @@
 #include <vector>
 
 #include "corecel/cont/Range.hh"
+#include "corecel/io/Logger.hh"
 #include "celeritas/Quantities.hh"
 
 #include "FieldDriverOptionsIO.json.hh"
@@ -27,6 +28,8 @@ namespace celeritas
 void from_json(nlohmann::json const& j, RZMapFieldInput& inp)
 {
 #define RZFI_LOAD(NAME) j.at(#NAME).get_to(inp.NAME)
+    using namespace celeritas::units;
+
     RZFI_LOAD(num_grid_z);
     RZFI_LOAD(num_grid_r);
     RZFI_LOAD(min_z);
@@ -40,38 +43,96 @@ void from_json(nlohmann::json const& j, RZMapFieldInput& inp)
         RZFI_LOAD(driver_options);
     }
 
-    // Convert from tesla if explicitly specified *OR* if no _units given
-    bool convert_from_tesla{true};
+    // Convert unit systems based on input
+    UnitSystem length_units{UnitSystem::cgs};  // cm
+    UnitSystem field_units{UnitSystem::si};  // tesla
     if (auto iter = j.find("_units"); iter != j.end())
     {
-        auto units = iter->get<std::string>();
-        if (units == "tesla" || units == "T")
+        auto const& ustr = iter->get<std::string>();
+        if (ustr == "tesla" || ustr == "T")
         {
-            convert_from_tesla = true;
+            CELER_LOG(warning)
+                << "Deprecated RZ field input units '" << ustr
+                << "': use SI units for length (m) and field (T) "
+                   "and set units to 'si'";
+            field_units = UnitSystem::si;
         }
-        else if (units == "native" || units == "gauss")
+        else if (ustr == "gauss" || ustr == Gauss::label() || ustr == "native")
         {
-            convert_from_tesla = false;
+            // TODO: Remove in 1.0
+            CELER_LOG(warning) << "Deprecated RZ field input units '" << ustr
+                               << "': replace with 'cgs' (Gauss + cm)";
+            field_units = UnitSystem::cgs;
         }
         else
         {
-            CELER_VALIDATE(false,
-                           << "unrecognized value '" << units
-                           << "' for \"_units\" field");
+            try
+            {
+                // Input should be si/cgs/clhep
+                length_units = to_unit_system(ustr);
+                field_units = length_units;
+            }
+            catch (RuntimeError const& e)
+            {
+                CELER_VALIDATE(false,
+                               << "unrecognized value '" << ustr
+                               << "' for \"_units\" field: " << e.what());
+            }
         }
     }
-    if (convert_from_tesla)
+    else
     {
-        using FieldTesla = Quantity<units::Tesla, double>;
+        auto msg = CELER_LOG(warning);
+        msg << "No units given in RZ field input: assuming CGS for length "
+               "(cm) and SI for strength (T)";
+    }
+
+    if (field_units != UnitSystem::native)
+    {
+        CELER_LOG(info) << "Converting magnetic field input strength from "
+                        << to_cstring(field_units) << " to ["
+                        << NativeTraits::BField::label() << "]";
+
+        double field_scale = visit_unit_system(
+            [](auto traits) {
+                using Unit = typename decltype(traits)::BField;
+                return native_value_from(Quantity<Unit, double>{1});
+            },
+            field_units);
+
+        CELER_LOG(debug) << "Scaling input magnetic field by " << field_scale;
+
         // Convert units from JSON tesla to input native
         for (auto* f : {&inp.field_z, &inp.field_r})
         {
             for (double& v : *f)
             {
-                v = native_value_from(FieldTesla(v));
+                v *= field_scale;
             }
         }
     }
+
+    if (length_units != UnitSystem::native)
+    {
+        CELER_LOG(info) << "Converting magnetic field input positions from "
+                        << to_cstring(length_units) << " to ["
+                        << NativeTraits::Length::label() << "]";
+
+        double length_scale = visit_unit_system(
+            [](auto traits) {
+                using Unit = typename decltype(traits)::Length;
+                return native_value_from(Quantity<Unit, double>{1});
+            },
+            length_units);
+
+        CELER_LOG(debug) << "Scaling input lengths by " << length_scale;
+
+        // Convert units from JSON tesla to input native
+        for (auto* v : {&inp.min_z, &inp.max_z, &inp.min_r, &inp.max_r})
+        {
+            *v *= length_scale;
+        }
+    }
 #undef RZFI_LOAD
 }
 
@@ -83,7 +144,7 @@ void to_json(nlohmann::json& j, RZMapFieldInput const& inp)
 {
 #define RZFI_KEY_VALUE(NAME) {#NAME, inp.NAME}
     j = {
-        {"_units", "gauss"},
+        {"_units", units::NativeTraits::label()},
         RZFI_KEY_VALUE(num_grid_z),
         RZFI_KEY_VALUE(num_grid_r),
         RZFI_KEY_VALUE(min_z),

src/celeritas/field/RZMapFieldParams.hh

@@ -22,6 +22,8 @@ struct RZMapFieldInput;
 //---------------------------------------------------------------------------//
 /*!
  * Set up a 2D RZMapFieldParams.
+ *
+ * The input values should be converted to the native unit system.
  */
 class RZMapFieldParams final : public ParamsDataInterface<RZMapFieldParamsData>
 {