celeritas/field/CylMapField.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/field/CylMapField.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Macros.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/cont/EnumArray.hh"
0013 #include "corecel/cont/Range.hh"
0014 #include "corecel/grid/FindInterp.hh"
0015 #include "corecel/grid/NonuniformGrid.hh"
0016 #include "corecel/math/Algorithms.hh"
0017 #include "corecel/math/Quantity.hh"
0018 #include "corecel/math/Turn.hh"
0019 #include "celeritas/Types.hh"
0020
0021 #include "CylMapFieldData.hh"
0022
0023 namespace celeritas
0024 {
0025 //---------------------------------------------------------------------------//
0026 /*!
0027  * Interpolate a magnetic field vector on an r/phi/z grid.
0028  *
0029  * The field vector is stored as a cartesian \f$(x,y,z)\f$ value on the
0030  * cylindrical mesh grid points, and trilinear interpolation is performed
0031  * within each grid cell. The value outside the grid is zero.
0032  *
0033  * Currently the grid requires a full \f$2\pi\f$ azimuthal grid.
0034  */
0035 class CylMapField
0036 {
0037   public:
0038     //!@{
0039     //! \name Type aliases
0040     using real_type = cylmap_real_type;
0041     using Real3 = Array<celeritas::real_type, 3>;
0042     using FieldParamsRef = NativeCRef<CylMapFieldParamsData>;
0043     //!@}
0044
0045   public:
0046     // Construct with the shared map data
0047     inline CELER_FUNCTION explicit CylMapField(FieldParamsRef const& shared);
0048
0049     // Evaluate the magnetic field value for the given position
0050     CELER_FUNCTION
0051     inline Real3 operator()(Real3 const& pos) const;
0052
0053   private:
0054     // Shared constant field map
0055     FieldParamsRef const& params_;
0056
0057     NonuniformGrid<real_type> const grid_r_;
0058     NonuniformGrid<real_type> const grid_phi_;
0059     NonuniformGrid<real_type> const grid_z_;
0060 };
0061
0062 //---------------------------------------------------------------------------//
0063 // INLINE DEFINITIONS
0064 //---------------------------------------------------------------------------//
0065 /*!
0066  * Construct with the shared magnetic field map data.
0067  */
0068 CELER_FUNCTION
0069 CylMapField::CylMapField(FieldParamsRef const& params)
0070     : params_{params}
0071     , grid_r_{params_.grids.axes[CylAxis::r], params_.grids.storage}
0072     , grid_phi_{params_.grids.axes[CylAxis::phi], params_.grids.storage}
0073     , grid_z_{params_.grids.axes[CylAxis::z], params_.grids.storage}
0074 {
0075 }
0076
0077 //---------------------------------------------------------------------------//
0078 /*!
0079  * Calculate the magnetic field vector for the given position.
0080  *
0081  * This does a 3-D interpolation on the input grid and reconstructs the
0082  * magnetic field vector from the stored R, Z, and Phi components of the field.
0083  * The result is in the native Celeritas unit system.
0084  */
0085 CELER_FUNCTION auto CylMapField::operator()(Real3 const& pos) const -> Real3
0086 {
0087     CELER_ENSURE(params_);
0088
0089     Array<real_type, 3> value{0, 0, 0};
0090
0091     // Convert Cartesian to cylindrical coordinates
0092     real_type r = hypot(pos[0], pos[1]);
0093     // Ensure phi is in [0, 2\f$\pi\f$)
0094     auto phi = atan2turn<real_type>(pos[1], pos[0]);
0095     if (phi < zero_quantity())
0096     {
0097         phi.value() += 1;
0098     }
0099     if (CELER_UNLIKELY(phi.value() == 1))
0100     {
0101         // Make sure phi is in [0, 1). If phi is a negative value smaller
0102         // than machine epsilon, adding 1 will result in phi equal to 1
0103         phi.value() = 0;
0104     }
0105     CELER_ASSERT(phi >= zero_quantity() && phi.value() < 1);
0106
0107     // Check if point is within field map bounds
0108     if (!params_.valid(r, phi, pos[2]))
0109         return {0, 0, 0};
0110
0111     // Find interpolation points for given r, phi, z
0112     auto [ir, wr1] = find_interp<NonuniformGrid<real_type>>(grid_r_, r);
0113     auto [iphi, wphi1]
0114         = find_interp<NonuniformGrid<real_type>>(grid_phi_, phi.value());
0115     auto [iz, wz1] = find_interp<NonuniformGrid<real_type>>(grid_z_, pos[2]);
0116
0117     auto get_field = [this](size_type ir, size_type iphi, size_type iz) {
0118         return params_.fieldmap[params_.id(ir, iphi, iz)];
0119     };
0120
0121     EnumArray<CylAxis, real_type> interp_field;
0122
0123     for (auto axis : range(CylAxis::size_))
0124     {
0125         // clang-format off
0126         real_type v000 = get_field(ir,     iphi    ,     iz)[axis];
0127         real_type v001 = get_field(ir,     iphi    , iz + 1)[axis];
0128         real_type v010 = get_field(ir,     iphi + 1,     iz)[axis];
0129         real_type v011 = get_field(ir,     iphi + 1, iz + 1)[axis];
0130         real_type v100 = get_field(ir + 1, iphi    ,     iz)[axis];
0131         real_type v101 = get_field(ir + 1, iphi    , iz + 1)[axis];
0132         real_type v110 = get_field(ir + 1, iphi + 1,     iz)[axis];
0133         real_type v111 = get_field(ir + 1, iphi + 1, iz + 1)[axis];
0134         // clang-format on
0135         // Trilinear interpolation formula for the current component
0136         interp_field[axis]
0137             = (1 - wr1)
0138                   * ((1 - wphi1) * ((1 - wz1) * v000 + wz1 * v001)
0139                      + wphi1 * ((1 - wz1) * v010 + wz1 * v011))
0140               + wr1
0141                     * ((1 - wphi1) * ((1 - wz1) * v100 + wz1 * v101)
0142                        + wphi1 * ((1 - wz1) * v110 + wz1 * v111));
0143     }
0144
0145     // Project cylindrical components to Cartesian coordinates
0146     // default for r == 0
0147     real_type cos_phi = 1;
0148     real_type sin_phi = 0;
0149
0150     if (r != 0)
0151     {
0152         cos_phi = pos[0] / r;
0153         sin_phi = pos[1] / r;
0154     }
0155     value[0] = interp_field[CylAxis::r] * cos_phi
0156                - interp_field[CylAxis::phi] * sin_phi;
0157     value[1] = interp_field[CylAxis::r] * sin_phi
0158                + interp_field[CylAxis::phi] * cos_phi;
0159     value[2] = interp_field[CylAxis::z];
0160
0161     return {value[0], value[1], value[2]};
0162 }
0163
0164 //---------------------------------------------------------------------------//
0165 }  // namespace celeritas