Gabriele Cosmo
238 lines
7.1 KiB
Plaintext
238 lines
7.1 KiB
Plaintext
//
|
|
// ********************************************************************
|
|
// * 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. *
|
|
// ********************************************************************
|
|
//
|
|
// G4RKIntegrationDriver inline implementation
|
|
//
|
|
// Author: Dmitry Sorokin (CERN, Google Summer of Code 2017), 20.10.2017
|
|
// --------------------------------------------------------------------
|
|
|
|
#include "globals.hh"
|
|
|
|
template <class T>
|
|
G4RKIntegrationDriver<T>::G4RKIntegrationDriver(T* pStepper)
|
|
{
|
|
RenewStepperAndAdjustImpl(pStepper);
|
|
fMaxStepBase = 250;
|
|
fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
|
|
}
|
|
|
|
// Step failed; compute the size of retrial Step.
|
|
//
|
|
template <class T> G4double G4RKIntegrationDriver<T>::
|
|
ShrinkStepSize(G4double h, G4double error) const
|
|
{
|
|
if (error > errorConstraintShrink)
|
|
{
|
|
return max_stepping_decrease * h;
|
|
}
|
|
return GetSafety() * h * std::pow(error, GetPshrnk());
|
|
}
|
|
|
|
template <class T> G4double G4RKIntegrationDriver<T>::
|
|
ShrinkStepSize2(G4double h, G4double error2) const
|
|
{
|
|
if (error2 > errorConstraintShrink * errorConstraintShrink)
|
|
{
|
|
return max_stepping_decrease * h;
|
|
}
|
|
return GetSafety() * h * std::pow(error2, 0.5 * GetPshrnk());
|
|
}
|
|
|
|
// Compute size of next Step
|
|
//
|
|
template<class T> G4double G4RKIntegrationDriver<T>::
|
|
GrowStepSize(G4double h, G4double error) const
|
|
{
|
|
if (error < errorConstraintGrow)
|
|
{
|
|
return max_stepping_increase * h;
|
|
}
|
|
return GetSafety() * h * std::pow(error, GetPgrow());
|
|
}
|
|
|
|
template<class T> G4double G4RKIntegrationDriver<T>::
|
|
GrowStepSize2(G4double h, G4double error2) const
|
|
{
|
|
if (error2 < errorConstraintGrow * errorConstraintGrow)
|
|
{
|
|
return max_stepping_increase * h;
|
|
}
|
|
return GetSafety() * h * std::pow(error2, 0.5 * GetPgrow());
|
|
}
|
|
|
|
template <class T>
|
|
G4double G4RKIntegrationDriver<T>::
|
|
ComputeNewStepSize( G4double errMaxNorm, // max error (normalised)
|
|
G4double hstepCurrent ) // current step size
|
|
{
|
|
if (errMaxNorm > 1)
|
|
{
|
|
return ShrinkStepSize(hstepCurrent, errMaxNorm);
|
|
}
|
|
if (errMaxNorm >= 0)
|
|
{
|
|
return GrowStepSize(hstepCurrent, errMaxNorm);
|
|
}
|
|
|
|
G4Exception("G4RKIntegrationDriver::ComputeNewStepSize", "GeomField0003",
|
|
FatalException, "Error is negative!");
|
|
|
|
return max_stepping_increase * hstepCurrent;
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::
|
|
GetDerivatives( const G4FieldTrack& track, G4double dydx[] ) const
|
|
{
|
|
G4double y[G4FieldTrack::ncompSVEC];
|
|
track.DumpToArray(y);
|
|
pIntStepper->RightHandSide(y, dydx);
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::GetDerivatives(const G4FieldTrack& track,
|
|
G4double dydx[],
|
|
G4double field[]) const
|
|
{
|
|
G4double y[G4FieldTrack::ncompSVEC];
|
|
track.DumpToArray(y);
|
|
pIntStepper->RightHandSide(y, dydx, field);
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::UpdateErrorConstraints()
|
|
{
|
|
errorConstraintShrink = std::pow(
|
|
max_stepping_decrease / GetSafety(), 1. / GetPshrnk());
|
|
|
|
errorConstraintGrow = std::pow(
|
|
max_stepping_increase / GetSafety(), 1. / GetPgrow());
|
|
}
|
|
|
|
template <class T>
|
|
inline G4double G4RKIntegrationDriver<T>::GetSafety() const
|
|
{
|
|
return safety;
|
|
}
|
|
|
|
template <class T>
|
|
inline G4double G4RKIntegrationDriver<T>::GetPshrnk() const
|
|
{
|
|
return pshrnk;
|
|
}
|
|
|
|
template <class T>
|
|
G4double G4RKIntegrationDriver<T>::GetPgrow() const
|
|
{
|
|
return pgrow;
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::ReSetParameters(G4double new_safety)
|
|
{
|
|
safety = new_safety;
|
|
pshrnk = -1.0 / pIntStepper->IntegratorOrder();
|
|
pgrow = -1.0 / (1.0 + pIntStepper->IntegratorOrder());
|
|
UpdateErrorConstraints();
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::SetSafety(G4double val)
|
|
{
|
|
safety = val;
|
|
UpdateErrorConstraints();
|
|
}
|
|
|
|
template <class T> void G4RKIntegrationDriver<T>::
|
|
RenewStepperAndAdjust(G4MagIntegratorStepper* stepper)
|
|
{
|
|
T* ourStepper = dynamic_cast<T*>(stepper);
|
|
if ( ourStepper )
|
|
{
|
|
RenewStepperAndAdjustImpl( ourStepper );
|
|
}
|
|
else
|
|
{
|
|
G4Exception("G4RKIntegrationDriver::RenewStepperAndAdjust()",
|
|
"GeomField0002", FatalException,
|
|
"The type of the stepper provided is incorrect for this templated driver");
|
|
}
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::RenewStepperAndAdjustImpl(T* stepper)
|
|
{
|
|
pIntStepper = stepper;
|
|
ReSetParameters();
|
|
}
|
|
|
|
template <class T>
|
|
const T* G4RKIntegrationDriver<T>::GetStepper() const
|
|
{
|
|
return pIntStepper;
|
|
}
|
|
|
|
template <class T>
|
|
T* G4RKIntegrationDriver<T>::GetStepper()
|
|
{
|
|
return pIntStepper;
|
|
}
|
|
|
|
template <class T>
|
|
G4int G4RKIntegrationDriver<T>::GetMaxNoSteps() const
|
|
{
|
|
return fMaxNoSteps;
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::SetMaxNoSteps(G4int val)
|
|
{
|
|
fMaxNoSteps = val;
|
|
}
|
|
|
|
template <class T>
|
|
G4EquationOfMotion* G4RKIntegrationDriver<T>::GetEquationOfMotion()
|
|
{
|
|
return pIntStepper->GetEquationOfMotion();
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::SetEquationOfMotion(G4EquationOfMotion* equation)
|
|
{
|
|
pIntStepper->SetEquationOfMotion(equation);
|
|
}
|
|
|
|
template <class T>
|
|
void G4RKIntegrationDriver<T>::StreamInfo( std::ostream& os ) const
|
|
{
|
|
os << "State of G4RKIntegrationDriver: " << std::endl;
|
|
os << " Max number of Steps = " << fMaxNoSteps << std::endl;
|
|
os << " Safety factor = " << safety << std::endl;
|
|
os << " Power - shrink = " << pshrnk << std::endl;
|
|
os << " Power - grow = " << pgrow << std::endl;
|
|
os << " threshold - shrink = " << errorConstraintShrink << std::endl;
|
|
os << " threshold - grow = " << errorConstraintGrow << std::endl;
|
|
}
|