include/Geant4/G4QSS2.hh

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 // G4QSS2
0027 //
0028 // G4QSS2 simulator
0029
0030 // Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
0031 // --------------------------------------------------------------------
0032 #ifndef G4QSS2_HH
0033 #define G4QSS2_HH
0034
0035 #include "G4Types.hh"
0036 #include "G4qss_misc.hh"
0037
0038 #include <cmath>
0039 #include <cassert>
0040
0041 #define  REPORT_CRITICAL_PROBLEM  1
0042
0043 #ifdef   REPORT_CRITICAL_PROBLEM
0044 #include <cassert>
0045 #include "G4Log.hh"
0046 #endif
0047
0048 /**
0049  * @brief G4QSS2 defines the QSS2 simulator engine used in QSS field stepper.
0050  */
0051
0052 class G4QSS2
0053 {
0054   public:
0055
0056     inline G4QSS2(QSS_simulator sim) : simulator(sim) {}
0057
0058     inline QSS_simulator getSimulator() const { return this->simulator; }
0059
0060     inline G4int order() const { return 2; }
0061
0062     inline void full_definition(G4double coeff)
0063     {
0064       G4double* const x = simulator->q;
0065       G4double* const dx = simulator->x;
0066       G4double* const alg = simulator->alg;
0067
0068       dx[1] = x[9];
0069       dx[2] = 0;
0070
0071       dx[4] = x[12];
0072       dx[5] = 0;
0073
0074       dx[7] = x[15];
0075       dx[8] = 0;
0076
0077       dx[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0078       dx[11] = 0;
0079
0080       dx[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0081       dx[14] = 0;
0082
0083       dx[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0084       dx[17] = 0;
0085     }
0086
0087     inline void dependencies(G4int i, G4double coeff)
0088     {
0089       G4double* const x = simulator->q;
0090       G4double* const der = simulator->x;
0091       G4double* const alg = simulator->alg;
0092
0093       switch (i)
0094       {
0095         case 0:
0096           der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0097           der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
0098
0099           der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0100           der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
0101
0102           der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0103           der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
0104           return;
0105         case 1:
0106           der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0107           der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
0108
0109           der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0110           der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
0111
0112           der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0113           der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
0114           return;
0115         case 2:
0116           der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0117           der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
0118
0119           der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0120           der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
0121
0122           der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0123           der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
0124           return;
0125         case 3:
0126           der[1] = x[9];
0127           der[2] = (x[10]) / 2;
0128
0129           der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0130           der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
0131
0132           der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0133           der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
0134           return;
0135         case 4:
0136           der[4] = x[12];
0137           der[5] = (x[13]) / 2;
0138
0139           der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0140           der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
0141
0142           der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
0143           der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
0144           return;
0145         case 5:
0146           der[7] = x[15];
0147           der[8] = (x[16]) / 2;
0148
0149           der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
0150           der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
0151
0152           der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
0153           der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
0154           return;
0155       }
0156     }
0157
0158     inline void recompute_next_times(G4int* inf, G4double t)
0159     {
0160       G4int i;
0161       G4double* x = simulator->x;
0162       G4double* q = simulator->q;
0163       G4double* lqu = simulator->lqu;
0164       G4double* time = simulator->nextStateTime;
0165
0166       for (i = 0; i < 3; i++)
0167       {
0168         const G4int var = inf[i];
0169         const G4int icf0 = 3 * var;
0170         const G4int icf1 = icf0 + 1;
0171         const G4int icf2 = icf1 + 1;
0172
0173         time[var] = t;
0174
0175         if (std::fabs(q[icf0] - x[icf0]) < lqu[var])
0176         {
0177           G4double mpr = -1, mpr2;
0178           G4double cf0 = q[icf0] + lqu[var] - x[icf0];
0179           G4double cf1 = q[icf1] - x[icf1];
0180           G4double cf2 = -x[icf2];
0181           G4double cf0Alt = q[icf0] - lqu[var] - x[icf0];
0182
0183           if (unlikely(cf2 == 0 || (1000 * std::fabs(cf2)) < std::fabs(cf1)))
0184           {
0185             if (cf1 == 0) {
0186               mpr = Qss_misc::INF;
0187             } else
0188             {
0189               mpr = -cf0 / cf1;
0190               mpr2 = -cf0Alt / cf1;
0191               if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
0192             }
0193
0194             if (mpr < 0) { mpr = Qss_misc::INF; }
0195           }
0196           else
0197           {
0198             static G4ThreadLocal unsigned long long okCalls=0LL, badCalls= 0LL;
0199             constexpr G4double dangerZone = 1.0e+30;
0200             static G4ThreadLocal G4double bigCf1_pr = dangerZone,
0201                                           bigCf2_pr = dangerZone;
0202             static G4ThreadLocal G4double bigCf1 = 0.0, bigCf2 = 0.0;
0203             if( std::abs(cf1) > dangerZone || std::fabs(cf2) > dangerZone )
0204             {
0205               badCalls++;
0206               if( badCalls == 1
0207                  || ( badCalls < 1000 && badCalls % 20 == 0 )
0208                  || (   1000 < badCalls && badCalls <   10000 && badCalls %  100 == 0 )
0209                  || (  10000 < badCalls && badCalls <  100000 && badCalls % 1000 == 0 )
0210                  || ( 100000 < badCalls &&                       badCalls % 10000 == 0 )
0211                  || ( std::fabs(cf1) > 1.5 * bigCf1_pr || std::fabs(cf2) > 1.5 * bigCf2_pr )
0212                 )
0213               {
0214                 std::cout << " cf1 = " << std::setw(15) << cf1 << " cf2= " << std::setw(15) << cf2
0215                           << "  badCall # " << badCalls << " of " << badCalls + okCalls
0216                           << "  fraction = " << double(badCalls) / double(badCalls+okCalls);
0217
0218                 if( std::fabs(cf1) > 1.5 * bigCf1_pr ) { bigCf1_pr = std::fabs(cf1); std::cout << " Bigger cf1 "; }
0219                 if( std::fabs(cf2) > 1.5 * bigCf2_pr ) { bigCf2_pr = std::fabs(cf2); std::cout << " Bigger cf2 "; }
0220                 std::cout << std::endl;
0221               }
0222               if( std::fabs(cf1) > 1.5 * bigCf1 ) { bigCf1 = std::fabs(cf1); }
0223               if( std::fabs(cf2) > 1.5 * bigCf2 ) { bigCf2 = std::fabs(cf2); }
0224             }
0225             else
0226             {
0227               okCalls++;
0228             }
0229
0230 #ifdef REPORT_CRITICAL_PROBLEM
0231             constexpr unsigned int exp_limit= 140;
0232             constexpr G4double limit= 1.0e+140; // std::pow(10,exp_limit));
0233             assert( std::fabs( std::pow(10, exp_limit) - limit ) < 1.0e-14*limit );
0234             G4bool bad_cf2fac= G4Log(std::fabs(cf2))
0235                              + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*limit;
0236             if( std::fabs(cf1) > limit
0237                || G4Log(std::fabs(cf2))
0238                 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*exp_limit )
0239             {
0240               G4ExceptionDescription ermsg;
0241               ermsg << "QSS2: Coefficients exceed tolerable values -- beyond " << limit << G4endl;
0242               if( std::fabs(cf1) > limit )
0243               {
0244                 ermsg << " |cf1| = " << cf1 << " is > " << limit << " (limit)";
0245               }
0246               if( bad_cf2fac)
0247               {
0248                 ermsg << " bad cf2-factor:  cf2 = " << cf2
0249                       << " product is > " << 2*limit << " (limit)";
0250               }
0251               G4Exception("QSS2::recompute_next_times",
0252                           "Field/Qss2-", FatalException, ermsg );
0253             }
0254 #endif
0255             G4double cf1_2 = cf1 * cf1;
0256             G4double cf2_4 = 4 * cf2;
0257             G4double disc1 = cf1_2 - cf2_4 * cf0;
0258             G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
0259             G4double cf2_d2 = 2 * cf2;
0260
0261             if (unlikely(disc1 < 0 && disc2 < 0))  // no real roots
0262             {
0263               mpr = Qss_misc::INF;
0264             }
0265             else if (disc2 < 0)
0266             {
0267               G4double sd, r1;
0268               sd = std::sqrt(disc1);
0269               r1 = (-cf1 + sd) / cf2_d2;
0270               if (r1 > 0) {
0271                 mpr = r1;
0272               } else {
0273                 mpr = Qss_misc::INF;
0274               }
0275               r1 = (-cf1 - sd) / cf2_d2;
0276               if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0277             }
0278             else if (disc1 < 0)
0279             {
0280               G4double sd, r1;
0281               sd = std::sqrt(disc2);
0282               r1 = (-cf1 + sd) / cf2_d2;
0283               if (r1 > 0) {
0284                 mpr = r1;
0285               } else {
0286                 mpr = Qss_misc::INF;
0287               }
0288               r1 = (-cf1 - sd) / cf2_d2;
0289               if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0290             }
0291             else
0292             {
0293               G4double sd1, r1, sd2, r2;
0294               sd1 = std::sqrt(disc1);
0295               sd2 = std::sqrt(disc2);
0296               r1 = (-cf1 + sd1) / cf2_d2;
0297               r2 = (-cf1 + sd2) / cf2_d2;
0298               if (r1 > 0) { mpr = r1; }
0299               else { mpr = Qss_misc::INF; }
0300               r1 = (-cf1 - sd1) / cf2_d2;
0301               if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0302               if (r2 > 0 && r2 < mpr) { mpr = r2; }
0303               r2 = (-cf1 - sd2) / cf2_d2;
0304               if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
0305             }
0306           }
0307           time[var] += mpr;
0308         }
0309       }
0310     }
0311
0312     inline void recompute_all_state_times(G4double t)
0313     {
0314       G4double mpr;
0315       G4double* const x = simulator->x;
0316       G4double* const lqu = simulator->lqu;
0317       G4double* const time = simulator->nextStateTime;
0318
0319       for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 3)
0320       {
0321         const G4int icf1 = icf0 + 1;
0322
0323         if (x[icf1] == 0)
0324         {
0325           time[var] = Qss_misc::INF;
0326         }
0327         else
0328         {
0329           mpr = lqu[var] / x[icf1];
0330           if (mpr < 0) { mpr *= -1; }
0331           time[var] = t + mpr;
0332         }
0333       }
0334     }
0335
0336     inline void next_time(G4int var, G4double t)
0337     {
0338       const G4int cf2 = var * 3 + 2;
0339       G4double* const x = simulator->x;
0340       G4double* const lqu = simulator->lqu;
0341       G4double* const time = simulator->nextStateTime;
0342
0343       if (x[cf2] != 0.0) {
0344         time[var] = t + std::sqrt(lqu[var] / std::fabs(x[cf2]));
0345       } else {
0346         time[var] = Qss_misc::INF;
0347       }
0348     }
0349
0350     inline void update_quantized_state(G4int i)
0351     {
0352       const G4int cf0 = i * 3, cf1 = cf0 + 1;
0353       G4double* const q = simulator->q;
0354       G4double* const x = simulator->x;
0355
0356       q[cf0] = x[cf0];
0357       q[cf1] = x[cf1];
0358     }
0359
0360     inline void reset_state(G4int i, G4double value)
0361     {
0362       G4double* const x = simulator->x;
0363       G4double* const q = simulator->q;
0364       G4double* const tq = simulator->tq;
0365       G4double* const tx = simulator->tx;
0366       const G4int idx = 3 * i;
0367
0368       x[idx] = value;
0369
0370       simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
0371       if (simulator->lqu[i] < simulator->dQMin[i])
0372       {
0373         simulator->lqu[i] = simulator->dQMin[i];
0374       }
0375
0376       q[idx] = value;
0377       q[idx + 1] = tq[i] = tx[i] = 0;
0378     }
0379
0380     inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
0381     {
0382       return (p[i + 2] * dt + p[i + 1]) * dt + p[i];
0383     }
0384
0385     inline void advance_time_q(G4int i, G4double dt)  // __attribute__((hot))
0386     {
0387       G4double* const p = simulator->q;
0388       p[i] = p[i] + dt * p[i + 1];
0389     }
0390
0391     inline void advance_time_x(G4int i, G4double dt)  // __attribute__((hot))
0392     {
0393       G4double* const p = simulator->x;
0394       const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
0395       p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
0396       p[i1] = p[i1] + 2 * dt * p[i2];
0397     }
0398
0399   private:
0400
0401     QSS_simulator simulator;
0402 };
0403
0404 #endif