diff --git a/src/accel/CMakeLists.txt b/src/accel/CMakeLists.txt index 3c65522c40..3facbb754b 100644 --- a/src/accel/CMakeLists.txt +++ b/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 diff --git a/src/accel/CartMapMagneticField.cc b/src/accel/CartMapMagneticField.cc new file mode 100644 index 0000000000..08669711cc --- /dev/null +++ b/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 +#include +#include +#include + +#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()}) + { + 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(Axis::size_) + * total_points); + + Array const dims{params.x.num, + params.y.num, + params.z.num, + static_cast(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{x, y, z, 0}; + }; + + // Field converter for Cartesian coordinates (no transformation needed) + auto field_converter = [](Array 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(result[i]); + field[i] = convert_to_geant(ft.value(), CLHEP::tesla); + } +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/accel/CartMapMagneticField.hh b/src/accel/CartMapMagneticField.hh new file mode 100644 index 0000000000..3eaaa04516 --- /dev/null +++ b/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 +#include + +#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 x{}; //!< X-axis grid specification + AxisGrid y{}; //!< Y-axis grid specification + AxisGrid 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; + //!@} + + // 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 pimpl_; +}; + +//---------------------------------------------------------------------------// + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/accel/CylMapMagneticField.cc b/src/accel/CylMapMagneticField.cc index 770433ec2e..1b03c58046 100644 --- a/src/accel/CylMapMagneticField.cc +++ b/src/accel/CylMapMagneticField.cc @@ -18,7 +18,6 @@ #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" @@ -26,9 +25,10 @@ #include "celeritas/field/CylMapFieldInput.hh" #include "celeritas/field/CylMapFieldParams.hh" +#include "detail/MagneticFieldUtils.hh" + namespace celeritas { -//---------------------------------------------------------------------------// namespace { //---------------------------------------------------------------------------// @@ -63,7 +63,7 @@ MakeCylMapFieldInput(std::vector 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 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(CylAxis::size_) * total_points); Array const dims{ nr, nphi, nz, static_cast(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 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 pos - = {r * std::cos(phi), r * std::sin(phi), z_grid[iz], 0}; - g4field->GetFieldValue(pos.data(), bfield.data()); - EnumArray 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{r * std::cos(phi), r * std::sin(phi), z, 0}; + }; + + // Field converter for cylindrical coordinates (requires coordinate + // transformation) + auto field_converter = [](Array const& bfield, + real_type* cur_bfield) { + EnumArray 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; } diff --git a/src/accel/detail/MagneticFieldUtils.hh b/src/accel/detail/MagneticFieldUtils.hh new file mode 100644 index 0000000000..cab749b537 --- /dev/null +++ b/src/accel/detail/MagneticFieldUtils.hh @@ -0,0 +1,80 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file accel/detail/MagneticFieldUtils.hh +//! \brief Utilities for magnetic field map creation +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Types.hh" +#include "corecel/cont/Array.hh" +#include "corecel/data/HyperslabIndexer.hh" +#include "geocel/GeantGeoUtils.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Common field sampling setup and execution. + * + * This function encapsulates the shared pattern of: + * 1. Getting the G4 field + * 2. Setting up the HyperslabIndexer + * 3. Sampling the field on a grid + * + * \param field_data Output parameter array to store field values (must be + * pre-allocated with size equal to the product of all dims) + * \param dims Grid dimensions + * \param calc_position Callable that computes position given + * grid indices. Must have signature: + * Array(size_type, size_type, + * size_type) returning [x, y, z, 0] coordinates + * \param convert_field Callable that converts field from G4 to + * native units in the correct coordinate space. Must + * have signature: void(Array const&, + * real_type*) taking G4 field [Bx, By, Bz] and writing + * converted values to output pointer + */ +template +inline void setup_and_sample_field(real_type* field_data, + Array const& dims, + PositionCalc const& calc_position, + FieldConverter const& convert_field) +{ + HyperslabIndexer const flat_index{dims}; + G4Field const* g4field = celeritas::geant_field(); + CELER_VALIDATE(g4field, + << "no Geant4 global field has been set: cannot build " + "magnetic field map"); + + Array bfield; + + for (size_type i = 0; i < dims[0]; ++i) + { + for (size_type j = 0; j < dims[1]; ++j) + { + for (size_type k = 0; k < dims[2]; ++k) + { + // Calculate position for this grid point + Array pos = calc_position(i, j, k); + + // Sample field at this position + g4field->GetFieldValue(pos.data(), bfield.data()); + + // Convert and store field values + auto* cur_bfield = field_data + flat_index(i, j, k, 0); + convert_field(bfield, cur_bfield); + } + } + } +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/field/CartMapFieldInput.hh b/src/celeritas/field/CartMapFieldInput.hh index dedcd3598f..2044c5e2ca 100644 --- a/src/celeritas/field/CartMapFieldInput.hh +++ b/src/celeritas/field/CartMapFieldInput.hh @@ -16,6 +16,19 @@ namespace celeritas { +//---------------------------------------------------------------------------// +//! Grid specification for a single axis +template +struct AxisGrid +{ + T min{}; //!< Minimum coordinate value + T max{}; //!< Maximum coordinate value + size_type num{}; //!< Number of grid points + + //! Check if parameters are valid + explicit operator bool() const { return max > min && num > 1; } +}; + //---------------------------------------------------------------------------// /*! * Input data for a magnetic X-Y-Z vector field stored on an X-Y-Z grid. @@ -30,17 +43,9 @@ namespace celeritas */ struct CartMapFieldInput { - real_type min_x; //!< Minimum X grid point [len] - real_type max_x; //!< Maximum X grid point [len] - size_type num_x; //!< Number of X grid points - - real_type min_y; //!< Minimum Y grid point [len] - real_type max_y; //!< Maximum Y grid point [len] - size_type num_y; //!< Number of Y grid points - - real_type min_z; //!< Minimum Z grid point [len] - real_type max_z; //!< Maximum Z grid point [len] - size_type num_z; //!< Number of Z grid points + AxisGrid x; //!< X-axis grid specification [len] + AxisGrid y; //!< Y-axis grid specification [len] + AxisGrid z; //!< Z-axis grid specification [len] std::vector field; //!< Flattened X-Y-Z field component //!< [bfield] @@ -51,16 +56,10 @@ struct CartMapFieldInput //! Whether all data are assigned and valid explicit operator bool() const { - // clang-format off - return - (max_x >= min_x) - && (num_x >= 2) - && (max_y >= min_y) - && (num_y >= 2) - && (max_z >= min_z) - && (num_z >= 2) - && (field.size() == static_cast(Axis::size_) * num_x * num_y * num_z); - // clang-format on + return x && y && z + && (field.size() + == static_cast(Axis::size_) * x.num * y.num + * z.num); } }; diff --git a/src/celeritas/field/CartMapFieldParams.covfie.cc b/src/celeritas/field/CartMapFieldParams.covfie.cc index 9420d3695f..78ad1b2bf9 100644 --- a/src/celeritas/field/CartMapFieldParams.covfie.cc +++ b/src/celeritas/field/CartMapFieldParams.covfie.cc @@ -34,9 +34,9 @@ struct CartMapFieldParams::Impl : host_{[&inp] { HostVal host; - Array const dims{inp.num_x, - inp.num_y, - inp.num_z, + Array const dims{inp.x.num, + inp.y.num, + inp.z.num, static_cast(Axis::size_)}; HyperslabIndexer const flat_index{dims}; @@ -45,14 +45,14 @@ struct CartMapFieldParams::Impl builder_t builder{covfie::make_parameter_pack( builder_t::backend_t::configuration_t{ - inp.num_x, inp.num_y, inp.num_z})}; + inp.x.num, inp.y.num, inp.z.num})}; builder_t::view_t builder_view{builder}; // fill the covfie field data - for (auto ix : range(inp.num_x)) + for (auto ix : range(inp.x.num)) { - for (auto iy : range(inp.num_y)) + for (auto iy : range(inp.y.num)) { - for (auto iz : range(inp.num_z)) + for (auto iz : range(inp.z.num)) { auto* fv = builder_view.at(ix, iy, iz).begin(); auto* finp = inp.field.data() @@ -66,17 +66,17 @@ struct CartMapFieldParams::Impl using field_real_type = CartMapField::real_type; auto affine_translate = covfie::algebra::affine<3>::translation( - static_cast(-inp.min_x), - static_cast(-inp.min_y), - static_cast(-inp.min_z)); + static_cast(-inp.x.min), + static_cast(-inp.y.min), + static_cast(-inp.z.min)); auto affine_scale = covfie::algebra::affine<3>::scaling( - static_cast((inp.num_x - 1) - / (inp.max_x - inp.min_x)), - static_cast((inp.num_y - 1) - / (inp.max_y - inp.min_y)), - static_cast((inp.num_z - 1) - / (inp.max_z - inp.min_z))); + static_cast((inp.x.num - 1) + / (inp.x.max - inp.x.min)), + static_cast((inp.y.num - 1) + / (inp.y.max - inp.y.min)), + static_cast((inp.z.num - 1) + / (inp.z.max - inp.z.min))); using field_t = detail::CovfieFieldTraits::field_t; host.field = std::make_unique(covfie::make_parameter_pack( diff --git a/test/celeritas/field/Fields.test.cc b/test/celeritas/field/Fields.test.cc index 9ddbe425fa..0cad868c11 100644 --- a/test/celeritas/field/Fields.test.cc +++ b/test/celeritas/field/Fields.test.cc @@ -289,29 +289,29 @@ using CartMapFieldTest = ::celeritas::test::Test; TEST_F(CartMapFieldTest, all) { CartMapFieldInput inp; - inp.min_x = -2750; - inp.max_x = 2750; - inp.num_x = static_cast(inp.max_x * 2 / 100); - inp.min_y = -2750; - inp.max_y = 2750; - inp.num_y = static_cast(inp.max_y * 2 / 100); - inp.min_z = -6350; - inp.max_z = 6350; - inp.num_z = static_cast(inp.max_z * 2 / 100); + inp.x.min = -2750; + inp.x.max = 2750; + inp.x.num = static_cast(inp.x.max * 2 / 100); + inp.y.min = -2750; + inp.y.max = 2750; + inp.y.num = static_cast(inp.y.max * 2 / 100); + inp.z.min = -6350; + inp.z.max = 6350; + inp.z.num = static_cast(inp.z.max * 2 / 100); Array const dims{ - inp.num_x, inp.num_y, inp.num_z, static_cast(Axis::size_)}; - size_type const total_points = inp.num_x * inp.num_y * inp.num_z; + inp.x.num, inp.y.num, inp.z.num, static_cast(Axis::size_)}; + size_type const total_points = inp.x.num * inp.y.num * inp.z.num; // Resize each component of the field inp.field.resize(static_cast(Axis::size_) * total_points); // Fill with a simple field pattern HyperslabIndexer const flat_index{dims}; - for (size_type x = 0; x < inp.num_x; ++x) + for (size_type x = 0; x < inp.x.num; ++x) { - for (size_type y = 0; y < inp.num_y; ++y) + for (size_type y = 0; y < inp.y.num; ++y) { - for (size_type z = 0; z < inp.num_z; ++z) + for (size_type z = 0; z < inp.z.num; ++z) { auto arr = inp.field.begin() + flat_index(x, y, z, 0); arr[static_cast(Axis::x)] = std::cos(x); @@ -335,19 +335,19 @@ TEST_F(CartMapFieldTest, all) for (size_type ix = 0; ix < nx_samples; ++ix) { - real_type x = inp.min_x - + ix * (inp.max_x - inp.min_x) / (nx_samples - 1); - x = std::min(x, inp.max_x - 1); + real_type x = inp.x.min + + ix * (inp.x.max - inp.x.min) / (nx_samples - 1); + x = std::min(x, inp.x.max - 1); for (size_type iy = 0; iy < ny_samples; ++iy) { - real_type y = inp.min_y - + iy * (inp.max_y - inp.min_y) / (ny_samples - 1); - y = std::min(y, inp.max_y - 1); + real_type y = inp.y.min + + iy * (inp.y.max - inp.y.min) / (ny_samples - 1); + y = std::min(y, inp.y.max - 1); for (size_type iz = 0; iz < nz_samples; ++iz) { - real_type z = inp.min_z - + iz * (inp.max_z - inp.min_z) / (nz_samples - 1); - z = std::min(z, inp.max_z - 1); + real_type z = inp.z.min + + iz * (inp.z.max - inp.z.min) / (nz_samples - 1); + z = std::min(z, inp.z.max - 1); Real3 field = calc_field({x, y, z}); for (real_type f : field)