From 12a6aee85a51da6881635e401e0dd703baf2c5cb Mon Sep 17 00:00:00 2001 From: John Apostolakis Date: Mon, 14 Jul 2025 18:35:55 +0200 Subject: [PATCH] Additional RK integration methods --- .../AdePT/magneticfield/BogackiShampine23.h | 114 ++++++++ .../AdePT/magneticfield/DoLoMcPrinceRK34.h | 138 ++++++++++ include/AdePT/magneticfield/NystromRk4.h | 249 ++++++++++++++++++ 3 files changed, 501 insertions(+) create mode 100644 include/AdePT/magneticfield/BogackiShampine23.h create mode 100644 include/AdePT/magneticfield/DoLoMcPrinceRK34.h create mode 100644 include/AdePT/magneticfield/NystromRk4.h diff --git a/include/AdePT/magneticfield/BogackiShampine23.h b/include/AdePT/magneticfield/BogackiShampine23.h new file mode 100644 index 00000000..3d7ae800 --- /dev/null +++ b/include/AdePT/magneticfield/BogackiShampine23.h @@ -0,0 +1,114 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 10 Jan 2024 +// +// Bogacki-Shampine - 4 - 3(2) non-FSAL implementation +// +// An implementation of the embedded RK method from the paper +// [1] P. Bogacki and L. F. Shampine, +// "A 3(2) pair of Runge - Kutta formulas" +// Appl. Math. Lett., vol. 2, no. 4, pp. 321-325, Jan. 1989. +// +// Non-FSAL implementation +// Based on G4BogackiShampine23 +// Created: Somnath Banerjee, Google Summer of Code 2015, 20 May 2015 +// Supervision: John Apostolakis, CERN +// -------------------------------------------------------------------- + +#pragma once + +template +class BogackiShampine23 // : public VScalarIntegrationStepper +{ + // using ThreeVector = vecgeom::Vector3D; + +public: + static constexpr unsigned int kMethodOrder = 3; + + static __host__ __device__ + void EvaluateDerivatives(const T_Field & field, const Real_t y[], int charge, Real_t dydx[]) ; + + static __host__ __device__ + void StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t * yOut, // Output: y values at end, + Real_t * yerr, // estimated errors, + Real_t * next_dydx); // next value of dydx +}; + +template +__host__ __device__ void +BogackiShampine23:: + EvaluateDerivatives( const T_Field& field, const Real_t yIn[], int charge, Real_t dy_ds[] ) +{ + Equation_t::EvaluateDerivatives( /* const T_Field& */ field, yIn, charge, dy_ds ); +} + +// The Bogacki shampine method has the following Butcher's tableau +// +// 0 | +// 1/2|1/2 +// 3/4|0 3/4 +// 1 |2/9 1/3 4/9 +// ------------------- +// |2/9 1/3 4/9 0 +// |7/24 1/4 1/3 1/8 + +template +inline +__host__ __device__ void +BogackiShampine23:: + StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t * yOut, + Real_t * yError, + Real_t * next_dydx ) +// yIn and yOut MUST NOT be aliases for same array +{ + assert( yIn != yOut ); + + static constexpr Real_t + b21 = 0.5 , + b31 = 0., b32 = 3.0 / 4.0, + b41 = 2.0 / 9.0, b42 = 1.0 / 3.0, b43 = 4.0 / 9.0; + + static constexpr Real_t dc1 = b41 - 7.0 / 24.0, dc2 = b42 - 1.0 / 4.0, + dc3 = b43 - 1.0 / 3.0, dc4 = - 1.0 / 8.0; + + Real_t ak2[Nvar]; + { + Real_t yTemp2[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp2[i] = yIn[i] + b21 * Step * dydx[i]; + } + EvaluateDerivatives( field, yTemp2, charge, ak2); // 2nd Step + } + + Real_t ak3[Nvar]; + { + Real_t yTemp3[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp3[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]); + } + EvaluateDerivatives( field, yTemp3, charge, ak3); // 3rd Step + } + + for(unsigned int i=0;i< Nvar;i++) + { + yOut[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i] ); + } + EvaluateDerivatives( field, yOut, charge, next_dydx); // 3rd Step + + for(unsigned int i=0;i< Nvar;i++) + { + yError[i] = Step * (dc1 * dydx[i] + dc2 * ak2[i] + + dc3 * ak3[i] + dc4 * next_dydx[i]); + } +} \ No newline at end of file diff --git a/include/AdePT/magneticfield/DoLoMcPrinceRK34.h b/include/AdePT/magneticfield/DoLoMcPrinceRK34.h new file mode 100644 index 00000000..40b89ec4 --- /dev/null +++ b/include/AdePT/magneticfield/DoLoMcPrinceRK34.h @@ -0,0 +1,138 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 10 Jan 2024 +// +// Implementation of the Dormand Lockyer McGorrigan Prince non-FSAL method for AdePT. +// +// Notes: +// - Current version is restricted to Magnetic fields (see EvaluateDerivatives.) +// - It provides the next value of dy/ds in 'next_dydx' +// - It uses a large number of registers and/or stack locations - 6 derivatives + In + Out + Err +// +// Dormand-Lockyer-McGorrigan-Prince-6-3-4 non-FSAL method +// ( 6 stage, 3rd & 4th order embedded RK method ) + +// Based on G4DoLoMcPrinceRK34, +// Created by Somnath Banerjee, Google Summer of Code 2015, 7 July 2015 +// Supervision: John Apostolakis, CERN +// -------------------------------------------------------------------- + +#pragma once + +#include "RkIntegrationDriver.h" + +template +class DoLoMcPrinceRK34 // : public VScalarIntegrationStepper +{ + public: + static constexpr unsigned int kMethodOrder = 4; + + static __host__ __device__ + void EvaluateDerivatives(const T_Field & field, const Real_t y[], int charge, Real_t dydx[]) ; + + static __host__ __device__ + void StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t yOut[Nvar], // Output: y values at end, + Real_t yerr[Nvar], // estimated errors, + Real_t * next_dydx ); // next value of dydx +}; + +template +__host__ __device__ void +DoLoMcPrinceRK34:: + EvaluateDerivatives( const T_Field& field, const Real_t yIn[], int charge, Real_t dy_ds[] ) +{ + Equation_t::EvaluateDerivatives( /* const T_Field& */ field, yIn, charge, dy_ds ); +} + +template +inline +__host__ __device__ void +DoLoMcPrinceRK34:: + StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t yOut[Nvar], + Real_t yErr[Nvar], + Real_t * next_dydx ) +// yIn and yOut MUST NOT be aliases for same array +{ + // The constants from the butcher tableu + // + static constexpr Real_t + b21 = 7.0/27.0 , + b31 = 7.0/72.0 , b32 = 7.0/24.0 , + b41 = 3043.0/3528.0 , b42 = -3757.0/1176. , b43 = 1445.0/441.0, + b51 = 17617.0/11662 , b52 = -4023.0/686.0 , b53 = 9372.0/1715. , b54 = -66.0/595.0 , + b61 = 29.0/238.0 , b62 = 0.0 , b63 = 216.0/385.0 , b64 = 54.0/85.0 , b65 = -7.0/22.0 , + + dc1 = 363.0/2975.0 - b61 , + dc2 = 0.0 - b62 , + dc3 = 981.0/1750.0 - b63, + dc4 = 2709.0/4250.0 - b64 , + dc5 = -3.0/10.0 - b65 , + dc6 = -1.0/50.0 ; // end of declaration + + assert( yIn != yOut); + + // EvaluateDerivatives( field, yIn, charge, dydx) ; // 1st Step + + Real_t ak2[Nvar]; + { + Real_t yTemp2[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp2[i] = yIn[i] + b21 * Step * dydx[i]; + } + EvaluateDerivatives( field, yTemp2, charge, ak2); // 2nd Step + } + + Real_t ak3[Nvar]; + { + Real_t yTemp3[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp3[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]); + } + EvaluateDerivatives( field, yTemp3, charge, ak3); // 3rd Step + } + + Real_t ak4[Nvar]; + { + Real_t yTemp4[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp4[i] = yIn[i] + Step * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]); + } + EvaluateDerivatives( field, yTemp4, charge, ak4); // 4th Step + } + + Real_t ak5[Nvar]; + { + Real_t yTemp5[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yTemp5[i] = yIn[i] + Step * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]); + } + EvaluateDerivatives( field, yTemp5, charge, ak5); // 5th Step + } + + Real_t ak6[Nvar]; + for (unsigned int i = 0; i < Nvar; i++) { + yOut[i] = + yIn[i] + Step * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]); + } + EvaluateDerivatives( field, yOut, charge, ak6); // 6th evaluation -- only for Derivative ... ? + + for (unsigned int i = 0; i < Nvar; i++) + { + yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] + + dc5*ak5[i] + dc6*ak6[i] ) ; + } + + return ; +} + diff --git a/include/AdePT/magneticfield/NystromRk4.h b/include/AdePT/magneticfield/NystromRk4.h new file mode 100644 index 00000000..cb7a6269 --- /dev/null +++ b/include/AdePT/magneticfield/NystromRk4.h @@ -0,0 +1,249 @@ +// SPDX-FileCopyrightText: 2021 CERN +// SPDX-License-Identifier: Apache-2.0 +// +// Author: J. Apostolakis, 11 Dec 2023 +// +// Implementation of the Nystrom Rk4 Runge-Kutta integrator for AdePT. +// Derived from the Geant4 G4NystromRK4 +// +// Notes: +// - Current version is restricted to Magnetic fields (see EvaluateDerivatives.) +// - It provides the next value of dy/ds in 'next_dydx' +// - It uses a large number of registers and/or stack locations - 7 derivatives + In + Out + Err + +template +class NystromRK4 // : public VScalarIntegrationStepper +{ + // using ThreeVector = vecgeom::Vector3D; + +public: + static constexpr unsigned int kMethodOrder = 4; + // inline NystromRK4(Equation_t *EqRhs, bool verbose = false); + // NystromRK4(const NystromRK4 &) = delete; + // ~NystromRK4() {} + + static __host__ __device__ + void EvaluateDerivatives(const T_Field & field, const Real_t y[], int charge, Real_t dydx[]) ; + + static __host__ __device__ + void StepWithErrorEstimate( const T_Field & field, + const Real_t * yIn, + const Real_t * dydx, + int charge, + Real_t Step, + Real_t * yOut, // Output: y values at end, + Real_t * yerr, // estimated errors, + Real_t * next_dydx); // next value of dydx +}; + +template +__host__ __device__ void +NystromRK4:: + EvaluateDerivatives( const T_Field& field, const Real_t yIn[], int charge, Real_t dy_ds[] ) +{ +/* #ifdef VERBOSE_RHS + using geant::units::tesla; + std::cout << "NystromRK4::EvaluateDerivatives called with q= " << charge + << " at Position = " << yIn[0] << " y= " << yIn[1] << " z= " << yIn[2] + << " with Momentum = " << yIn[3] << " y= " << yIn[4] << " z= " << yIn[5] << " "; +#endif */ + + // Vector3D Bfield; + // Equation_t::EvaluateDerivativesReturnB( field, yIn, charge, dy_ds, Bfield ); + + Equation_t::EvaluateDerivatives( /* const T_Field& */ field, yIn, charge, dy_ds ); + + /********* + using copcore::units::tesla; + using std::setw; + constexpr int prec= 5; + constexpr int nf= prec+5; + int old_prec = std::cout.precision(prec); + std::cout << " DoPri5: Evaluate Derivatives - using B-field, Bx= " << Bfield.x() / tesla << " By= " << Bfield.y() / tesla << " Bz= " << Bfield.z() / tesla << " "; + std::cout << " gives Derivs dy_ds= : " + << " x = " << setw(nf) << dy_ds[0] << " y = " << setw(nf) << dy_ds[1] << " z = " << setw(nf) << dy_ds[2] + << " px= " << setw(nf) << dy_ds[3] << " py= " << setw(nf) << dy_ds[4] << " pz= " << setw(nf) << dy_ds[5] + << std::endl; + std::cout.precision(old_prec); + ********/ +} + +template +inline +__host__ __device__ void +NystromRK4:: + StepWithErrorEstimate( const T_Field & field, + const Real_t P[], // yIn + const Real_t dPdS[], // dydx + int charge, + Real_t Step, + Real_t * Po, // yOut + Real_t * Err, // yErr + Real_t * next_dydx ) +// yIn and yOut MUST NOT be aliases for same array +{ + assert( yIn != yOut ); + const Real_t perMillion = 1.0e-6; + Real_t R[4] = { P[0], P[1] , P[2], P[7] }; // x, y, z, t + Real_t A[3] = {dPdS[0], dPdS[1], dPdS[2]}; + + // m_iPoint[0]=R[0]; m_iPoint[1]=R[1]; m_iPoint[2]=R[2]; + + constexpr Real_t one_sixth= 1./6.; + const Real_t S = Step ; + const Real_t S5 = .5*Step ; + const Real_t S4 = .25*Step ; + const Real_t S6 = Step * one_sixth; // Step / 6.; + + // Ensure that the location and cached field value are correct + EvaluateDerivatives( field, R ); + + // Ensure that the momentum is set correctly. + + // - Quick check momentum magnitude (squared) against previous value + Real_t newmom2 = (P[3]*P[3]+P[4]*P[4]+P[5]*P[5]); + Real_t oldmom2 = m_mom * m_mom; + if( std::fabs(newmom2 - oldmom2) > perMillion * oldmom2 ) + { + m_mom = std::sqrt(newmom2) ; + m_imom = 1./m_mom; + m_cof = m_fEq->FCof()*m_imom; + } + +#ifdef G4DEBUG_FIELD + CheckCachedMomemtum( P, m_mom ); + CheckFieldPosition( P, m_fldPosition ); +#endif + + // Point 1 + // + Real_t K1[3] = { m_imom*dPdS[3], m_imom*dPdS[4], m_imom*dPdS[5] }; + + // Point2 + // + Real_t p[4] = { R[0]+S5*(A[0]+S4*K1[0]), + R[1]+S5*(A[1]+S4*K1[1]), + R[2]+S5*(A[2]+S4*K1[2]), + P[7] }; + EvaluateDerivatives( field, p); + + Real_t A2[3] = {A[0]+S5*K1[0], A[1]+S5*K1[1], A[2]+S5*K1[2]}; + Real_t K2[3] = {(A2[1]*m_lastField[2]-A2[2]*m_lastField[1])*m_cof, + (A2[2]*m_lastField[0]-A2[0]*m_lastField[2])*m_cof, + (A2[0]*m_lastField[1]-A2[1]*m_lastField[0])*m_cof}; + + m_mPoint[0]=p[0]; m_mPoint[1]=p[1]; m_mPoint[2]=p[2]; + + // Point 3 with the same magnetic field + // + Real_t A3[3] = {A[0]+S5*K2[0],A[1]+S5*K2[1],A[2]+S5*K2[2]}; + Real_t K3[3] = {(A3[1]*m_lastField[2]-A3[2]*m_lastField[1])*m_cof, + (A3[2]*m_lastField[0]-A3[0]*m_lastField[2])*m_cof, + (A3[0]*m_lastField[1]-A3[1]*m_lastField[0])*m_cof}; + + // Point 4 + // + p[0] = R[0]+S*(A[0]+S5*K3[0]); + p[1] = R[1]+S*(A[1]+S5*K3[1]); + p[2] = R[2]+S*(A[2]+S5*K3[2]); + + EvaluateDerivatives( field, p); + + Real_t A4[3] = {A[0]+S*K3[0],A[1]+S*K3[1],A[2]+S*K3[2]}; + Real_t K4[3] = {(A4[1]*m_lastField[2]-A4[2]*m_lastField[1])*m_cof, + (A4[2]*m_lastField[0]-A4[0]*m_lastField[2])*m_cof, + (A4[0]*m_lastField[1]-A4[1]*m_lastField[0])*m_cof}; + + // New position + // + Po[0] = P[0]+S*(A[0]+S6*(K1[0]+K2[0]+K3[0])); + Po[1] = P[1]+S*(A[1]+S6*(K1[1]+K2[1]+K3[1])); + Po[2] = P[2]+S*(A[2]+S6*(K1[2]+K2[2]+K3[2])); + + m_fPoint[0]=Po[0]; m_fPoint[1]=Po[1]; m_fPoint[2]=Po[2]; + + // New direction + // + Po[3] = A[0]+S6*(K1[0]+K4[0]+2.*(K2[0]+K3[0])); + Po[4] = A[1]+S6*(K1[1]+K4[1]+2.*(K2[1]+K3[1])); + Po[5] = A[2]+S6*(K1[2]+K4[2]+2.*(K2[2]+K3[2])); + + // Errors + // + Err[3] = S*std::fabs(K1[0]-K2[0]-K3[0]+K4[0]); + Err[4] = S*std::fabs(K1[1]-K2[1]-K3[1]+K4[1]); + Err[5] = S*std::fabs(K1[2]-K2[2]-K3[2]+K4[2]); + Err[0] = S*Err[3] ; + Err[1] = S*Err[4] ; + Err[2] = S*Err[5] ; + Err[3]*= m_mom ; + Err[4]*= m_mom ; + Err[5]*= m_mom ; + + // Normalize momentum + // + Real_t normF = m_mom/std::sqrt(Po[3]*Po[3]+Po[4]*Po[4]+Po[5]*Po[5]); + Po [3]*=normF; Po[4]*=normF; Po[5]*=normF; + + // Pass Energy, time unchanged -- time is not integrated !! + // Po[6]=P[6]; Po[7]=P[7]; +#if ENABLE_CHORD_DIST + for (unsigned int i = 0; i < Nvar; i++) { + // Store Input and Final values, for possible use in calculating chord + fLastInitialVector[i] = yIn[i]; + fLastFinalVector[i] = yOut[i]; + fInitialDyDx[i] = dydx[i]; // At initial point + } +#endif + // fLastStepLength = Step; + + // std::cout << " Exiting StepWithErrorEstimate of scalar " << std::endl; + + return; +} + +/********************************************************************** +#if ENABLE_CHORD_DIST +template +inline Real_t NystromRK4::DistChord() const +{ + // Coefficients were taken from Some Practical Runge-Kutta Formulas by Lawrence F. Shampine, page 149, c* + static constexpr Real_t hf1 = 6025192743.0 / 30085553152.0, + hf2 = 0.0, + hf3 = 51252292925.0 / 65400821598.0, + hf4 = - 2691868925.0 / 45128329728.0, + hf5 = 187940372067.0 / 1594534317056.0, + hf6 = - 1776094331.0 / 19743644256.0, + hf7 = 11237099.0 / 235043384.0; + + Real_t midVector[3]; + + for(int i = 0; i < 3; ++i) { + midVector[i] = fLastInitialVector[i] + 0.5 * fLastStepLength * + (hf1 * fInitialDyDx[i] + hf2 * ak2[i] + hf3 * ak3[i] + + hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * next_dydx[i]); + } + Real_t distChord; + ThreeVector initialPoint, finalPoint, midPoint; + + initialPoint = ThreeVector(fLastInitialVector[0], fLastInitialVector[1], fLastInitialVector[2]); + finalPoint = ThreeVector(fLastFinalVector[0], fLastFinalVector[1], fLastFinalVector[2]); + midPoint = ThreeVector(midVector[0], midVector[1], midVector[2]); + + // Use stored values of Initial and Endpoint + new Midpoint to evaluate + // distance of Chord + distChord = GULineSection::Distline(midPoint, initialPoint, finalPoint); + + return distChord; +} +#endif +**********************************************************************/ + +// template +// inline void NystromRK4::PrintField(const char *label, const Real_t y[Nvar], +// const vecgeom::Vector3D &Bfield) const + + +// template +// inline void NystromRK4::PrintDyDx(const char *label, const Real_t dydx[Nvar], +// const Real_t y[Nvar]) const