EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-01-18 09:54:49
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 corecel/math/ArrayUtils.hh 0007 //! \brief Math functions using celeritas::Array 0008 //---------------------------------------------------------------------------// 0009 #pragma once 0010 0011 #include <cmath> 0012 0013 #include "corecel/Assert.hh" 0014 #include "corecel/Types.hh" 0015 #include "corecel/cont/Array.hh" 0016 0017 #include "Algorithms.hh" 0018 #include "ArraySoftUnit.hh" 0019 #include "SoftEqual.hh" 0020 0021 #include "detail/ArrayUtilsImpl.hh" 0022 0023 namespace celeritas 0024 { 0025 //---------------------------------------------------------------------------// 0026 // Perform y <- ax + y 0027 template<class T, size_type N> 0028 inline CELER_FUNCTION void axpy(T a, Array<T, N> const& x, Array<T, N>* y); 0029 0030 //---------------------------------------------------------------------------// 0031 // Calculate product of two vectors 0032 template<class T, size_type N> 0033 [[nodiscard]] inline CELER_FUNCTION T dot_product(Array<T, N> const& x, 0034 Array<T, N> const& y); 0035 0036 //---------------------------------------------------------------------------// 0037 // Calculate product of two vectors 0038 template<class T> 0039 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0040 cross_product(Array<T, 3> const& x, Array<T, 3> const& y); 0041 0042 //---------------------------------------------------------------------------// 0043 // Calculate the Euclidian (2) norm of a vector 0044 template<class T, size_type N> 0045 [[nodiscard]] inline CELER_FUNCTION T norm(Array<T, N> const& vec); 0046 0047 //---------------------------------------------------------------------------// 0048 // Construct a vector with unit magnitude 0049 template<class T, size_type N> 0050 [[nodiscard]] inline CELER_FUNCTION Array<T, N> 0051 make_unit_vector(Array<T, N> const& v); 0052 0053 //---------------------------------------------------------------------------// 0054 // Calculate the Euclidian (2) distance between two points 0055 template<class T, size_type N> 0056 [[nodiscard]] inline CELER_FUNCTION T distance(Array<T, N> const& x, 0057 Array<T, N> const& y); 0058 0059 //---------------------------------------------------------------------------// 0060 // Calculate a cartesian unit vector from spherical coordinates 0061 template<class T> 0062 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0063 from_spherical(T costheta, T phi); 0064 0065 //---------------------------------------------------------------------------// 0066 // Rotate the direction 'dir' according to the reference rotation axis 'rot' 0067 template<class T> 0068 [[nodiscard]] inline CELER_FUNCTION Array<T, 3> 0069 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot); 0070 0071 //---------------------------------------------------------------------------// 0072 // INLINE DEFINITIONS 0073 //---------------------------------------------------------------------------// 0074 /*! 0075 * Increment a vector by another vector multiplied by a scalar. 0076 */ 0077 template<class T, size_type N> 0078 CELER_FUNCTION void axpy(T a, Array<T, N> const& x, Array<T, N>* y) 0079 { 0080 CELER_EXPECT(y); 0081 for (size_type i = 0; i != N; ++i) 0082 { 0083 (*y)[i] = fma(a, x[i], (*y)[i]); 0084 } 0085 } 0086 0087 //---------------------------------------------------------------------------// 0088 /*! 0089 * Dot product of two vectors. 0090 */ 0091 template<class T, size_type N> 0092 CELER_FUNCTION T dot_product(Array<T, N> const& x, Array<T, N> const& y) 0093 { 0094 T result{}; 0095 for (size_type i = 0; i != N; ++i) 0096 { 0097 result = fma(x[i], y[i], result); 0098 } 0099 return result; 0100 } 0101 0102 //---------------------------------------------------------------------------// 0103 /*! 0104 * Cross product of two space vectors. 0105 */ 0106 template<class T> 0107 CELER_FUNCTION Array<T, 3> 0108 cross_product(Array<T, 3> const& x, Array<T, 3> const& y) 0109 { 0110 return {x[1] * y[2] - x[2] * y[1], 0111 x[2] * y[0] - x[0] * y[2], 0112 x[0] * y[1] - x[1] * y[0]}; 0113 } 0114 0115 //---------------------------------------------------------------------------// 0116 /*! 0117 * Calculate the Euclidian (2) norm of a vector. 0118 */ 0119 template<class T, size_type N> 0120 CELER_FUNCTION T norm(Array<T, N> const& v) 0121 { 0122 return std::sqrt(dot_product(v, v)); 0123 } 0124 0125 //---------------------------------------------------------------------------// 0126 /*! 0127 * Construct a unit vector. 0128 * 0129 * Unit vectors have an Euclidian norm magnitude of 1. 0130 */ 0131 template<class T, size_type N> 0132 CELER_FUNCTION Array<T, N> make_unit_vector(Array<T, N> const& v) 0133 { 0134 Array<T, N> result{v}; 0135 T const scale_factor = 1 / norm(result); 0136 for (auto& el : result) 0137 { 0138 el *= scale_factor; 0139 } 0140 return result; 0141 } 0142 0143 //---------------------------------------------------------------------------// 0144 /*! 0145 * Calculate the Euclidian (2) distance between two points. 0146 */ 0147 template<class T, size_type N> 0148 CELER_FUNCTION T distance(Array<T, N> const& x, Array<T, N> const& y) 0149 { 0150 T dist_sq = 0; 0151 for (size_type i = 0; i != N; ++i) 0152 { 0153 dist_sq += ipow<2>(y[i] - x[i]); 0154 } 0155 return std::sqrt(dist_sq); 0156 } 0157 0158 //---------------------------------------------------------------------------// 0159 /*! 0160 * Calculate a Cartesian vector from spherical coordinates. 0161 * 0162 * Theta is the angle between the Z axis and the outgoing vector, and phi is 0163 * the angle between the x axis and the projection of the vector onto the x-y 0164 * plane. 0165 */ 0166 template<class T> 0167 inline CELER_FUNCTION Array<T, 3> from_spherical(T costheta, T phi) 0168 { 0169 CELER_EXPECT(costheta >= -1 && costheta <= 1); 0170 0171 T const sintheta = std::sqrt(1 - costheta * costheta); 0172 return {sintheta * std::cos(phi), sintheta * std::sin(phi), costheta}; 0173 } 0174 0175 //---------------------------------------------------------------------------// 0176 /*! 0177 * Rotate the direction about the given Z-based scatter direction. 0178 * 0179 * The equivalent to calling the Shift transport code's \code 0180 void cartesian_vector_transform( 0181 double costheta, 0182 double phi, 0183 Vector_View vector); 0184 * \endcode 0185 * is the call 0186 * \code 0187 vector = rotate(from_spherical(costheta, phi), vector); 0188 * \endcode 0189 * 0190 * This code effectively decomposes the given rotation vector \c rot into two 0191 * sequential transform matrices, one with an angle \em theta about the \em y 0192 * axis and one about \em phi rotating around the \em z axis. These two angles 0193 * are the spherical coordinate transform of the given \c rot cartesian 0194 * direction vector. 0195 * 0196 * There is some extra code in here to deal with loss of precision when the 0197 * incident direction is along the \em z axis. As \c rot approaches \em z, the 0198 * azimuthal angle \em phi must be calculated carefully from both the x and y 0199 * components of the vector, not independently. If \c rot actually equals \em z 0200 * then the azimuthal angle is completely indeterminate so we arbitrarily 0201 * choose \c phi = 0. 0202 * 0203 * This function is often used for calculating exiting scattering angles. In 0204 * that case, \c dir is the exiting angle from the scattering calculation, and 0205 * \c rot is the original direction of the particle. The direction vectors are 0206 * defined as 0207 * \f[ 0208 \Omega = \sin\theta\cos\phi\mathbf{i} 0209 + \sin\theta\sin\phi\mathbf{j} 0210 + \cos\theta\mathbf{k} \,. 0211 * \f] 0212 */ 0213 template<class T> 0214 inline CELER_FUNCTION Array<T, 3> 0215 rotate(Array<T, 3> const& dir, Array<T, 3> const& rot) 0216 { 0217 CELER_EXPECT(is_soft_unit_vector(dir)); 0218 CELER_EXPECT(is_soft_unit_vector(rot)); 0219 0220 // Direction enumeration 0221 enum 0222 { 0223 X = 0, 0224 Y = 1, 0225 Z = 2 0226 }; 0227 0228 // Transform direction vector into theta, phi so we can use it as a 0229 // rotation matrix 0230 T sintheta = std::sqrt(1 - ipow<2>(rot[Z])); 0231 T cosphi; 0232 T sinphi; 0233 0234 if (sintheta >= detail::RealVecTraits<T>::min_accurate_sintheta()) 0235 { 0236 // Typical case: far enough from z axis to assume the X and Y 0237 // components have a hypotenuse of 1 within epsilon tolerance 0238 T const inv_sintheta = 1 / (sintheta); 0239 cosphi = rot[X] * inv_sintheta; 0240 sinphi = rot[Y] * inv_sintheta; 0241 } 0242 else if (sintheta > 0) 0243 { 0244 // Avoid catastrophic roundoff error by normalizing x/y components 0245 cosphi = rot[X] / std::sqrt(ipow<2>(rot[X]) + ipow<2>(rot[Y])); 0246 sinphi = std::sqrt(1 - ipow<2>(cosphi)); 0247 } 0248 else 0249 { 0250 // NaN or 0: choose an arbitrary azimuthal angle for the incident dir 0251 cosphi = 1; 0252 sinphi = 0; 0253 } 0254 0255 Array<T, 3> result 0256 = {(rot[Z] * dir[X] + sintheta * dir[Z]) * cosphi - sinphi * dir[Y], 0257 (rot[Z] * dir[X] + sintheta * dir[Z]) * sinphi + cosphi * dir[Y], 0258 -sintheta * dir[X] + rot[Z] * dir[Z]}; 0259 0260 // Always normalize to prevent roundoff error from propagating 0261 return make_unit_vector(result); 0262 } 0263 0264 //---------------------------------------------------------------------------// 0265 } // namespace celeritas
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|