diff --git a/src/celeritas/Types.cc b/src/celeritas/Types.cc index 8e0e7b5971..cc2682a335 100644 --- a/src/celeritas/Types.cc +++ b/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(UnitSystem::cgs) == CELERITAS_UNITS_CGS); + static_assert(static_cast(UnitSystem::si) == CELERITAS_UNITS_SI); + static_assert(static_cast(UnitSystem::clhep) == CELERITAS_UNITS_CLHEP); + static_assert(static_cast(UnitSystem::native) == CELERITAS_UNITS); + + static EnumStringMapper 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::from_cstring_func(to_cstring, + "unit system"); + return from_string(s); +} + //---------------------------------------------------------------------------// /*! * Get a string corresponding to a state of matter. diff --git a/src/celeritas/Types.hh b/src/celeritas/Types.hh index 367d3e9ab8..05e69a1b85 100644 --- a/src/celeritas/Types.hh +++ b/src/celeritas/Types.hh @@ -10,6 +10,7 @@ #include +#include "celeritas_config.h" // IWYU pragma: begin_exports #include "corecel/Assert.hh" #include "corecel/Macros.hh" @@ -77,6 +78,18 @@ using SubshellId = OpaqueId; //---------------------------------------------------------------------------// // 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); diff --git a/src/celeritas/UnitTypes.hh b/src/celeritas/UnitTypes.hh index acc5624d1a..2d6fc27be9 100644 --- a/src/celeritas/UnitTypes.hh +++ b/src/celeritas/UnitTypes.hh @@ -8,10 +8,13 @@ //---------------------------------------------------------------------------// #pragma once +#include + #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 -struct UnitTraits; +template +struct UnitSystemTraits; //! CGS unit traits template<> -struct UnitTraits +struct UnitSystemTraits { 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 +struct UnitSystemTraits { using Length = Meter; using Mass = Kilogram; using Time = Second; using BField = Tesla; + + static char const* label() { return "si"; } }; //! CLHEP unit traits template<> -struct UnitTraits +struct UnitSystemTraits { using Length = Millimeter; using Mass = ClhepUnitMass; using Time = Nanosecond; using BField = ClhepUnitBField; + + static char const* label() { return "clhep"; } }; -using CgsTraits = UnitTraits; -using SiTraits = UnitTraits; -using ClhepTraits = UnitTraits; -using NativeTraits = UnitTraits; +using CgsTraits = UnitSystemTraits; +using SiTraits = UnitSystemTraits; +using ClhepTraits = UnitSystemTraits; +using NativeTraits = UnitSystemTraits; + +//---------------------------------------------------------------------------// +/*! + * 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 +constexpr decltype(auto) visit_unit_system(F&& func, UnitSystem sys) +{ +#define CELER_US_VISIT_CASE(TYPE) \ + case UnitSystem::TYPE: \ + return std::forward(func)(UnitSystemTraits{}); + + 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 diff --git a/src/celeritas/field/RZMapFieldInput.hh b/src/celeritas/field/RZMapFieldInput.hh index e54d3467ba..d31f1a4ff4 100644 --- a/src/celeritas/field/RZMapFieldInput.hh +++ b/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 field_z; //!< Flattened Z field component [bfield] std::vector field_r; //!< Flattened R field component [bfield] diff --git a/src/celeritas/field/RZMapFieldInputIO.json.cc b/src/celeritas/field/RZMapFieldInputIO.json.cc index 61faa2358f..4796715ab1 100644 --- a/src/celeritas/field/RZMapFieldInputIO.json.cc +++ b/src/celeritas/field/RZMapFieldInputIO.json.cc @@ -13,6 +13,7 @@ #include #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(); - if (units == "tesla" || units == "T") + auto const& ustr = iter->get(); + 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; + 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{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{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), diff --git a/src/celeritas/field/RZMapFieldParams.hh b/src/celeritas/field/RZMapFieldParams.hh index 312dfc505c..90d5a3491f 100644 --- a/src/celeritas/field/RZMapFieldParams.hh +++ b/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 { diff --git a/src/corecel/cont/Array.hh b/src/corecel/cont/Array.hh index e1e94f98d7..299e7cb8da 100644 --- a/src/corecel/cont/Array.hh +++ b/src/corecel/cont/Array.hh @@ -12,8 +12,6 @@ namespace celeritas { -#define CFIF_ CELER_FORCEINLINE_FUNCTION - //---------------------------------------------------------------------------// /*! * Fixed-size simple array for storage. @@ -52,24 +50,33 @@ struct Array //!@{ //! \name Element access - CFIF_ const_reference operator[](size_type i) const { return data_[i]; } - CFIF_ reference operator[](size_type i) { return data_[i]; } - CFIF_ const_reference front() const { return data_[0]; } - CFIF_ reference front() { return data_[0]; } - CFIF_ const_reference back() const { return data_[N - 1]; } - CFIF_ reference back() { return data_[N - 1]; } - CFIF_ const_pointer data() const { return data_; } - CFIF_ pointer data() { return data_; } + CELER_CONSTEXPR_FUNCTION const_reference operator[](size_type i) const + { + return data_[i]; + } + CELER_CONSTEXPR_FUNCTION reference operator[](size_type i) + { + return data_[i]; + } + CELER_CONSTEXPR_FUNCTION const_reference front() const { return data_[0]; } + CELER_CONSTEXPR_FUNCTION reference front() { return data_[0]; } + CELER_CONSTEXPR_FUNCTION const_reference back() const + { + return data_[N - 1]; + } + CELER_CONSTEXPR_FUNCTION reference back() { return data_[N - 1]; } + CELER_CONSTEXPR_FUNCTION const_pointer data() const { return data_; } + CELER_CONSTEXPR_FUNCTION pointer data() { return data_; } //!@} //!@{ //! \name Iterators - CFIF_ iterator begin() { return data_; } - CFIF_ iterator end() { return data_ + N; } - CFIF_ const_iterator begin() const { return data_; } - CFIF_ const_iterator end() const { return data_ + N; } - CFIF_ const_iterator cbegin() const { return data_; } - CFIF_ const_iterator cend() const { return data_ + N; } + CELER_CONSTEXPR_FUNCTION iterator begin() { return data_; } + CELER_CONSTEXPR_FUNCTION iterator end() { return data_ + N; } + CELER_CONSTEXPR_FUNCTION const_iterator begin() const { return data_; } + CELER_CONSTEXPR_FUNCTION const_iterator end() const { return data_ + N; } + CELER_CONSTEXPR_FUNCTION const_iterator cbegin() const { return data_; } + CELER_CONSTEXPR_FUNCTION const_iterator cend() const { return data_ + N; } //!@} //!@{ @@ -81,7 +88,7 @@ struct Array //!@{ //! \name Operations //! Fill the array with a constant value - CFIF_ void fill(const_reference value) + CELER_CONSTEXPR_FUNCTION void fill(const_reference value) { for (size_type i = 0; i != N; ++i) data_[i] = value; @@ -118,7 +125,5 @@ operator!=(Array const& lhs, Array const& rhs) return !(lhs == rhs); } -#undef CFIF_ - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/test/celeritas/UnitUtils.hh b/test/celeritas/UnitUtils.hh new file mode 100644 index 0000000000..bab280db3d --- /dev/null +++ b/test/celeritas/UnitUtils.hh @@ -0,0 +1,48 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/UnitUtils.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Types.hh" +#include "corecel/cont/Array.hh" +#include "celeritas/Quantities.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +//! Convert a value to centimeters from the native system +constexpr inline real_type to_cm(real_type v) +{ + return native_value_to(v).value(); +} + +//---------------------------------------------------------------------------// +//! Convert a value *from* centimeters to the native system +constexpr inline real_type from_cm(real_type v) +{ + return native_value_from(units::CmLength{v}); +} + +//---------------------------------------------------------------------------// +//! Convert an array to centimeters from the native system +constexpr inline Array to_cm(Array const& v) +{ + return {to_cm(v[0]), to_cm(v[1]), to_cm(v[2])}; +} + +//---------------------------------------------------------------------------// +//! Convert an array *from* centimeters to the native system +constexpr inline Array from_cm(Array const& v) +{ + return {from_cm(v[0]), from_cm(v[1]), from_cm(v[2])}; +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/Units.test.cc b/test/celeritas/Units.test.cc index a1726e610e..189e88ff9b 100644 --- a/test/celeritas/Units.test.cc +++ b/test/celeritas/Units.test.cc @@ -68,6 +68,19 @@ TEST(UnitsTest, traits) EXPECT_REAL_EQ(ClhepTraits::BField::value(), ClhepUnitBField::value()); } +//---------------------------------------------------------------------------// +TEST(UnitsTest, trait_visitor) +{ + auto get_length_str = [](auto utraits) { + using Length = typename decltype(utraits)::Length; + return Length::label(); + }; + + EXPECT_STREQ("cm", visit_unit_system(get_length_str, UnitSystem::cgs)); + EXPECT_STREQ("m", visit_unit_system(get_length_str, UnitSystem::si)); + EXPECT_STREQ("mm", visit_unit_system(get_length_str, UnitSystem::clhep)); +} + //---------------------------------------------------------------------------// } // namespace test } // namespace units diff --git a/test/celeritas/field/Fields.test.cc b/test/celeritas/field/Fields.test.cc index 03350700bc..54881e74ee 100644 --- a/test/celeritas/field/Fields.test.cc +++ b/test/celeritas/field/Fields.test.cc @@ -8,6 +8,7 @@ #include #include "corecel/cont/Range.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/field/RZMapField.hh" #include "celeritas/field/RZMapFieldInput.hh" #include "celeritas/field/RZMapFieldParams.hh" @@ -30,7 +31,8 @@ TEST(UniformZFieldTest, all) { UniformZField calc_field(123); - EXPECT_VEC_SOFT_EQ((Real3{0, 0, 123}), calc_field({100, -1, 0.5})); + EXPECT_VEC_SOFT_EQ((Real3{0, 0, 123}), + calc_field(from_cm(Real3{100, -1, 0.5}))); } TEST(UniformFieldTest, all) @@ -38,7 +40,7 @@ TEST(UniformFieldTest, all) Real3 field_vec{1, 2, 3}; UniformField calc_field(field_vec); - EXPECT_VEC_SOFT_EQ(field_vec, calc_field({100, -1, 0.5})); + EXPECT_VEC_SOFT_EQ(field_vec, calc_field(from_cm(Real3{100, -1, 0.5}))); } TEST(CMSParameterizedFieldTest, all) @@ -47,8 +49,8 @@ TEST(CMSParameterizedFieldTest, all) CMSParameterizedField calc_field; int const nsamples = 8; - real_type delta_z = 25.0; - real_type delta_r = 12.0; + real_type const delta_z = from_cm(25.0); + real_type const delta_r = from_cm(12.0); std::vector actual; @@ -57,7 +59,7 @@ TEST(CMSParameterizedFieldTest, all) Real3 field = calc_field(Real3{i * delta_r, i * delta_r, i * delta_z}); for (real_type f : field) { - actual.push_back(f / units::tesla); + actual.push_back(native_value_to(f).value()); } } @@ -107,8 +109,8 @@ TEST_F(RZMapFieldTest, all) RZMapField calc_field(field_map.host_ref()); int const nsamples = 8; - real_type delta_z = 25.0; - real_type delta_r = 12.0; + real_type delta_z = from_cm(25.0); + real_type delta_r = from_cm(12.0); std::vector actual; @@ -118,7 +120,7 @@ TEST_F(RZMapFieldTest, all) for (real_type f : field) { // Reference result is in [T]: convert from native units - actual.push_back(f / units::tesla); + actual.push_back(native_value_to(f).value()); } }