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