include/Geant4/G4BulirschStoerDriver.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4BulirschStoer driver inline methods implementation
0027
0028 // Author: Dmitry Sorokin, Google Summer of Code 2016
0029 // Supervision: John Apostolakis, CERN
0030 // --------------------------------------------------------------------
0031
0032 #include <cassert>
0033
0034 #include "G4LineSection.hh"
0035 #include "G4FieldUtils.hh"
0036
0037 G4IntegrationDriver<G4BulirschStoer>::
0038 G4IntegrationDriver( G4double hminimum, G4BulirschStoer* stepper,
0039                      G4int numberOfComponents, G4int statisticsVerbosity )
0040   : fMinimumStep(hminimum),
0041     fVerbosity(statisticsVerbosity),
0042     fMidpointMethod(stepper->GetEquationOfMotion(),
0043                     stepper->GetNumberOfVariables()),
0044     bulirschStoer(stepper),
0045     interval_sequence{2,4}
0046 {
0047     assert(stepper->GetNumberOfVariables() == numberOfComponents);
0048
0049     if(stepper->GetNumberOfVariables() != numberOfComponents)
0050     {
0051       std::ostringstream msg;
0052       msg << "Disagreement in number of variables = "
0053           << stepper->GetNumberOfVariables()
0054           << " vs no of components = " << numberOfComponents;
0055       G4Exception("G4IntegrationDriver<G4BulirschStoer> Constructor:",
0056                   "GeomField1001", FatalException, msg);
0057     }
0058 }
0059
0060 G4bool G4IntegrationDriver<G4BulirschStoer>::
0061 AccurateAdvance( G4FieldTrack& track, G4double hstep,
0062                  G4double eps, G4double hinitial)
0063 {
0064     G4int fNoTotalSteps = 0;
0065     G4int fMaxNoSteps = 10000;
0066     G4double fNoBadSteps = 0.0;
0067     G4double fSmallestFraction = 1.0e-12;
0068
0069     // Driver with adaptive stepsize control. Integrate starting
0070     // values at y_current over hstep x2 with accuracy eps.
0071     // On output ystart is replaced by values at the end of the integration
0072     // interval. RightHandSide is the right-hand side of ODE system.
0073     // The source is similar to odeint routine from NRC p.721-722 .
0074
0075     //  Ensure that hstep > 0
0076     if(hstep == 0)
0077     {
0078       std::ostringstream message;
0079       message << "Proposed step is zero; hstep = " << hstep << " !";
0080       G4Exception("G4IntegrationDriver<G4BulirschStoer>::AccurateAdvance()",
0081                   "GeomField1001", JustWarning, message);
0082
0083       return true;
0084     }
0085     if(hstep < 0)
0086     {
0087         std::ostringstream message;
0088         message << "Invalid run condition." << G4endl
0089                 << "Proposed step is negative; hstep = "
0090                 << hstep << "." << G4endl
0091                 << "Requested step cannot be negative! Aborting event.";
0092         G4Exception("G4IntegrationDriver<G4BulirschStoer>::AccurateAdvance()",
0093                     "GeomField0003", EventMustBeAborted, message);
0094
0095         return false;
0096     }
0097
0098     // init first step size
0099     //
0100     G4double h;
0101     if ( (hinitial > 0) && (hinitial < hstep)
0102       && (hinitial > perMillion * hstep) )
0103     {
0104        h = hinitial;
0105     }
0106     else  //  Initial Step size "h" defaults to the full interval
0107     {
0108        h = hstep;
0109     }
0110
0111     // integration variables
0112     //
0113     track.DumpToArray(yCurrent);
0114
0115     // copy non-integration variables to out array
0116     //
0117     std::memcpy(yOut + GetNumberOfVarialbles(),
0118                 yCurrent + GetNumberOfVarialbles(),
0119          sizeof(G4double)*(G4FieldTrack::ncompSVEC-GetNumberOfVarialbles()));
0120
0121     G4double startCurveLength = track.GetCurveLength();
0122     G4double curveLength = startCurveLength;
0123     G4double endCurveLength = startCurveLength + hstep;
0124
0125     // loop variables
0126     //
0127     G4int nstp = 1, no_warnings = 0;
0128     G4double hnext, hdid;
0129
0130     G4bool succeeded = true, lastStepSucceeded;
0131
0132     G4int  noFullIntegr = 0, noSmallIntegr = 0 ;
0133     static G4ThreadLocal G4int  noGoodSteps = 0 ;  // Bad = chord > curve-len
0134
0135     G4bool lastStep = false;
0136
0137     // BulirschStoer->reset();
0138
0139     G4FieldTrack yFldTrk(track);
0140
0141     do
0142     {
0143         G4ThreeVector StartPos(yCurrent[0], yCurrent[1], yCurrent[2]);
0144         GetEquationOfMotion()->RightHandSide(yCurrent, dydxCurrent);
0145         fNoTotalSteps++;
0146
0147         // Perform the Integration
0148         //
0149         if(h == 0)
0150         {
0151           G4Exception("G4IntegrationDriver<G4BulirschStoer>::AccurateAdvance()",
0152                       "GeomField0003", FatalException,
0153                       "Integration Step became Zero!");
0154         }
0155         else if(h > fMinimumStep)
0156         {
0157             // step size if Ok
0158             //
0159             OneGoodStep(yCurrent,dydxCurrent,curveLength,h,eps,hdid,hnext);
0160             lastStepSucceeded = (hdid == h);
0161         }
0162         else  // for small steps call QuickAdvance for speed
0163         {
0164             G4double dchord_step, dyerr, dyerr_len;  // What to do with these ?
0165             yFldTrk.LoadFromArray(yCurrent, G4FieldTrack::ncompSVEC);
0166             yFldTrk.SetCurveLength(curveLength);
0167
0168             QuickAdvance(yFldTrk, dydxCurrent, h, dchord_step, dyerr_len);
0169
0170             yFldTrk.DumpToArray(yCurrent);
0171
0172
0173             dyerr = dyerr_len / h;
0174             hdid = h;
0175             curveLength += hdid;
0176
0177             // Compute suggested new step
0178             // hnext = ComputeNewStepSize(dyerr/eps, h);
0179             hnext = h;
0180
0181             // hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
0182             lastStepSucceeded = (dyerr <= eps);
0183         }
0184
0185         lastStepSucceeded ? ++noFullIntegr : ++noSmallIntegr;
0186
0187         G4ThreeVector EndPos(yCurrent[0], yCurrent[1], yCurrent[2]);
0188
0189         // Check the endpoint
0190         //
0191         G4double endPointDist = (EndPos - StartPos).mag();
0192         if (endPointDist >= hdid*(1. + perMillion))
0193         {
0194             ++fNoBadSteps;
0195
0196             // Issue a warning only for gross differences -
0197             // we understand how small difference occur.
0198             //
0199             if (endPointDist >= hdid*(1.+perThousand))
0200             {
0201                 ++no_warnings;
0202             }
0203         }
0204         else
0205         {
0206             ++noGoodSteps;
0207         }
0208
0209         //  Avoid numerous small last steps
0210         //
0211         if((h < eps * hstep) || (h < fSmallestFraction * startCurveLength))
0212         {
0213             // No more integration -- the next step will not happen
0214             //
0215             lastStep = true;
0216         }
0217         else
0218         {
0219             // Check the proposed next stepsize
0220             //
0221             if(std::fabs(hnext) < fMinimumStep)
0222             {
0223               // Make sure that the next step is at least Hmin
0224               //
0225               h = fMinimumStep;
0226             }
0227             else
0228             {
0229               h = hnext;
0230             }
0231
0232             //  Ensure that the next step does not overshoot
0233             //
0234             if (curveLength + h > endCurveLength)
0235             {
0236               h = endCurveLength - curveLength;
0237             }
0238
0239             if (h == 0)
0240             {
0241               // Cannot progress - accept this as last step - by default
0242               //
0243               lastStep = true;
0244             }
0245         }
0246     } while (((nstp++) <= fMaxNoSteps)
0247            && (curveLength < endCurveLength) && (!lastStep));
0248        //
0249        // Have we reached the end ?
0250        // --> a better test might be x-x2 > an_epsilon
0251
0252     succeeded = (curveLength >= endCurveLength);
0253       // If it was a "forced" last step
0254
0255     // Copy integrated vars to output array
0256     //
0257     std::memcpy(yOut, yCurrent, sizeof(G4double) * GetNumberOfVarialbles());
0258
0259     // upload new state
0260     track.LoadFromArray(yOut, G4FieldTrack::ncompSVEC);
0261     track.SetCurveLength(curveLength);
0262
0263     if(nstp > fMaxNoSteps)
0264     {
0265       ++no_warnings;
0266       succeeded = false;
0267     }
0268
0269     return succeeded;
0270 }
0271
0272 G4bool G4IntegrationDriver<G4BulirschStoer>::
0273 QuickAdvance( G4FieldTrack& track, const G4double dydx[],
0274               G4double hstep, G4double& missDist, G4double& dyerr)
0275 {
0276     const auto nvar = fMidpointMethod.GetNumberOfVariables();
0277
0278     track.DumpToArray(yIn);
0279     const G4double curveLength = track.GetCurveLength();
0280
0281     fMidpointMethod.SetSteps(interval_sequence[0]);
0282     fMidpointMethod.DoStep(yIn, dydx, yOut, hstep, yMid, derivs[0]);
0283
0284     fMidpointMethod.SetSteps(interval_sequence[1]);
0285     fMidpointMethod.DoStep(yIn, dydx, yOut2, hstep, yMid2, derivs[1]);
0286
0287     // extrapolation
0288     //
0289     static const G4double coeff =
0290         1. / (sqr(static_cast<G4double>(interval_sequence[1]) /
0291                      static_cast<G4double>(interval_sequence[0])) - 1.);
0292
0293     for (G4int i = 0; i < nvar; ++i)
0294     {
0295         yOut[i] = yOut2[i] + (yOut2[i] - yOut[i]) * coeff;
0296         yMid[i] = yMid2[i] + (yMid2[i] - yMid[i]) * coeff;
0297     }
0298
0299     // calculate chord length
0300     //
0301     const auto mid = field_utils::makeVector(yMid,
0302                      field_utils::Value3D::Position);
0303     const auto in  = field_utils::makeVector(yIn,
0304                      field_utils::Value3D::Position);
0305     const auto out = field_utils::makeVector(yOut,
0306                      field_utils::Value3D::Position);
0307
0308     missDist = G4LineSection::Distline(mid, in, out);
0309
0310     // calculate error
0311     //
0312     for (G4int i = 0; i < nvar; ++i)
0313     {
0314         yError[i] = yOut[i] - yOut2[i];
0315     }
0316
0317     dyerr = field_utils::absoluteError(yOut, yError, hstep);
0318
0319     // copy non-integrated variables to output array
0320     //
0321     std::memcpy(yOut + nvar, yIn + nvar,
0322                 sizeof(G4double) * (G4FieldTrack::ncompSVEC - nvar));
0323
0324     // set new state
0325     //
0326     track.LoadFromArray(yOut, G4FieldTrack::ncompSVEC);
0327     track.SetCurveLength(curveLength + hstep);
0328
0329     return true;
0330 }
0331
0332 void G4IntegrationDriver<G4BulirschStoer>::
0333 OneGoodStep( G4double y[], const G4double dydx[], G4double& curveLength,
0334              G4double htry, G4double eps, G4double& hdid, G4double& hnext)
0335 {
0336     hnext = htry;
0337     G4double curveLengthBegin = curveLength;
0338
0339     // set maximum allowed error
0340     //
0341     bulirschStoer->set_max_relative_error(eps);
0342
0343     while (true)
0344     {
0345         auto res = bulirschStoer->try_step(y, dydx, curveLength, yOut, hnext);
0346         if (res == G4BulirschStoer::step_result::success)
0347         {
0348             break;
0349         }
0350     }
0351
0352     std::memcpy(y, yOut, sizeof(G4double) * GetNumberOfVarialbles());
0353     hdid = curveLength - curveLengthBegin;
0354 }
0355
0356 void G4IntegrationDriver<G4BulirschStoer>::
0357 GetDerivatives( const G4FieldTrack& track, G4double dydx[]) const
0358 {
0359     G4double y[G4FieldTrack::ncompSVEC];
0360     track.DumpToArray(y);
0361     GetEquationOfMotion()->RightHandSide(y, dydx);
0362 }
0363
0364 void G4IntegrationDriver<G4BulirschStoer>::
0365 GetDerivatives( const G4FieldTrack& track, G4double dydx[],
0366                 G4double field[]) const
0367 {
0368     G4double y[G4FieldTrack::ncompSVEC];
0369     track.DumpToArray(y);
0370     GetEquationOfMotion()->EvaluateRhsReturnB(y, dydx, field);
0371 }
0372
0373 void G4IntegrationDriver<G4BulirschStoer>::SetVerboseLevel(G4int level)
0374 {
0375     fVerbosity = level;
0376 }
0377
0378 G4int G4IntegrationDriver<G4BulirschStoer>::GetVerboseLevel() const
0379 {
0380     return fVerbosity;
0381 }
0382
0383 G4double G4IntegrationDriver<G4BulirschStoer>::
0384 ComputeNewStepSize( G4double /* errMaxNorm*/, G4double  hstepCurrent)
0385 {
0386     return hstepCurrent;
0387 }
0388
0389 G4EquationOfMotion* G4IntegrationDriver<G4BulirschStoer>::GetEquationOfMotion()
0390 {
0391     assert(bulirschStoer->GetEquationOfMotion() ==
0392            fMidpointMethod.GetEquationOfMotion());
0393
0394     return bulirschStoer->GetEquationOfMotion();
0395 }
0396
0397 const G4EquationOfMotion*
0398 G4IntegrationDriver<G4BulirschStoer>::GetEquationOfMotion() const
0399 {
0400     return const_cast<G4IntegrationDriver<G4BulirschStoer>*>(this)->
0401             GetEquationOfMotion();
0402 }
0403
0404 void G4IntegrationDriver<G4BulirschStoer>::
0405 SetEquationOfMotion( G4EquationOfMotion* equation )
0406 {
0407     bulirschStoer->SetEquationOfMotion(equation);
0408     fMidpointMethod.SetEquationOfMotion(equation);
0409 }
0410
0411 G4int G4IntegrationDriver<G4BulirschStoer>::GetNumberOfVarialbles() const
0412 {
0413     assert(bulirschStoer->GetNumberOfVariables() ==
0414            fMidpointMethod.GetNumberOfVariables());
0415
0416     return bulirschStoer->GetNumberOfVariables();
0417 }
0418
0419 const G4MagIntegratorStepper*
0420 G4IntegrationDriver<G4BulirschStoer>::GetStepper() const
0421 {
0422     return nullptr;
0423 }
0424
0425 G4MagIntegratorStepper*
0426 G4IntegrationDriver<G4BulirschStoer>::GetStepper()
0427 {
0428     return nullptr;
0429 }
0430
0431 void
0432 G4IntegrationDriver<G4BulirschStoer>::StreamInfo( std::ostream& os ) const
0433 {
0434   os << "State of G4IntegrationDriver<G4BulirschStoer>: " << std::endl;
0435   os << "   Method is implemented, but gives no information. " << std::endl;
0436 }