Back to home page

EIC code displayed by LXR
 
 

    

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

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