G4InterpolationDriver.icc

Go to the documentation of this file.

//
// ********************************************************************
// * 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.          *
// ********************************************************************
//
// G4InterpolationDriver inline implementation
//
// Created: D.Sorokin, 2018
// --------------------------------------------------------------------
 
#include "G4FieldUtils.hh"
#include "G4LineSection.hh"
#include "G4Exception.hh"
 
#include <algorithm>
 
 
template <class T>
G4InterpolationDriver<T>::
G4InterpolationDriver ( G4double hminimum, T* pStepper,
                        G4int numComponents, G4int statisticsVerbose )
    : G4RKIntegrationDriver<T>(pStepper),
      fMinimumStep(hminimum),
      fVerboseLevel(statisticsVerbose)
{
    if (numComponents != Base::GetStepper()->GetNumberOfVariables())
    {
        std::ostringstream message;
        message << "Driver's number of integrated components "
                << numComponents
                << " != Stepper's number of components "
                << pStepper->GetNumberOfVariables();
        G4Exception("G4InterpolationDriver","GeomField0002",
                    FatalException, message);
    }
 
    for (G4int i = 0; i < Base::GetMaxNoSteps(); ++i) 
    {
        fSteppers.push_back({
            std::unique_ptr<T>(new T(pStepper->GetSpecificEquation(),   // Interpolating stepper must have this!
                                     pStepper->GetNumberOfVariables()) ),
            DBL_MAX, -DBL_MAX, 0.0
        });
    }
 
    fLastStepper = fSteppers.end();
}
 
template <class T>
G4InterpolationDriver<T>::~G4InterpolationDriver()
{
#ifdef G4VERBOSE
    if (fVerboseLevel > 0)
    {
        G4cout << "G4ChordFinder statistics report: \n"
               << "  No trials: " << fTotalNoTrials
               << "  No Calls: "  << fNoCalls
               << "  Max-trial: " <<  fmaxTrials
               << G4endl; 
    }
#endif
}
 
template <class T>
void G4InterpolationDriver<T>::OnStartTracking()
{
    fChordStepEstimate = DBL_MAX;
    fhnext = DBL_MAX;
    fTotalStepsForTrack = 0;
}
 
template <class T>
void G4InterpolationDriver<T>::OnComputeStep()
{ 
    fKeepLastStepper = false;
    fFirstStep = true;
    fLastStepper = fSteppers.end();
}
 
template <class T>
void G4InterpolationDriver<T>::SetVerboseLevel(G4int level)
{
    fVerboseLevel = level;
}
 
template <class T>
G4int G4InterpolationDriver<T>::GetVerboseLevel() const
{
    return fVerboseLevel;
}
 
template <class T>
void G4InterpolationDriver<T>::
Interpolate(G4double curveLength, field_utils::State& y) const
{
    if (fLastStepper == fSteppers.end()) 
    {
        std::ostringstream message;
        message << "LOGICK ERROR: fLastStepper == end";
        G4Exception("G4InterpolationDriver::Interpolate()", 
                    "GeomField1001", FatalException, message);
        return;
    }
 
    ConstStepperIterator end = fLastStepper + 1;
 
    auto it = std::lower_bound(fSteppers.cbegin(), end, curveLength, 
                               [](const InterpStepper& stepper, G4double value)
                               {
                                 return stepper.end < value;
                               }
                              );
    if (it == end)
    {
        if (curveLength - fLastStepper->end > CLHEP::perMillion)
        {
            std::ostringstream message;
            message << "curveLength = " << curveLength << " > "
                    << fLastStepper->end;
            G4Exception("G4InterpolationDriver::Interpolate()", 
                        "GeomField1001", JustWarning, message);
        }
 
        return fLastStepper->stepper->Interpolate(1, y);
    }
 
    if (curveLength < it->begin) 
    {
        if (it->begin - curveLength > CLHEP::perMillion)
        {
            std::ostringstream message;
            message << "curveLength = " << curveLength << " < " << it->begin;
            G4Exception("G4InterpolationDriver::Interpolate()", 
                        "GeomField1001", JustWarning, message);
        }
 
        return it->stepper->Interpolate(0, y);
    }
 
    return InterpolateImpl(curveLength, it, y);
}
 
template <class T>
void G4InterpolationDriver<T>::InterpolateImpl(G4double curveLength, 
                                               ConstStepperIterator it, 
                                               field_utils::State& y) const
{
    const G4double tau = (curveLength - it->begin) * it->inverseLength;
    return it->stepper->Interpolate(field_utils::clamp(tau, 0., 1.), y);  
}
 
template <class T>
G4double G4InterpolationDriver<T>::DistChord(const field_utils::State& yBegin, 
                                             G4double curveLengthBegin, 
                                             const field_utils::State& yEnd, 
                                             G4double curveLengthEnd) const
{
    // optimization check if it worth
    //
    if (curveLengthBegin == fLastStepper->begin && 
        curveLengthEnd == fLastStepper->end) 
    {
        return fLastStepper->stepper->DistChord();
    }
 
    const G4double curveLengthMid = 0.5 * (curveLengthBegin + curveLengthEnd);
    field_utils::State yMid;
    
    Interpolate(curveLengthMid, yMid);
 
    return G4LineSection::Distline(
        field_utils::makeVector(yMid, field_utils::Value3D::Position),
        field_utils::makeVector(yBegin, field_utils::Value3D::Position),
        field_utils::makeVector(yEnd, field_utils::Value3D::Position)
    );
}
 
template <class T>
G4double G4InterpolationDriver<T>::AdvanceChordLimited(G4FieldTrack& track,
                                                       G4double hstep,
                                                       G4double epsStep,
                                                       G4double chordDistance)
{
    ++fTotalStepsForTrack;
 
    const G4double curveLengthBegin = track.GetCurveLength();
    const G4double hend = std::min(hstep, fChordStepEstimate);
    G4double hdid = 0.0;
    auto it = fSteppers.begin();
    G4double dChordStep = 0.0;
 
    field_utils::State yBegin, y;
    track.DumpToArray(yBegin);
    track.DumpToArray(y);
 
    if (fFirstStep) 
    {
        Base::GetEquationOfMotion()->RightHandSide(y, fdydx);
        fFirstStep = false;
    }
 
    if (fKeepLastStepper) 
    {
        std::swap(*fSteppers.begin(), *fLastStepper);
        it = fSteppers.begin(); //new begin, update iterator
        fLastStepper = it;
        hdid = it->end - curveLengthBegin;
        if (hdid > hend)
        {
            hdid = hend;
            InterpolateImpl(curveLengthBegin + hdid, it, y);
        }
        else
        {
            field_utils::copy(y, it->stepper->GetYOut());
        }
 
        dChordStep = DistChord(yBegin, curveLengthBegin, y,
                               curveLengthBegin + hdid);
 
        ++it;
    }
 
    // accurate advance & check chord distance
    G4double h = fhnext;
    for (; hdid<hend && dChordStep<chordDistance && it!=fSteppers.end(); ++it)
    {
        h = std::min(h, hstep - hdid);
 
        // make one step
        hdid += OneGoodStep(it, y, fdydx, h, epsStep, curveLengthBegin + hdid);
 
        // update last stepper
        fLastStepper = it;
 
        // estimate chord distance
        dChordStep = std::max( dChordStep,
            DistChord(yBegin, curveLengthBegin, y, curveLengthBegin + hdid) );
    }
 
    // Now, either
    //   - full integration ( hdid >= hend )
    //   - estimated chord has exceeded limit 'chordDistance'
    //   - reached maximum number of steps (from number of steppers.)
    
    // update step estimation
    if (h > fMinimumStep)
    {
        fhnext = h;
    }
 
    // CheckState();
 
    // update chord step estimate
    //
    hdid = FindNextChord(yBegin, curveLengthBegin, y, curveLengthBegin + hdid, 
                         dChordStep, chordDistance);
 
    const G4double curveLengthEnd = curveLengthBegin + hdid;
    fKeepLastStepper = fLastStepper->end - curveLengthEnd > CLHEP::perMillion;
    track.LoadFromArray(y, fLastStepper->stepper->GetNumberOfVariables());
    track.SetCurveLength(curveLengthBegin + hdid);
 
    return hdid;
}
 
template <class T>
G4double G4InterpolationDriver<T>::
FindNextChord(const field_utils::State& yBegin,
                    G4double curveLengthBegin,
                    field_utils::State& yEnd,
                    G4double curveLengthEnd,
                    G4double dChord,
                    G4double chordDistance)
{
    G4double hstep = curveLengthEnd - curveLengthBegin;
    G4double curveLength = curveLengthEnd;
 
    G4int i = 1;
    for (; i < fMaxTrials && dChord > chordDistance
                          && curveLength > fLastStepper->begin; ++i) 
    {
        // crop step size
        hstep = CalcChordStep(hstep, dChord, chordDistance);
 
        // hstep should be in the last stepper
        hstep = std::max(hstep, fLastStepper->begin - curveLengthBegin);
        curveLength = curveLengthBegin + hstep;
 
        // use fLastStepper!
        InterpolateImpl(curveLength, fLastStepper, yEnd);
 
        // update chord distance
        dChord = DistChord(yBegin, curveLengthBegin, yEnd, curveLength);
    }
 
    // dChord may be zero
    //
    if (dChord > 0.0)
    {
        fChordStepEstimate = hstep * std::sqrt(chordDistance / dChord);
    }
 
    if (i == fMaxTrials) 
    {
        G4Exception("G4InterpolationDriver::FindNextChord()", 
                    "GeomField1001", JustWarning, "cannot converge");
    }
 
    AccumulateStatistics(i);
 
    return hstep;
}
 
// Is called to estimate the next step size, even for successful steps,
// in order to predict an accurate 'chord-sensitive' first step
// which is likely to assist in more performant 'stepping'.
//
template <class T>
G4double G4InterpolationDriver<T>::CalcChordStep(G4double stepTrialOld, 
                                                 G4double dChordStep,
                                                 G4double chordDistance)  
{
    const G4double chordStepEstimate = stepTrialOld
                                     * std::sqrt(chordDistance / dChordStep);
    G4double stepTrial = fFractionNextEstimate * chordStepEstimate;
 
    if (stepTrial <= 0.001 * stepTrialOld)
    {
        if (dChordStep > 1000.0 * chordDistance)
        {
            stepTrial = stepTrialOld * 0.03;   
        }
        else
        {
            if (dChordStep > 100. * chordDistance)
            {
                stepTrial = stepTrialOld * 0.1;   
            }
            else   // Try halving the length until dChordStep OK
            {
                stepTrial = stepTrialOld * 0.5;   
            }
        }
    }
    else if (stepTrial > 1000.0 * stepTrialOld)
    {
        stepTrial = 1000.0 * stepTrialOld;
    }
 
    if (stepTrial == 0.0)
    {
        stepTrial = 0.000001;
    }
 
    // A more sophisticated chord-finder could figure out a better
    // stepTrial, from dChordStep and the required d_geometry
    //   e.g.
    //      Calculate R, r_helix (eg at orig point)
    //      if( stepTrial < 2 pi  R )
    //          stepTrial = R arc_cos( 1 - chordDistance / r_helix )
    //      else    
    //          ??
 
    return stepTrial;
}
 
template <class T>
G4bool G4InterpolationDriver<T>::
AccurateAdvance(G4FieldTrack& track, G4double hstep,
                G4double /*eps*/,
                G4double /*hinitial*/
   )
{
    if (hstep == 0.0)
    {
        std::ostringstream message;
        message << "Proposed step is zero; hstep = " << hstep << " !";
        G4Exception("G4InterpolationDriver::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("G4InterpolationDriver::AccurateAdvance()", 
                    "GeomField0003", EventMustBeAborted, message);
        return false;
    }
 
    const G4double curveLength = track.GetCurveLength();
    const G4double curveLengthEnd = curveLength + hstep;
 
    field_utils::State y;
    Interpolate(curveLengthEnd, y);
 
    track.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
    track.SetCurveLength(curveLengthEnd);
 
    return true;
}
 
// Driver for one Runge-Kutta Step with monitoring of local truncation error
// to ensure accuracy and adjust stepsize. Input are dependent variable
// array y[0,...,5] and its derivative dydx[0,...,5] at the
// starting value of the independent variable x . Also input are stepsize
// to be attempted htry, and the required accuracy eps. On output y and x
// are replaced by their new values, hdid is the stepsize that was actually
// accomplished, and hnext is the estimated next stepsize. 
// This is similar to the function rkqs from the book:
// Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
// Edition, by William H. Press, Saul A. Teukolsky, William T.
// Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
// 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
//
template <class T>
G4double G4InterpolationDriver<T>::OneGoodStep(StepperIterator it,
                                               field_utils::State& y,
                                               field_utils::State& dydx,
                                               G4double& hstep,
                                               G4double epsStep,
                                               G4double curveLength)
 
{
    G4double error2 = DBL_MAX;
    field_utils::State yerr, ytemp, dydxtemp;
    G4double h = hstep;
 
    G4int i = 0;
    for (; i < fMaxTrials; ++i)
    {
        it->stepper->Stepper(y, dydx, h, ytemp, yerr, dydxtemp); 
        error2 = field_utils::relativeError2(y, yerr, h, epsStep);
 
        if (error2 <= 1.0)
        {
            hstep = std::max(Base::GrowStepSize2(h, error2), fMinimumStep);
            break; 
        }
 
        // don't control error for small steps
        if (h <= fMinimumStep)
        {
            hstep = fMinimumStep;
            break;
        }
 
        h = std::max(Base::ShrinkStepSize2(h, error2), fMinimumStep);
    }
 
    if (i == fMaxTrials) 
    {
        G4Exception("G4InterpolationDriver::OneGoodStep()", 
                    "GeomField1001", JustWarning, "cannot converge");
        hstep = std::max(Base::ShrinkStepSize2(h, error2), fMinimumStep);
    }
 
    // set interpolation inverval
    it->begin = curveLength;
    it->end = curveLength + h;
    it->inverseLength = 1. / h;
 
    // setup interpolation
    it->stepper->SetupInterpolation();
 
    field_utils::copy(dydx, dydxtemp);
    field_utils::copy(y, ytemp);
 
    return h;
}
 
template <class T>
void G4InterpolationDriver<T>::PrintState() const
{
    using namespace field_utils;
    State prevEnd, currBegin;
    auto prev = fSteppers.begin();
 
    G4cout << "====== curr state ========" << G4endl;
    for (auto i = fSteppers.begin(); i <= fLastStepper; ++i)
    {
        i->stepper->Interpolate(0, currBegin);
 
        G4cout << "cl_begin: " <<i->begin << " "
               << "cl_end: " << i->end << " ";
 
        if (prev != i) {
            prev->stepper->Interpolate(1, prevEnd);
            auto prevPos = makeVector(prevEnd, Value3D::Position);
            auto currPos = makeVector(currBegin, Value3D::Position);
            G4cout << "diff_begin: " << (prevPos - currPos).mag();
        }
                   
        G4cout << G4endl;
        prev = i;
    }
 
    const G4double clBegin = fSteppers.begin()->begin;
    const G4double clEnd = fLastStepper->end;
    const G4double hstep = (clEnd - clBegin) / 10.;
    State yBegin, yCurr;
    Interpolate(0, yBegin);
    for (G4double cl = clBegin; cl <= clEnd + 1e-12; cl += hstep)
    {
        Interpolate(cl, yCurr);
        auto d = DistChord(yBegin, clBegin, yCurr, cl);
        G4cout << "cl: " << cl << " chord_distance: " << d << G4endl;
    }
 
    G4cout << "==========================" << G4endl;
}
 
 
template <class T>
void G4InterpolationDriver<T>::CheckState() const
{
    G4int smallSteps = 0;
    for (auto i = fSteppers.begin(); i <= fLastStepper; ++i) 
    {
        G4double stepLength = i->end - i->begin;
        if (stepLength < fMinimumStep)
        {
            ++smallSteps;
        }
    }
 
    if (smallSteps > 1) 
    {
        std::ostringstream message;
        message << "====== curr state ========\n";
        for (auto i = fSteppers.begin(); i <= fLastStepper; ++i)
        {
            message << "cl_begin: " <<i->begin << " "
                    << "cl_end: " << i->end << "\n";
        }
 
        G4Exception("G4InterpolationDriver::CheckState()", 
                    "GeomField0003", FatalException, message);
    }
}
 
template <class T>
void G4InterpolationDriver<T>::AccumulateStatistics(G4int noTrials) 
{
    fTotalNoTrials += noTrials; 
    ++fNoCalls; 
      
    if (noTrials > fmaxTrials) 
    { 
        fmaxTrials = noTrials; 
    }
}
 
template <class T>
void G4InterpolationDriver<T>::StreamInfo( std::ostream& os ) const
{
    os << "State of G4InterpolationDriver: " << std::endl;
    os << "--Base state (G4RKIntegrationDriver): " << std::endl;
    Base::StreamInfo( os );    
    os << "  fMinimumStep   =      " << fMinimumStep  << 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;
 
    os << "  Max num of Trials   = " << fMaxTrials << std::endl;
    os << "  Fract Next Estimate = " << fFractionNextEstimate << std::endl;
    os << "  Smallest Curve Fract= " << fSmallestCurveFraction << std::endl;
    
    os << "  VerboseLevel        = " << fVerboseLevel << std::endl;
    os << "  KeepLastStepper     = " << fKeepLastStepper << std::endl;
}