field/detail/FieldUtils.hh

0001 //----------------------------------*-C++-*----------------------------------//
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/detail/FieldUtils.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011 #include <iostream>
0012
0013 #include "corecel/Assert.hh"
0014 #include "corecel/cont/Array.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/math/ArrayOperators.hh"
0017 #include "corecel/math/ArrayUtils.hh"
0018 #include "celeritas/Types.hh"
0019
0020 #include "../Types.hh"
0021
0022 namespace celeritas
0023 {
0024 namespace detail
0025 {
0026 //---------------------------------------------------------------------------//
0027 // HELPER STRUCTS
0028 //---------------------------------------------------------------------------//
0029
0030 //! Calculated line segment between two points
0031 struct Chord
0032 {
0033     real_type length;
0034     Real3 dir;
0035 };
0036
0037 //---------------------------------------------------------------------------//
0038 // INLINE DEFINITIONS
0039 //---------------------------------------------------------------------------//
0040 /*!
0041  * Return y <- ax for a real variable
0042  */
0043 template<class T>
0044 inline CELER_FUNCTION Array<T, 3> ax(T a, Array<T, 3> const& x)
0045 {
0046     Array<T, 3> result;
0047     for (int i = 0; i < 3; ++i)
0048     {
0049         result[i] = a * x[i];
0050     }
0051     return result;
0052 }
0053
0054 //---------------------------------------------------------------------------//
0055 /*!
0056  * Calculate the direction between the source and destination.
0057  */
0058 inline CELER_FUNCTION Chord make_chord(Real3 const& src, Real3 const& dst)
0059 {
0060     Chord result;
0061     result.dir = dst - src;
0062     result.length = norm(result.dir);
0063     result.dir /= result.length;
0064     return result;
0065 }
0066
0067 //---------------------------------------------------------------------------//
0068 /*!
0069  * Whether the straight-line position is within a distance of the target.
0070  *
0071  * This is equivalent to:
0072  * \code
0073      Real3 temp = pos;
0074      axpy(distance, dir, &pos);
0075
0076      return distance(pos, target) <= tolerance;
0077  * \endcode
0078  */
0079 inline CELER_FUNCTION bool is_intercept_close(Real3 const& pos,
0080                                               Real3 const& dir,
0081                                               real_type distance,
0082                                               Real3 const& target,
0083                                               real_type tolerance)
0084 {
0085     real_type delta_sq = 0;
0086     for (int i = 0; i < 3; ++i)
0087     {
0088         delta_sq += ipow<2>(pos[i] - target[i] + distance * dir[i]);
0089     }
0090     return delta_sq <= ipow<2>(tolerance);
0091 }
0092
0093 //---------------------------------------------------------------------------//
0094 /*!
0095  * Evaluate the square of the relative stepper truncation error.
0096  *
0097  * \f$ \max(\delta_\textrm{pos}^{2}, \epsilon \delta_\textrm{mom}^{2}) \f$
0098  *
0099  * The return value is the square of \c dyerr in
0100  * \c G4MagIntegratorDriver::AccurateAdvance .
0101  */
0102 inline CELER_FUNCTION real_type rel_err_sq(OdeState const& err_state,
0103                                            real_type step,
0104                                            Real3 const& mom)
0105 {
0106     CELER_EXPECT(step > 0);
0107
0108     // Evaluate square of the position and momentum accuracy
0109     real_type errpos2 = dot_product(err_state.pos, err_state.pos);
0110     real_type errvel2 = dot_product(err_state.mom, err_state.mom);
0111
0112     // Scale position error relative to step
0113     errpos2 /= ipow<2>(step);
0114     // Scale momentum error relative to starting momentum
0115     errvel2 /= dot_product(mom, mom);
0116
0117     real_type result = max(errpos2, errvel2);
0118     CELER_ENSURE(result >= 0);
0119     return result;
0120 }
0121
0122 //---------------------------------------------------------------------------//
0123 /*!
0124  * Closest distance between the segment from beg.pos (\em A) to end.pos(\em B)
0125  * and the mid.pos (\em M):
0126  * \f[
0127  *   d = |\vec{AM}| \sin(\theta) = \frac{\vec{AM} \times \vec{AB}}{|\vec{AB}|}
0128  * \f]
0129  */
0130 inline CELER_FUNCTION real_type distance_chord(Real3 const& beg_pos,
0131                                                Real3 const& mid_pos,
0132                                                Real3 const& end_pos)
0133 {
0134     Real3 beg_mid;
0135     Real3 beg_end;
0136
0137     for (int i = 0; i < 3; ++i)
0138     {
0139         beg_mid[i] = mid_pos[i] - beg_pos[i];
0140         beg_end[i] = end_pos[i] - beg_pos[i];
0141     }
0142
0143     Real3 cross = cross_product(beg_end, beg_mid);
0144     return std::sqrt(dot_product(cross, cross) / dot_product(beg_end, beg_end));
0145 }
0146
0147 //---------------------------------------------------------------------------//
0148 }  // namespace detail
0149 }  // namespace celeritas