celeritas/field/MagFieldEquation.hh

0001 //---------------------------------*-CUDA-*----------------------------------//
0002 // Copyright 2020-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/field/MagFieldEquation.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/Types.hh"
0013 #include "corecel/math/Algorithms.hh"
0014 #include "corecel/math/ArrayOperators.hh"
0015 #include "celeritas/Constants.hh"
0016 #include "celeritas/Quantities.hh"
0017
0018 #include "Types.hh"
0019
0020 #include "detail/FieldUtils.hh"
0021
0022 namespace celeritas
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Evaluate the right hand side of the Lorentz equation.
0027  *
0028  * The templated \c FieldT must be a function-like object with the signature
0029  * \code
0030  * Real3 (*)(const Real3&)
0031  * \endcode
0032  * which returns a magnetic field vector at a given position. The field
0033  * strength is in Celeritas native units, not Tesla.
0034  *
0035  * Calling an instance of this class calculates the local derivatives of
0036  * position and momentum (i.e.  direction and force) based on the given
0037  * magnetic field state.
0038  *
0039  * \f[
0040     m \frac{d^2 \vec{x}}{d t^2} = (q/c)(\vec{v} \times  \vec{B})
0041     s = |v|t
0042     \vec{y} = d\vec{x}/ds
0043     \frac{d\vec{x}}{ds} = \vec{v}/|v|
0044     \frac{d\vec{y}}{ds} = (q/pc)(\vec{y} \times \vec{B})
0045    \f]
0046  *
0047  */
0048 template<class FieldT>
0049 class MagFieldEquation
0050 {
0051   public:
0052     //!@{
0053     //! \name Type aliases
0054     using Field_t = FieldT;
0055     //!@}
0056
0057   public:
0058     // Construct with a magnetic field
0059     inline CELER_FUNCTION
0060     MagFieldEquation(FieldT&& field, units::ElementaryCharge q);
0061
0062     // Evaluate the right hand side of the field equation
0063     inline CELER_FUNCTION OdeState operator()(OdeState const& y) const;
0064
0065   private:
0066     // Field evaluator
0067     Field_t calc_field_;
0068
0069     // The (Lorentz) coefficent in 1/OdeState::MomentumUnits
0070     real_type coeffi_;
0071 };
0072
0073 //---------------------------------------------------------------------------//
0074 // DEDUCTION GUIDES
0075 //---------------------------------------------------------------------------//
0076 template<class FieldT>
0077 CELER_FUNCTION
0078 MagFieldEquation(FieldT&&, units::ElementaryCharge) -> MagFieldEquation<FieldT>;
0079
0080 //---------------------------------------------------------------------------//
0081 // INLINE DEFINITIONS
0082 //---------------------------------------------------------------------------//
0083 /*!
0084  * Construct with a magnetic field equation and particle charge.
0085  *
0086  * The internal coefficient is based on Celeritas native units and the
0087  * "natural" unit system used by the \c ParticleTrackView.
0088  */
0089 template<class FieldT>
0090 CELER_FUNCTION
0091 MagFieldEquation<FieldT>::MagFieldEquation(FieldT&& field,
0092                                            units::ElementaryCharge charge)
0093     : calc_field_(::celeritas::forward<FieldT>(field))
0094     , coeffi_{native_value_from(charge)
0095               / native_value_from(OdeState::MomentumUnits{1})}
0096 {
0097 }
0098
0099 //---------------------------------------------------------------------------//
0100 /*!
0101  * Evaluate the right hand side of the Lorentz equation.
0102  */
0103 template<class FieldT>
0104 CELER_FUNCTION auto
0105 MagFieldEquation<FieldT>::operator()(OdeState const& y) const -> OdeState
0106 {
0107     CELER_EXPECT(y.mom[0] != 0 || y.mom[1] != 0 || y.mom[2] != 0);
0108     real_type momentum_inv = celeritas::rsqrt(dot_product(y.mom, y.mom));
0109
0110     // Evaluate the rate of change in particle's position per unit length: this
0111     // is just the direction
0112     OdeState result;
0113     result.pos = momentum_inv * y.mom;
0114
0115     // Calculate the magnetic field value at the current position
0116     // to calculate the force on the particle
0117     result.mom = (coeffi_ * momentum_inv)
0118                  * cross_product(y.mom, calc_field_(y.pos));
0119
0120     return result;
0121 }
0122
0123 //---------------------------------------------------------------------------//
0124 }  // namespace celeritas