From fb4b85b2fd544d164b948d2e963f094aebfa86cb Mon Sep 17 00:00:00 2001 From: divyanshtiwari237 <200020049@iitb.ac.in> Date: Sun, 11 Sep 2022 07:04:13 +0530 Subject: [PATCH] Implemented BorisSDC Implemented the borisSDC method which involves the boris method along with collocation techniques solved via spectral deferred corrections --- .../magneticfield/include/G4BorisDriverSDC.hh | 151 ++++++++++ .../include/G4BorisDriverSDC.icc | 257 ++++++++++++++++ .../magneticfield/include/G4BorisSDC.hh | 107 +++++++ .../magneticfield/include/G4BorisSDC.icc | 46 +++ source/geometry/magneticfield/sources.cmake | 5 + .../geometry/magneticfield/src/G4BorisSDC.cc | 280 ++++++++++++++++++ 6 files changed, 846 insertions(+) create mode 100644 source/geometry/magneticfield/include/G4BorisDriverSDC.hh create mode 100644 source/geometry/magneticfield/include/G4BorisDriverSDC.icc create mode 100644 source/geometry/magneticfield/include/G4BorisSDC.hh create mode 100644 source/geometry/magneticfield/include/G4BorisSDC.icc create mode 100644 source/geometry/magneticfield/src/G4BorisSDC.cc diff --git a/source/geometry/magneticfield/include/G4BorisDriverSDC.hh b/source/geometry/magneticfield/include/G4BorisDriverSDC.hh new file mode 100644 index 00000000000..235d048fb0b --- /dev/null +++ b/source/geometry/magneticfield/include/G4BorisDriverSDC.hh @@ -0,0 +1,151 @@ +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// G4BorisSDC +// +// Class description: +// +// G4BorisDriverSDC is a driver class using +// borisSDC method to integrate the equation of motion. + +// Author: Divyansh Tiwari, Google Summer of Code 2022 +// Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun +// -------------------------------------------------------------------- +#ifndef G4BORIS_DRIVER_SDC_HH +#define G4BORIS_DRIVER_SDC_HH + +#include "G4VIntegrationDriver.hh" +#include "G4BorisSDC.hh" +#include "G4ChordFinderDelegate.hh" + + +class G4BorisDriverSDC: + public G4VIntegrationDriver, + public G4ChordFinderDelegate +{ + public: + + inline G4BorisDriverSDC( G4double hminimum, + G4BorisSDC* Boris, + G4int numberOfComponents = 6, + G4int statisticsVerbosity = 1); + + inline ~G4BorisDriverSDC() = default; + + inline G4BorisDriverSDC(const G4BorisDriverSDC&) = delete; + inline G4BorisDriverSDC& operator=(const G4BorisDriverSDC&) = delete; + + inline virtual G4double AdvanceChordLimited(G4FieldTrack& track, + G4double hstep, + G4double eps, + G4double chordDistance) override + { + return ChordFinderDelegate:: + AdvanceChordLimitedImpl(track, hstep, eps, chordDistance); + } + + inline virtual void OnStartTracking() override + { + ChordFinderDelegate::ResetStepEstimate(); + } + + inline virtual void OnComputeStep() override {}; + + inline virtual G4bool DoesReIntegrate() const override { return false; } /// ???? + + inline virtual G4bool AccurateAdvance( G4FieldTrack& track, + G4double stepLen, + G4double eps, + G4double beginStep = 0) override; + + inline virtual G4bool QuickAdvance( G4FieldTrack& y_val, + const G4double dydx[], + G4double hstep, + G4double& missDist, + G4double& dyerr) override; + + // inline void OneGoodStep(G4FieldTrack& track, + // G4double y[], + // const G4double dydx[], + // G4double& curveLength, + // G4double htry, + // G4double eps, + // G4double& hdid, + // G4double& hnext); + + inline virtual void GetDerivatives( const G4FieldTrack& track, + G4double dydx[]) const override; + + inline virtual void GetDerivatives( const G4FieldTrack& track, + G4double dydx[], + G4double field[]) const override; + + inline virtual void SetVerboseLevel(G4int level) override; + inline virtual G4int GetVerboseLevel() const override; + + inline virtual G4double ComputeNewStepSize( + G4double errMaxNorm, // normalised error + G4double hstepCurrent) override; // current step size + + inline virtual G4EquationOfMotion* GetEquationOfMotion() override; + inline const G4EquationOfMotion* GetEquationOfMotion() const; + inline virtual void SetEquationOfMotion(G4EquationOfMotion* equation) override; + + inline virtual const G4MagIntegratorStepper* GetStepper() const override; + inline virtual G4MagIntegratorStepper* GetStepper() override; + + inline virtual void StreamInfo( std::ostream& os ) const override; + // Write out the parameters / state of the driver + + private: + + inline G4int GetNumberOfVarialbles() const; + + G4double fMinimumStep; + G4double fVerbosity; + + G4BorisSDC* boris; + + G4double yIn[G4FieldTrack::ncompSVEC], + yMid[G4FieldTrack::ncompSVEC], + yMid2[G4FieldTrack::ncompSVEC], + yOut[G4FieldTrack::ncompSVEC], + yOut2[G4FieldTrack::ncompSVEC], + yError[G4FieldTrack::ncompSVEC]; + + + G4double dydxCurrent[G4FieldTrack::ncompSVEC]; + G4double yCurrent[G4FieldTrack::ncompSVEC]; + + G4double derivs[2][6][G4FieldTrack::ncompSVEC]; + + //const G4int interval_sequence[2]; + + using ChordFinderDelegate = + G4ChordFinderDelegate; +}; + +#include "G4BorisDriverSDC.icc" + +#endif diff --git a/source/geometry/magneticfield/include/G4BorisDriverSDC.icc b/source/geometry/magneticfield/include/G4BorisDriverSDC.icc new file mode 100644 index 00000000000..a548bafa722 --- /dev/null +++ b/source/geometry/magneticfield/include/G4BorisDriverSDC.icc @@ -0,0 +1,257 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// G4BorisDriverSDC inline methods implementation + +// Author: Divyansh Tiwari, Google Summer of Code 2022 +// Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun +// -------------------------------------------------------------------- + +#include + +#include "G4LineSection.hh" +#include "G4FieldUtils.hh" +#include "G4SystemOfUnits.hh" + +G4BorisDriverSDC:: +G4BorisDriverSDC( G4double hminimum, G4BorisSDC* Boris, + G4int numberOfComponents, G4int statisticsVerbosity ) + : fMinimumStep(hminimum), + fVerbosity(statisticsVerbosity), + boris(Boris) + + +{ + + assert(boris->GetNumberOfVariables() == numberOfComponents); + + if(boris->GetNumberOfVariables() != numberOfComponents) + { + std::ostringstream msg; + msg << "Disagreement in number of variables = " + << boris->GetNumberOfVariables() + << " vs no of components = " << numberOfComponents; + G4Exception("G4BorisDriverSDC Constructor:", + "GeomField1001", FatalException, msg); + } +} + +G4bool G4BorisDriverSDC:: +AccurateAdvance( G4FieldTrack& track, G4double hstep, + G4double eps, G4double hinitial) + { + + + // Driver with adaptive stepsize control. Integrate starting + // values at y_current over hstep x2 with accuracy eps. + // On output ystart is replaced by values at the end of the integration + // interval. RightHandSide is the right-hand side of ODE system. + // The source is similar to odeint routine from NRC p.721-722 . + + // Ensure that hstep > 0 + if(hstep == 0) + { + std::ostringstream message; + message << "Proposed step is zero; hstep = " << hstep << " !"; + G4Exception("G4BorisDriverSDC::AccurateAdvance()", + "GeomField1001", JustWarning, message); + + return true; + } + if(hstep < 0) + { + std::ostringstream message; + message << "Invalid run condition." << G4endl + << "Proposed step is negative; hstep = " + << hstep << "." << G4endl + << "Requested step cannot be negative! Aborting event."; + G4Exception("G4BorisDriverSDC::AccurateAdvance()", + "GeomField0003", EventMustBeAborted, message); + + return false; + } + + // init first step size + // + G4double h = hstep; + // G4cout<<" BorisDriver::AccurateAdvance: Step Request "<GetNumberOfVariables(); + // G4cout<<" here "<DoStep(restMass, charge, yIn, yOut, hstep); + boris->DoStep(restMass, charge, yIn, yMid, hstep*0.5); + + // G4cout<GetEquationOfMotion(); + + return eq; +} + +void G4BorisDriverSDC:: +SetEquationOfMotion( G4EquationOfMotion* equation ) +{ + boris->SetEquationOfMotion(equation); + +} + +G4int G4BorisDriverSDC::GetNumberOfVarialbles() const +{ + return boris->GetNumberOfVariables(); +} + +const G4MagIntegratorStepper* +G4BorisDriverSDC::GetStepper() const +{ + return nullptr; +} + +G4MagIntegratorStepper* +G4BorisDriverSDC::GetStepper() +{ + return nullptr; +} + +void +G4BorisDriverSDC::StreamInfo( std::ostream& os ) const +{ + os << "State of G4BorisDriverSDC: " << std::endl; + os << " Method is implemented, but gives no information. " << std::endl; +} diff --git a/source/geometry/magneticfield/include/G4BorisSDC.hh b/source/geometry/magneticfield/include/G4BorisSDC.hh new file mode 100644 index 00000000000..38288ec0bdb --- /dev/null +++ b/source/geometry/magneticfield/include/G4BorisSDC.hh @@ -0,0 +1,107 @@ +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// G4BorisSDC +// +// Class description: +// +// Implimentation of the boris algorithm for advancing +// charged particles in an elctromagnetic field. + +// Author: Divyansh Tiwari, Google Summer of Code 2022 +// Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun +// -------------------------------------------------------------------- +#ifndef G4BORIS_SDC_HH +#define G4BORIS_SDC_HH + +#include "G4Types.hh" +#include "G4EquationOfMotion.hh" +#include "G4FieldTrack.hh" +#include "G4PhysicalConstants.hh" +#include"G4SystemOfUnits.hh" + +class G4BorisSDC +{ + public: + + G4BorisSDC() = default; + G4BorisSDC( G4EquationOfMotion* equation, + G4int nvar = 6); + ~G4BorisSDC() = default; + + G4ThreeVector GetLorentzForce(const G4ThreeVector position, const G4ThreeVector velocity) const; + + void DoStep(const G4double restMass, const G4double charge, const G4double yIn[], + G4double yOut[], G4double hstep); + + void UpdatePosition(G4int k, G4int m) ; + + void UpdateVelocity( G4int k, G4int m) ; + + + + + + inline void SetEquationOfMotion(G4EquationOfMotion* equation); + inline G4EquationOfMotion* GetEquationOfMotion(); + + inline G4int GetNumberOfVariables() const; + + private: + + void copy(G4double dst[], const G4double src[]) const; + + private: + + G4EquationOfMotion* fEquation = nullptr; + G4int fnvar = 0; + static constexpr G4double c_l = CLHEP::c_light/CLHEP::m*CLHEP::second; + G4double alpha; // charge/mass ratio (SI) + G4double mass_si; + G4double restMass_c2; + G4double charge_si; + static constexpr G4int M = 3; // no. of nodes(3) for the order 2M-2 = 4 + static constexpr G4int K = 4; // no. of iterations + static constexpr double sqrt15 = 3.872983346207417; + static constexpr G4double nodes[M + 1] ={0, 0.5 - sqrt15/10.0, 0.5, 0.5 + sqrt15/10.0 }; + G4double delta_t_m[M + 1]; + static constexpr G4double Q[M+1][M+1] ={0, 0, 0, 0, + 0, 5.0/36, 2.0/9 - sqrt15/15, 5.0/36 - sqrt15/30, + 0, 5.0/36 + sqrt15/24, 2.0/9, 5.0/36 - sqrt15/24, + 0, 5.0/36 + sqrt15/30, 2.0/9 + sqrt15/15, 5.0/36 }; + G4double S[M+1][M+1]; + // Initialize the starting values of velocity and position at various nodes and iterations + G4ThreeVector Velocity[M+1][K+1]; + G4ThreeVector Position[M+1][K+1]; + + + + + + + +}; + +#include "G4BorisSDC.icc" +#endif \ No newline at end of file diff --git a/source/geometry/magneticfield/include/G4BorisSDC.icc b/source/geometry/magneticfield/include/G4BorisSDC.icc new file mode 100644 index 00000000000..71ee1bbcfc0 --- /dev/null +++ b/source/geometry/magneticfield/include/G4BorisSDC.icc @@ -0,0 +1,46 @@ +// +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// G4BorisSDC inline methods implementation +// +// Author: Divyansh Tiwari, Google Summer of Code 2022 +// Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun +// -------------------------------------------------------------------- + + +inline void G4BorisSDC::SetEquationOfMotion(G4EquationOfMotion* eq) +{ + fEquation = eq; +} + +inline G4EquationOfMotion* G4BorisSDC::GetEquationOfMotion() +{ + return fEquation; +} + +inline G4int G4BorisSDC::GetNumberOfVariables() const +{ + return fnvar; +} \ No newline at end of file diff --git a/source/geometry/magneticfield/sources.cmake b/source/geometry/magneticfield/sources.cmake index a611a79273c..89e6cda4494 100644 --- a/source/geometry/magneticfield/sources.cmake +++ b/source/geometry/magneticfield/sources.cmake @@ -6,6 +6,10 @@ geant4_add_module(G4magneticfield G4BFieldIntegrationDriver.hh G4BogackiShampine23.hh G4BogackiShampine45.hh + G4BorisSDC.hh + G4BorisSDC.icc + G4BorisDriverSDC.hh + G4BorisDriverSDC.icc G4BulirschStoer.hh G4BulirschStoer.icc G4BulirschStoerDriver.hh @@ -115,6 +119,7 @@ geant4_add_module(G4magneticfield G4BFieldIntegrationDriver.cc G4BogackiShampine23.cc G4BogackiShampine45.cc + G4BorisSDC.cc G4BulirschStoer.cc G4CachedMagneticField.cc G4CashKarpRKF45.cc diff --git a/source/geometry/magneticfield/src/G4BorisSDC.cc b/source/geometry/magneticfield/src/G4BorisSDC.cc new file mode 100644 index 00000000000..c4af8c29916 --- /dev/null +++ b/source/geometry/magneticfield/src/G4BorisSDC.cc @@ -0,0 +1,280 @@ +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// G4BorisSDC implementation +// +// Author: Divyansh Tiwari, Google Summer of Code 2022 +// Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun +// -------------------------------------------------------------------- + +#include "G4BorisSDC.hh" +#include "G4FieldUtils.hh" +#include"G4SystemOfUnits.hh" +#include "globals.hh" +#include "G4PhysicalConstants.hh" +using namespace field_utils; + +G4BorisSDC::G4BorisSDC( G4EquationOfMotion* equation, + G4int nvar ) + : fEquation(equation), fnvar(nvar) +{ + if (nvar <= 0) + { + G4Exception("G4BorisSDC::G4BorisSDC()", + "GeomField0002", FatalException, + "Invalid number of variables; must be greater than zero!"); + } +} + + + +G4ThreeVector G4BorisSDC::GetLorentzForce(const G4ThreeVector position, const G4ThreeVector velocity) const +{ + // initialize variables + G4double dydx[G4FieldTrack::ncompSVEC]; + G4double y[G4FieldTrack::ncompSVEC]; + + //convert velocity into momentum + G4double v_mag = velocity.mag(); + G4ThreeVector v_dir = velocity/v_mag; + G4double momen_mag = (restMass_c2*v_mag)/(std::sqrt(c_l*c_l - v_mag*v_mag)); + G4ThreeVector momen = momen_mag*v_dir; + + // obtain the y-array + for(int i = 0; i <3; i++) + { + y[i] = position[i]; + y[i+2] = momen[i]; + } + + //Calculate field + G4double fieldValue[6] ={0,0,0,0,0,0}; + fEquation->EvaluateRhsReturnB(y, dydx, fieldValue); + + //initialize E and B + G4ThreeVector B; + G4ThreeVector E; + for( G4int i = 0; i < 3; i++) + { + E[i] = fieldValue[i+3]/CLHEP::volt*CLHEP::meter;// FIXME - Check Units + B[i] = fieldValue[i]/CLHEP::tesla; + } + + G4ThreeVector force = alpha*(E + velocity.cross(B)); + + // G4cout<EvaluateRhsReturnB(y, dydx, fieldValue); + + //initialize E and B + for( G4int i = 0; i < 3; i++) + { + E[i] = fieldValue[i+3]/CLHEP::volt*CLHEP::meter;// FIXME - Check Units + B[i] = fieldValue[i]/CLHEP::tesla; + } + + + // Calculating Velocity via the rotation scheme defined for a normal boris algorithm + G4ThreeVector v_minus = Velocity[m_i-1][k] + delta_t_m[m_i]*(alpha*E + c_k)/2; + G4double gamma; // how to find gamma + G4ThreeVector t = alpha*B*delta_t_m[m_i]/2; + G4double t_l = t[0]*t[0] + t[1]*t[1]+ t[2]*t[2]; + G4ThreeVector s_i = 2*t/(1 + t_l); + G4ThreeVector v_plus = v_minus + (v_minus + v_minus.cross(t)).cross(s_i); + G4ThreeVector v = v_plus + delta_t_m[m_i]*(alpha*E + c_k)/2; + //G4cout<<"c_k : "<