Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:59:01 

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 // G4QSStepper
0027 //
0028 // QSS Integrator Stepper
0029 
0030 // Authors: Lucio Santi, Rodrigo Castro - 2018-2021.
0031 // --------------------------------------------------------------------
0032 #ifndef QSS_Stepper_HH
0033 #define QSS_Stepper_HH
0034 
0035 #include "G4FieldTrack.hh"
0036 #include "G4FieldUtils.hh"
0037 #include "G4LineSection.hh"
0038 #include "G4MagIntegratorStepper.hh"
0039 #include "G4QSS2.hh"
0040 #include "G4QSS3.hh"
0041 #include "G4QSSDriver.hh"
0042 #include "G4QSSMessenger.hh"
0043 #include "G4VIntegrationDriver.hh"
0044 #include "G4qss_misc.hh"
0045 
0046 #include <cmath>
0047 #include <cassert>
0048 
0049 // Maximum allowed number of QSS substeps per integration step
0050 #define QSS_MAX_SUBSTEPS 1000
0051 
0052 template <class QSS>
0053 class G4QSStepper : public G4MagIntegratorStepper
0054 {
0055   public:
0056 
0057     G4QSStepper(G4EquationOfMotion* EqRhs,
0058                 G4int numberOfVariables = 6,
0059                 G4bool primary = true);
0060    ~G4QSStepper() override;
0061 
0062     void Stepper(const G4double y[],
0063                  const G4double dydx[],
0064                        G4double h,
0065                        G4double yout[],
0066                        G4double yerr[]) override;
0067 
0068     void Stepper(const G4double yInput[],
0069                  const G4double dydx[],
0070                        G4double hstep,
0071                        G4double yOutput[],
0072                        G4double yError[],
0073                        G4double dydxOutput[]);
0074 
0075     // For calculating the output at the tau fraction of Step
0076     //
0077     inline void SetupInterpolation() {}
0078     inline void Interpolate(G4double tau, G4double yOut[]);
0079 
0080     G4double DistChord() const override;
0081 
0082     G4int IntegratorOrder() const override { return method->order(); }
0083 
0084     void reset(const G4FieldTrack* track);
0085 
0086     void SetPrecision(G4double dq_rel, G4double dq_min);
0087       // precision parameters for QSS method
0088 
0089     static G4QSStepper<G4QSS2>* build_QSS2(G4EquationOfMotion* EqRhs,
0090                                            G4int numberOfVariables = 6,
0091                                            G4bool primary = true);
0092 
0093     static G4QSStepper<G4QSS3>* build_QSS3(G4EquationOfMotion* EqRhs,
0094                                            G4int numberOfVariables = 6,
0095                                            G4bool primary = true);
0096 
0097     inline G4EquationOfMotion* GetSpecificEquation() { return GetEquationOfMotion(); }
0098 
0099     inline const field_utils::State& GetYOut() const { return fyOut; }
0100 
0101     inline G4double GetLastStepLength() { return fLastStepLength; }
0102 
0103   private:
0104 
0105     G4QSStepper(QSS* method,
0106                 G4EquationOfMotion* EqRhs,
0107                 G4int numberOfVariables = 6,
0108                 G4bool primary = true);
0109 
0110     void initialize_data_structs();
0111     static QSS_simulator build_simulator();
0112 
0113     inline void update_field();
0114     inline void save_substep(G4double time, G4double length);
0115 
0116     inline void realloc_substeps();
0117     inline void get_state_from_poly(G4double* x, G4double* tx,
0118                                     G4double time, G4double* state);
0119 
0120     inline void recompute_derivatives(int index);
0121     inline void update_time();
0122 
0123     inline G4double get_coeff() { return fCoeff_local; }
0124 
0125     inline void set_coeff(G4double coeff) { fCoeff_local = coeff; }
0126 
0127     inline void set_charge(G4double q)
0128     {
0129       f_charge_c2 = q * cLight_local * cLight_local;  // 89875.5178737;
0130     }
0131 
0132     inline G4double get_qc2() { return f_charge_c2; }
0133 
0134     inline void set_mg() { fMassGamma = f_mass * fGamma2; }
0135 
0136     inline void set_gamma2(G4double gamma2) { fGamma2 = gamma2; }
0137     inline void set_velocity(G4double v) { fVelocity = v; }
0138 
0139     inline void velocity_to_momentum(G4double* state);
0140 
0141     inline void set_gamma(G4double p_sq)
0142     {
0143       set_gamma2(std::sqrt(p_sq / (f_mass * f_mass) + 1));
0144       set_mg();
0145       set_coeff(get_qc2() / fMassGamma);
0146     }
0147 
0148   private:
0149 
0150     QSS_simulator simulator;
0151     QSS* method;
0152 
0153     // State
0154     //
0155     G4double fLastStepLength;
0156     field_utils::State fyIn, fyOut;
0157 
0158     G4double f_mass;
0159     static constexpr G4double cLight_local = 299.792458;  // should use CLHEP
0160     G4double f_charge_c2;
0161     G4double fMassGamma;
0162     G4double fGamma2;
0163     G4double fCoeff_local;
0164     G4double fVelocity;
0165 };
0166 
0167 using G4QSStepper_QSS2 = G4QSStepper<G4QSS2>;
0168 using G4QSStepper_QSS3 = G4QSStepper<G4QSS3>;
0169 
0170 template <class QSS>
0171 inline G4QSStepper<QSS>::G4QSStepper(QSS* qss, G4EquationOfMotion* EqRhs,
0172                               G4int noIntegrationVariables, G4bool)
0173   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
0174     simulator(qss->getSimulator()),
0175     method(qss)
0176 {
0177   SetIsQSS(true);  //  Replaces virtual method IsQSS
0178   fLastStepLength = -1.0;
0179 
0180   f_mass = 0;
0181   f_charge_c2 = 0;
0182   fMassGamma = 0;
0183   fGamma2 = 0;
0184   fCoeff_local = 0;
0185   fVelocity = 0;
0186 
0187   this->initialize_data_structs();
0188   this->SetPrecision(1e-4, 1e-7);  // Default values
0189 }
0190 
0191 template <class QSS>
0192 inline G4QSStepper<QSS>::~G4QSStepper()
0193 {
0194   for (auto & i : simulator->SD) { free(i); }
0195 
0196   free(SUBSTEPS(this->simulator));
0197   free(this->simulator);
0198 }
0199 
0200 template <class QSS>
0201 inline void G4QSStepper<QSS>::Stepper(const G4double yInput[],
0202                                       const G4double dydx[],
0203                                             G4double hstep,
0204                                             G4double yOutput[],
0205                                             G4double yError[],
0206                                             G4double /*dydxOutput*/[])
0207 {
0208   Stepper(yInput, dydx, hstep, yOutput, yError);
0209 }
0210 
0211 template <class QSS>
0212 inline void G4QSStepper<QSS>::update_time()
0213 {
0214   auto* const sim = this->simulator;
0215 
0216   sim->time = sim->nextStateTime[0];
0217   sim->minIndex = 0;
0218 
0219   if (sim->nextStateTime[1] < sim->time) {
0220     sim->time = sim->nextStateTime[1];
0221     sim->minIndex = 1;
0222   }
0223   if (sim->nextStateTime[2] < sim->time) {
0224     sim->time = sim->nextStateTime[2];
0225     sim->minIndex = 2;
0226   }
0227   if (sim->nextStateTime[3] < sim->time) {
0228     sim->time = sim->nextStateTime[3];
0229     sim->minIndex = 3;
0230   }
0231   if (sim->nextStateTime[4] < sim->time) {
0232     sim->time = sim->nextStateTime[4];
0233     sim->minIndex = 4;
0234   }
0235   if (sim->nextStateTime[5] < sim->time) {
0236     sim->time = sim->nextStateTime[5];
0237     sim->minIndex = 5;
0238   }
0239 }
0240 
0241 template <class QSS>
0242 inline void G4QSStepper<QSS>::Stepper(const G4double yInput[],
0243                                       const G4double /*DyDx*/[],
0244                                             G4double max_length,
0245                                             G4double yOut[],
0246                                             G4double[] /*yErr[]*/)
0247 {
0248   G4double elapsed;
0249   G4double t, prev_time = 0;
0250   G4double length = 0.;
0251   G4int index;
0252 
0253   const G4int coeffs = method->order() + 1;
0254   G4double* tq = simulator->tq;
0255   G4double* tx = simulator->tx;
0256   G4double* dQRel = simulator->dQRel;
0257   G4double* dQMin = simulator->dQMin;
0258   G4double* lqu = simulator->lqu;
0259   G4double* x = simulator->x;
0260   G4int** SD = simulator->SD;
0261   G4int cf0, infCf0;
0262 
0263   CUR_SUBSTEP(simulator) = 0;
0264 
0265   this->save_substep(0, length);
0266 
0267   this->update_time();
0268   t = simulator->time;
0269   index = simulator->minIndex;
0270 
0271   while (length < max_length && t < Qss_misc::INF && CUR_SUBSTEP(simulator) < QSS_MAX_SUBSTEPS) {
0272     cf0 = index * coeffs;
0273     elapsed = t - tx[index];
0274     method->advance_time_x(cf0, elapsed);
0275     tx[index] = t;
0276     lqu[index] = dQRel[index] * std::fabs(x[cf0]);
0277     if (lqu[index] < dQMin[index]) {
0278       lqu[index] = dQMin[index];
0279     }
0280     method->update_quantized_state(index);
0281     tq[index] = t;
0282     method->next_time(index, t);
0283     for (G4int i = 0; i < 3; i++) {
0284       G4int j = SD[index][i];
0285       elapsed = t - tx[j];
0286       infCf0 = j * coeffs;
0287       if (elapsed > 0) {
0288         x[infCf0] = method->evaluate_x_poly(infCf0, elapsed, x);
0289         tx[j] = t;
0290       }
0291     }
0292 
0293     this->update_field();
0294     this->recompute_derivatives(index);
0295     method->recompute_next_times(SD[index], t);
0296 
0297     if (t > prev_time) {
0298       length += fVelocity * (t - prev_time);
0299       if (length <= max_length) { this->save_substep(t, length); }
0300       else { break; }
0301     }
0302 
0303     this->update_time();
0304     prev_time = t;
0305     t = simulator->time;
0306     index = simulator->minIndex;
0307   }
0308 
0309   if(CUR_SUBSTEP(simulator) >= QSS_MAX_SUBSTEPS) {
0310     max_length = length;
0311   }
0312 
0313   auto* const substep = &LAST_SUBSTEP_STRUCT(simulator);
0314   t = substep->start_time + (max_length - substep->len) / fVelocity;
0315 
0316   this->get_state_from_poly(substep->x, substep->tx, t, yOut);
0317 
0318   velocity_to_momentum(yOut);
0319 
0320   const G4int numberOfVariables = GetNumberOfVariables();
0321   for (G4int i = 0; i < numberOfVariables; ++i) {
0322     // Store Input and Final values, for possible use in calculating chord
0323     fyIn[i] = yInput[i];
0324     fyOut[i] = yOut[i];
0325   }
0326 
0327   fLastStepLength = max_length;
0328 }
0329 
0330 template<class QSS>
0331 inline G4double G4QSStepper<QSS>::DistChord() const
0332 {
0333   G4double yMid[6];
0334   const_cast<G4QSStepper<QSS>*>(this)->Interpolate(0.5, yMid);
0335 
0336   const G4ThreeVector begin = makeVector(fyIn, field_utils::Value3D::Position);
0337   const G4ThreeVector end = makeVector(fyOut, field_utils::Value3D::Position);
0338   const G4ThreeVector mid = makeVector(yMid, field_utils::Value3D::Position);
0339 
0340   return G4LineSection::Distline(mid, begin, end);
0341 }
0342 
0343 template <class QSS>
0344 inline void G4QSStepper<QSS>::Interpolate(G4double tau, G4double yOut[])
0345 {
0346   G4double length = tau * fLastStepLength;
0347   G4int idx = 0, j = LAST_SUBSTEP(simulator);
0348   G4double end_time;
0349 
0350   if (j >= 15) {
0351     G4int i = 0, k = j;
0352     idx = j >> 1;
0353     while (idx < k && i < j - 1) {
0354       if (length < SUBSTEP_LEN(simulator, idx)) {
0355         j = idx;
0356       } else if (length >= SUBSTEP_LEN(simulator, idx + 1)) {
0357         i = idx;
0358       } else {
0359         break;
0360       }
0361 
0362       idx = (i + j) >> 1;
0363     }
0364   }
0365   else {
0366     for (; idx < j && length >= SUBSTEP_LEN(simulator, idx + 1); idx++) {;}
0367   }
0368 
0369   auto* const substep = &SUBSTEP_STRUCT(simulator, idx);
0370   end_time = substep->start_time + (length - substep->len) / fVelocity;
0371 
0372   this->get_state_from_poly(substep->x, substep->tx, end_time, yOut);
0373 
0374   velocity_to_momentum(yOut);
0375 }
0376 
0377 template <class QSS>
0378 inline void G4QSStepper<QSS>::reset(const G4FieldTrack* track)
0379 {
0380   using Qss_misc::PXidx;
0381   using Qss_misc::PYidx;
0382   using Qss_misc::PZidx;
0383   using Qss_misc::VXidx;
0384   using Qss_misc::VYidx;
0385   using Qss_misc::VZidx;
0386 
0387   G4ThreeVector pos = track->GetPosition();
0388   G4ThreeVector momentum = track->GetMomentum();
0389 
0390   f_mass = track->GetRestMass();
0391   set_charge(track->GetCharge());
0392   set_gamma(momentum.mag2());
0393   G4double c_mg = cLight_local / fMassGamma;
0394   set_velocity(momentum.mag() * c_mg);
0395 
0396   method->reset_state(PXidx, pos.getX());
0397   method->reset_state(PYidx, pos.getY());
0398   method->reset_state(PZidx, pos.getZ());
0399 
0400   method->reset_state(VXidx, momentum.getX() * c_mg);
0401   method->reset_state(VYidx, momentum.getY() * c_mg);
0402   method->reset_state(VZidx, momentum.getZ() * c_mg);
0403 
0404   this->update_field();
0405   method->full_definition(get_coeff());
0406 
0407   method->recompute_all_state_times(0);
0408 
0409   simulator->time = 0;
0410 }
0411 
0412 template <class QSS>
0413 inline void G4QSStepper<QSS>::SetPrecision(G4double dq_rel, G4double dq_min)
0414 {
0415   G4double* dQMin = simulator->dQMin;
0416   G4double* dQRel = simulator->dQRel;
0417   G4int n_vars = simulator->states;
0418 
0419   if (dq_min <= 0) { dq_min = dq_rel * 1e-3; }
0420 
0421   for (G4int i = 0; i < n_vars; ++i) {
0422     dQRel[i] = dq_rel;
0423     dQMin[i] = dq_min;
0424   }
0425 }
0426 
0427 template <class QSS>
0428 inline void G4QSStepper<QSS>::initialize_data_structs()
0429 {
0430   auto sim = this->simulator;
0431   auto  states = (G4int*)calloc(Qss_misc::VAR_IDX_END, sizeof(G4int));
0432 
0433   sim->states = Qss_misc::VAR_IDX_END;
0434   sim->it = 0.;
0435 
0436   for (unsigned int i = 0; i < Qss_misc::VAR_IDX_END; i++) {
0437     sim->SD[i] = (G4int*)malloc(3 * sizeof(G4int));
0438   }
0439 
0440   sim->SD[0][states[0]++] = 3;
0441   sim->SD[0][states[0]++] = 4;
0442   sim->SD[0][states[0]++] = 5;
0443 
0444   sim->SD[1][states[1]++] = 3;
0445   sim->SD[1][states[1]++] = 4;
0446   sim->SD[1][states[1]++] = 5;
0447 
0448   sim->SD[2][states[2]++] = 3;
0449   sim->SD[2][states[2]++] = 4;
0450   sim->SD[2][states[2]++] = 5;
0451 
0452   sim->SD[3][states[3]++] = 0;
0453   sim->SD[3][states[3]++] = 4;
0454   sim->SD[3][states[3]++] = 5;
0455 
0456   sim->SD[4][states[4]++] = 1;
0457   sim->SD[4][states[4]++] = 3;
0458   sim->SD[4][states[4]++] = 5;
0459 
0460   sim->SD[5][states[5]++] = 2;
0461   sim->SD[5][states[5]++] = 3;
0462   sim->SD[5][states[5]++] = 4;
0463 
0464   free(states);
0465 }
0466 
0467 template <class QSS>
0468 inline QSS_simulator G4QSStepper<QSS>::build_simulator()
0469 {
0470   QSS_simulator sim = (QSS_simulator)malloc(sizeof(*sim));
0471   MAX_SUBSTEP(sim) = Qss_misc::MIN_SUBSTEPS;
0472   SUBSTEPS(sim) = (QSSSubstep)malloc(Qss_misc::MIN_SUBSTEPS * sizeof(*SUBSTEPS(sim)));
0473   return sim;
0474 }
0475 
0476 template <class QSS>
0477 inline void G4QSStepper<QSS>::recompute_derivatives(G4int index)
0478 {
0479   const G4int coeffs = method->order() + 1;
0480   G4double e;
0481   G4int idx = 0;
0482 
0483   e = simulator->time - simulator->tq[0];
0484   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0485   simulator->tq[0] = simulator->time;
0486 
0487   idx += coeffs;
0488   e = simulator->time - simulator->tq[1];
0489   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0490   simulator->tq[1] = simulator->time;
0491 
0492   idx += coeffs;
0493   e = simulator->time - simulator->tq[2];
0494   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0495   simulator->tq[2] = simulator->time;
0496 
0497   idx += coeffs;
0498   e = simulator->time - simulator->tq[3];
0499   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0500   simulator->tq[3] = simulator->time;
0501 
0502   idx += coeffs;
0503   e = simulator->time - simulator->tq[4];
0504   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0505   simulator->tq[4] = simulator->time;
0506 
0507   idx += coeffs;
0508   e = simulator->time - simulator->tq[5];
0509   if (likely(e > 0)) { method->advance_time_q(idx, e); }
0510   simulator->tq[5] = simulator->time;
0511 
0512   method->dependencies(index, get_coeff());
0513 }
0514 
0515 template <class QSS>
0516 inline void G4QSStepper<QSS>::update_field()
0517 {
0518   using Qss_misc::PXidx;
0519   using Qss_misc::PYidx;
0520   using Qss_misc::PZidx;
0521 
0522   const G4int order1 = method->order() + 1;
0523   G4double* const _field = simulator->alg;
0524   G4double* const _point = _field + order1;
0525 
0526   _point[PXidx] = simulator->x[PXidx];
0527   _point[PYidx] = simulator->x[PYidx * order1];
0528   _point[PZidx] = simulator->x[PZidx * order1];
0529 
0530   this->GetEquationOfMotion()->GetFieldValue(_point, _field);
0531 }
0532 
0533 template <class QSS>
0534 inline void G4QSStepper<QSS>::save_substep(G4double time, G4double length)
0535 {
0536   memcpy(CUR_SUBSTEP_X(simulator), simulator->x,
0537     (Qss_misc::VAR_IDX_END * (Qss_misc::MAX_QSS_STEPPER_ORDER + 2)) * sizeof(G4double));
0538 
0539   CUR_SUBSTEP_START(simulator) = time;
0540   CUR_SUBSTEP_LEN(simulator) = length;
0541   CUR_SUBSTEP(simulator)++;
0542 
0543   if (unlikely(CUR_SUBSTEP(simulator) == MAX_SUBSTEP(simulator))) {
0544     this->realloc_substeps();
0545   }
0546 }
0547 
0548 template <class QSS>
0549 inline void G4QSStepper<QSS>::realloc_substeps()
0550 {
0551   const G4int prev_index = MAX_SUBSTEP(simulator), new_index = 2 * prev_index;
0552 
0553   MAX_SUBSTEP(simulator) = new_index;
0554   SUBSTEPS(simulator) =
0555     (QSSSubstep)realloc(SUBSTEPS(simulator), new_index * sizeof(*SUBSTEPS(simulator)));
0556 }
0557 
0558 template <class QSS>
0559 inline void G4QSStepper<QSS>::get_state_from_poly(
0560   G4double* x, G4double* tx, G4double time, G4double* state)
0561 {
0562   unsigned int coeff_index = 0, i;
0563   const unsigned int x_order = method->order(), x_order1 = x_order + 1;
0564 
0565   for (i = 0; i < Qss_misc::VAR_IDX_END; ++i) {
0566     assert(tx[i] <= time);
0567     state[i] = method->evaluate_x_poly(coeff_index, time - tx[i], x);
0568     coeff_index += x_order1;
0569   }
0570 }
0571 
0572 template <class QSS>
0573 inline void G4QSStepper<QSS>::velocity_to_momentum(G4double* state)
0574 {
0575   using Qss_misc::VXidx;
0576   using Qss_misc::VYidx;
0577   using Qss_misc::VZidx;
0578   G4double coeff = fMassGamma / cLight_local;
0579 
0580   state[VXidx] *= coeff;
0581   state[VYidx] *= coeff;
0582   state[VZidx] *= coeff;
0583 }
0584 
0585 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!