Skip to content
2 changes: 2 additions & 0 deletions src/accel/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand All @@ -29,6 +30,7 @@ list(APPEND SOURCES
IntegrationBase.cc
Logger.cc
LocalTransporter.cc
CartMapMagneticField.cc
CylMapMagneticField.cc
SetupOptions.cc
SetupOptionsMessenger.cc
Expand Down
156 changes: 156 additions & 0 deletions src/accel/CartMapMagneticField.cc
Original file line number Diff line number Diff line change
@@ -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
68 changes: 68 additions & 0 deletions src/accel/CartMapMagneticField.hh
Original file line number Diff line number Diff line change
@@ -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
60 changes: 28 additions & 32 deletions src/accel/CylMapMagneticField.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -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),
Expand All @@ -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;
}
Expand Down
Loading
Loading