Add in-situ cartesian field map construction from Geant4

src/accel/CMakeLists.txt

@@ -6,6 +6,7 @@
 set(SOURCES)
 set(PRIVATE_DEPS
   nlohmann_json::nlohmann_json
+  Celeritas::Extcovfie
   Celeritas::ExtThrust # TODO: add thrust exception conversion to corecel?
   Celeritas::ExtOpenMP # TODO: provide thread utils in corecel
 )
@@ -29,6 +30,7 @@ list(APPEND SOURCES
   IntegrationBase.cc
   Logger.cc
   LocalTransporter.cc
+  CartMapMagneticField.cc
   CylMapMagneticField.cc
   SetupOptions.cc
   SetupOptionsMessenger.cc

src/accel/CartMapMagneticField.cc

@@ -0,0 +1,156 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file accel/CartMapMagneticField.cc
+//---------------------------------------------------------------------------//
+
+#include "CartMapMagneticField.hh"
+
+#include <algorithm>
+#include <CLHEP/Units/SystemOfUnits.h>
+#include <G4MagneticField.hh>
+#include <corecel/Assert.hh>
+
+#include "corecel/Types.hh"
+#include "corecel/cont/Array.hh"
+#include "corecel/math/ArrayUtils.hh"
+#include "geocel/g4/Convert.hh"
+#include "celeritas/Types.hh"
+#include "celeritas/ext/GeantUnits.hh"
+#include "celeritas/field/CartMapField.hh"
+#include "celeritas/field/CartMapFieldInput.hh"
+#include "celeritas/field/CartMapFieldParams.hh"
+
+#include "detail/MagneticFieldUtils.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+// PIMPL IMPLEMENTATION
+//---------------------------------------------------------------------------//
+
+/*!
+ * Implementation struct for CartMapMagneticField.
+ *
+ * This hides the C++20 dependency (CartMapField) from the header file.
+ */
+struct CartMapMagneticField::Impl
+{
+    CartMapMagneticField::SPConstFieldParams params;
+    CartMapField calc_field;
+
+    explicit Impl(CartMapMagneticField::SPConstFieldParams field_params)
+        : params(std::move(field_params))
+        , calc_field(CartMapField{params->ref<MemSpace::native>()})
+    {
+        CELER_EXPECT(params);
+    }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Generates input for CartMapField params with configurable uniform grid
+ * dimensions in native Geant4 units. This must be called after
+ * G4RunManager::Initialize as it will retrieve the G4FieldManager's field
+ * to sample it.
+ */
+CartMapFieldParams::Input
+MakeCartMapFieldInput(CartMapFieldGridParams const& params)
+{
+    // Validate input parameters
+    CELER_VALIDATE(params, << "invalid CartMapFieldGridParams provided");
+
+    CartMapFieldParams::Input field_input;
+
+    // Convert from Geant4 units to native units
+    field_input.x.min = convert_from_geant(params.x.min, clhep_length);
+    field_input.x.max = convert_from_geant(params.x.max, clhep_length);
+    field_input.x.num = params.x.num;
+
+    field_input.y.min = convert_from_geant(params.y.min, clhep_length);
+    field_input.y.max = convert_from_geant(params.y.max, clhep_length);
+    field_input.y.num = params.y.num;
+
+    field_input.z.min = convert_from_geant(params.z.min, clhep_length);
+    field_input.z.max = convert_from_geant(params.z.max, clhep_length);
+    field_input.z.num = params.z.num;
+
+    // Prepare field data storage
+    size_type const total_points = params.x.num * params.y.num * params.z.num;
+    field_input.field.resize(static_cast<size_type>(Axis::size_)
+                             * total_points);
+
+    Array<size_type, 4> const dims{params.x.num,
+                                   params.y.num,
+                                   params.z.num,
+                                   static_cast<size_type>(Axis::size_)};
+
+    // Calculate grid spacing
+    G4double const dx = (params.x.max - params.x.min) / (params.x.num - 1);
+    G4double const dy = (params.y.max - params.y.min) / (params.y.num - 1);
+    G4double const dz = (params.z.max - params.z.min) / (params.z.num - 1);
+
+    // Position calculator for Cartesian grid
+    auto position_calculator = [&](size_type ix, size_type iy, size_type iz) {
+        G4double x = params.x.min + ix * dx;
+        G4double y = params.y.min + iy * dy;
+        G4double z = params.z.min + iz * dz;
+        return Array<G4double, 4>{x, y, z, 0};
+    };
+
+    // Field converter for Cartesian coordinates (no transformation needed)
+    auto field_converter = [](Array<G4double, 3> const& bfield,
+                              real_type* cur_bfield) {
+        auto bfield_native = convert_from_geant(bfield.data(), clhep_field);
+        std::copy(bfield_native.cbegin(), bfield_native.cend(), cur_bfield);
+    };
+
+    // Sample field using common utility
+    detail::setup_and_sample_field(
+        field_input.field.data(), dims, position_calculator, field_converter);
+
+    CELER_ENSURE(field_input);
+    return field_input;
+}
+
+//---------------------------------------------------------------------------//
+// CARTMAPMAGNETIC FIELD IMPLEMENTATION
+//---------------------------------------------------------------------------//
+
+/*!
+ * Custom deleter implementation for PIMPL idiom.
+ */
+void CartMapMagneticField::ImplDeleter::operator()(Impl* ptr) const
+{
+    delete ptr;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with the Celeritas shared CartMapFieldParams.
+ */
+CartMapMagneticField::CartMapMagneticField(SPConstFieldParams field_params)
+    : pimpl_(new Impl(std::move(field_params)))
+{
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the magnetic field vector at the given position.
+ */
+void CartMapMagneticField::GetFieldValue(G4double const pos[3],
+                                         G4double* field) const
+{
+    // Calculate the magnetic field value in the native Celeritas unit system
+    Real3 result = pimpl_->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
+        auto ft = native_value_to<units::FieldTesla>(result[i]);
+        field[i] = convert_to_geant(ft.value(), CLHEP::tesla);
+    }
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/accel/CartMapMagneticField.hh

@@ -0,0 +1,68 @@
+//------------------------------- -*- C++ -*- -------------------------------//
+// Copyright Celeritas contributors: see top-level COPYRIGHT file for details
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file accel/CartMapMagneticField.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include <memory>
+#include <G4MagneticField.hh>
+
+#include "corecel/Types.hh"
+#include "celeritas/field/CartMapFieldInput.hh"
+#include "celeritas/field/CartMapFieldParams.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+//! POD struct for CartMap field grid parameters
+struct CartMapFieldGridParams
+{
+    AxisGrid<real_type> x{};  //!< X-axis grid specification
+    AxisGrid<real_type> y{};  //!< Y-axis grid specification
+    AxisGrid<real_type> z{};  //!< Z-axis grid specification
+
+    //! Check if parameters are valid for field generation
+    explicit operator bool() const { return x && y && z; }
+};
+
+//---------------------------------------------------------------------------//
+// Generate field input with user-defined uniform grid
+CartMapFieldParams::Input
+MakeCartMapFieldInput(CartMapFieldGridParams const& params);
+
+//---------------------------------------------------------------------------//
+/*!
+ * A user magnetic field equivalent to celeritas::CartMapField.
+ */
+class CartMapMagneticField : public G4MagneticField
+{
+  public:
+    //!@{
+    //! \name Type aliases
+    using SPConstFieldParams = std::shared_ptr<CartMapFieldParams const>;
+    //!@}
+
+    // Construct with CartMapFieldParams
+    explicit CartMapMagneticField(SPConstFieldParams field_params);
+
+    // Calculate values of the magnetic field vector
+    void GetFieldValue(G4double const point[3], G4double* field) const override;
+
+  private:
+    // Forward declaration for pImpl
+    struct Impl;
+
+    // Custom deleter for pImpl
+    struct ImplDeleter
+    {
+        void operator()(Impl* ptr) const;
+    };
+    std::unique_ptr<Impl, ImplDeleter> pimpl_;
+};
+
+//---------------------------------------------------------------------------//
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas

src/accel/CylMapMagneticField.cc

@@ -18,17 +18,17 @@
 #include "corecel/cont/EnumArray.hh"
 #include "corecel/math/Quantity.hh"
 #include "corecel/math/Turn.hh"
-#include "geocel/GeantGeoUtils.hh"
 #include "geocel/g4/Convert.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/ext/GeantUnits.hh"
 #include "celeritas/field/CylMapFieldInput.hh"
 #include "celeritas/field/CylMapFieldParams.hh"
 
+#include "detail/MagneticFieldUtils.hh"
+
 namespace celeritas
 {
-//---------------------------------------------------------------------------//
 namespace
 {
 //---------------------------------------------------------------------------//
@@ -63,7 +63,7 @@ MakeCylMapFieldInput(std::vector<G4double> const& r_grid,
     field_input.grid_phi.reserve(phi_values.size());
     field_input.grid_z.reserve(z_grid.size());
 
-    // Convert from geant
+    // Convert from geant units to native units
     std::transform(r_grid.cbegin(),
                    r_grid.cend(),
                    std::back_inserter(field_input.grid_r),
@@ -83,41 +83,37 @@ MakeCylMapFieldInput(std::vector<G4double> const& r_grid,
     size_type const nz = field_input.grid_z.size();
     size_type const total_points = nr * nphi * nz;
 
+    // Prepare field data storage
     field_input.field.resize(static_cast<size_type>(CylAxis::size_)
                              * total_points);
 
     Array<size_type, 4> const dims{
         nr, nphi, nz, static_cast<size_type>(CylAxis::size_)};
-    HyperslabIndexer const flat_index{dims};
-
-    G4Field const* g4field = celeritas::geant_field();
-    CELER_VALIDATE(g4field,
-                   << "no Geant4 global field has been set: cannot build "
-                      "CylMapMagneticField");
-    Array<G4double, 3> bfield;
-    for (size_type ir = 0; ir < nr; ++ir)
-    {
+
+    // Position calculator for cylindrical grid
+    auto position_calculator = [&](size_type ir, size_type iphi, size_type iz) {
         auto r = r_grid[ir];
-        for (size_type iphi = 0; iphi < nphi; ++iphi)
-        {
-            auto phi = phi_values[iphi];
-            for (size_type iz = 0; iz < nz; ++iz)
-            {
-                auto* cur_bfield = field_input.field.data()
-                                   + flat_index(ir, iphi, iz, 0);
-                Array<G4double, 4> pos
-                    = {r * std::cos(phi), r * std::sin(phi), z_grid[iz], 0};
-                g4field->GetFieldValue(pos.data(), bfield.data());
-                EnumArray<CylAxis, G4double> bfield_cyl;
-                cartesian_to_cylindrical(bfield, bfield_cyl);
-                auto bfield_cyl_native
-                    = convert_from_geant(bfield_cyl.data(), clhep_field);
-                std::copy(bfield_cyl_native.cbegin(),
-                          bfield_cyl_native.cend(),
-                          cur_bfield);
-            }
-        }
-    }
+        auto phi = phi_values[iphi];
+        auto z = z_grid[iz];
+        return Array<G4double, 4>{r * std::cos(phi), r * std::sin(phi), z, 0};
+    };
+
+    // Field converter for cylindrical coordinates (requires coordinate
+    // transformation)
+    auto field_converter = [](Array<G4double, 3> const& bfield,
+                              real_type* cur_bfield) {
+        EnumArray<CylAxis, G4double> bfield_cyl;
+        cartesian_to_cylindrical(bfield, bfield_cyl);
+        auto bfield_cyl_native
+            = convert_from_geant(bfield_cyl.data(), clhep_field);
+        std::copy(
+            bfield_cyl_native.cbegin(), bfield_cyl_native.cend(), cur_bfield);
+    };
+
+    // Sample field using common utility
+    detail::setup_and_sample_field(
+        field_input.field.data(), dims, position_calculator, field_converter);
+
     CELER_ENSURE(field_input);
     return field_input;
 }