G4QSS2 Class Reference

G4QSS2 defines the QSS2 simulator engine used in QSS field stepper. More...

#include <G4QSS2.hh>

Public Member Functions
	G4QSS2 (QSS_simulator sim)
QSS_simulator	getSimulator () const
G4int	order () const
void	full_definition (G4double coeff)
void	dependencies (G4int i, G4double coeff)
void	recompute_next_times (G4int *inf, G4double t)
void	recompute_all_state_times (G4double t)
void	next_time (G4int var, G4double t)
void	update_quantized_state (G4int i)
void	reset_state (G4int i, G4double value)
G4double	evaluate_x_poly (G4int i, G4double dt, G4double *p)
void	advance_time_q (G4int i, G4double dt)
void	advance_time_x (G4int i, G4double dt)

Detailed Description

G4QSS2 defines the QSS2 simulator engine used in QSS field stepper.

Definition at line 52 of file G4QSS2.hh.

Constructor & Destructor Documentation

G4QSS2()

G4QSS2::G4QSS2 ( QSS_simulator sim )

inline

Definition at line 56 of file G4QSS2.hh.

: simulator(sim) {}

Member Function Documentation

advance_time_q()

void G4QSS2::advance_time_q ( G4int i, G4double dt )

inline

Definition at line 385 of file G4QSS2.hh.

    {
      G4double* const p = simulator->q;
      p[i] = p[i] + dt * p[i + 1];
    }

advance_time_x()

void G4QSS2::advance_time_x ( G4int i, G4double dt )

inline

Definition at line 391 of file G4QSS2.hh.

    {
      G4double* const p = simulator->x;
      const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
      p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
      p[i1] = p[i1] + 2 * dt * p[i2];
    }

dependencies()

void G4QSS2::dependencies ( G4int i, G4double coeff )

inline

Definition at line 87 of file G4QSS2.hh.

    {
      G4double* const x = simulator->q;
      G4double* const der = simulator->x;
      G4double* const alg = simulator->alg;
 
      switch (i)
      {
        case 0:
          der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
          der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
 
          der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
          der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
 
          der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
          der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
          return;
        case 1:
          der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
          der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
 
          der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
          der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
 
          der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
          der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
          return;
        case 2:
          der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
          der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
 
          der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
          der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
 
          der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
          der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
          return;
        case 3:
          der[1] = x[9];
          der[2] = (x[10]) / 2;
 
          der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
          der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
 
          der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
          der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
          return;
        case 4:
          der[4] = x[12];
          der[5] = (x[13]) / 2;
 
          der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
          der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
 
          der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
          der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
          return;
        case 5:
          der[7] = x[15];
          der[8] = (x[16]) / 2;
 
          der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
          der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
 
          der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
          der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
          return;
      }
    }

evaluate_x_poly()

G4double G4QSS2::evaluate_x_poly ( G4int i, G4double dt, G4double * p )

inline

Definition at line 380 of file G4QSS2.hh.

    {
      return (p[i + 2] * dt + p[i + 1]) * dt + p[i];
    }

full_definition()

void G4QSS2::full_definition ( G4double coeff )

inline

Definition at line 62 of file G4QSS2.hh.

    {
      G4double* const x = simulator->q;
      G4double* const dx = simulator->x;
      G4double* const alg = simulator->alg;
 
      dx[1] = x[9];
      dx[2] = 0;
 
      dx[4] = x[12];
      dx[5] = 0;
 
      dx[7] = x[15];
      dx[8] = 0;
 
      dx[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
      dx[11] = 0;
 
      dx[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
      dx[14] = 0;
 
      dx[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
      dx[17] = 0;
    }

getSimulator()

QSS_simulator G4QSS2::getSimulator ( ) const

inline

Definition at line 58 of file G4QSS2.hh.

{ return this->simulator; }

next_time()

void G4QSS2::next_time ( G4int var, G4double t )

inline

Definition at line 336 of file G4QSS2.hh.

    {
      const G4int cf2 = var * 3 + 2;
      G4double* const x = simulator->x;
      G4double* const lqu = simulator->lqu;
      G4double* const time = simulator->nextStateTime;
 
      if (x[cf2] != 0.0) {
        time[var] = t + std::sqrt(lqu[var] / std::fabs(x[cf2]));
      } else {
        time[var] = Qss_misc::INF;
      }
    }

order()

G4int G4QSS2::order ( ) const

inline

Definition at line 60 of file G4QSS2.hh.

{ return 2; }

recompute_all_state_times()

void G4QSS2::recompute_all_state_times ( G4double t )

inline

Definition at line 312 of file G4QSS2.hh.

    {
      G4double mpr;
      G4double* const x = simulator->x;
      G4double* const lqu = simulator->lqu;
      G4double* const time = simulator->nextStateTime;
 
      for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 3)
      {
        const G4int icf1 = icf0 + 1;
 
        if (x[icf1] == 0)
        {
          time[var] = Qss_misc::INF;
        }
        else
        {
          mpr = lqu[var] / x[icf1];
          if (mpr < 0) { mpr *= -1; }
          time[var] = t + mpr;
        }
      }
    }

recompute_next_times()

void G4QSS2::recompute_next_times ( G4int * inf, G4double t )

inline

Definition at line 158 of file G4QSS2.hh.

    {
      G4int i;
      G4double* x = simulator->x;
      G4double* q = simulator->q;
      G4double* lqu = simulator->lqu;
      G4double* time = simulator->nextStateTime;
 
      for (i = 0; i < 3; i++)
      {
        const G4int var = inf[i];
        const G4int icf0 = 3 * var;
        const G4int icf1 = icf0 + 1;
        const G4int icf2 = icf1 + 1;
 
        time[var] = t;
 
        if (std::fabs(q[icf0] - x[icf0]) < lqu[var])
        {
          G4double mpr = -1, mpr2;
          G4double cf0 = q[icf0] + lqu[var] - x[icf0];
          G4double cf1 = q[icf1] - x[icf1];
          G4double cf2 = -x[icf2];
          G4double cf0Alt = q[icf0] - lqu[var] - x[icf0];
 
          if (unlikely(cf2 == 0 || (1000 * std::fabs(cf2)) < std::fabs(cf1)))
          {
            if (cf1 == 0) {
              mpr = Qss_misc::INF;
            } else
            {
              mpr = -cf0 / cf1;
              mpr2 = -cf0Alt / cf1;
              if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
            }
 
            if (mpr < 0) { mpr = Qss_misc::INF; }
          }
          else
          {
            static G4ThreadLocal unsigned long long okCalls=0LL, badCalls= 0LL;
            constexpr G4double dangerZone = 1.0e+30;
            static G4ThreadLocal G4double bigCf1_pr = dangerZone,
                                          bigCf2_pr = dangerZone;
            static G4ThreadLocal G4double bigCf1 = 0.0, bigCf2 = 0.0;
            if( std::abs(cf1) > dangerZone || std::fabs(cf2) > dangerZone )
            {
              badCalls++;
              if( badCalls == 1 
                 || ( badCalls < 1000 && badCalls % 20 == 0 )
                 || (   1000 < badCalls && badCalls <   10000 && badCalls %  100 == 0 )
                 || (  10000 < badCalls && badCalls <  100000 && badCalls % 1000 == 0 )
                 || ( 100000 < badCalls &&                       badCalls % 10000 == 0 )
                 || ( std::fabs(cf1) > 1.5 * bigCf1_pr || std::fabs(cf2) > 1.5 * bigCf2_pr )
                )
              {
                std::cout << " cf1 = " << std::setw(15) << cf1 << " cf2= " << std::setw(15) << cf2
                          << "  badCall # " << badCalls << " of " << badCalls + okCalls
                          << "  fraction = " << double(badCalls) / double(badCalls+okCalls);
 
                if( std::fabs(cf1) > 1.5 * bigCf1_pr ) { bigCf1_pr = std::fabs(cf1); std::cout << " Bigger cf1 "; }
                if( std::fabs(cf2) > 1.5 * bigCf2_pr ) { bigCf2_pr = std::fabs(cf2); std::cout << " Bigger cf2 "; }
                std::cout << std::endl;
              }
              if( std::fabs(cf1) > 1.5 * bigCf1 ) { bigCf1 = std::fabs(cf1); }
              if( std::fabs(cf2) > 1.5 * bigCf2 ) { bigCf2 = std::fabs(cf2); }
            }
            else
            {
              okCalls++;
            }
 
#ifdef REPORT_CRITICAL_PROBLEM
            constexpr unsigned int exp_limit= 140;
            constexpr G4double limit= 1.0e+140; // std::pow(10,exp_limit));
            assert( std::fabs( std::pow(10, exp_limit) - limit ) < 1.0e-14*limit );
            G4bool bad_cf2fac= G4Log(std::fabs(cf2))
                             + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*limit;
            if( std::fabs(cf1) > limit
               || G4Log(std::fabs(cf2))
                + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*exp_limit )
            {
              G4ExceptionDescription ermsg;
              ermsg << "QSS2: Coefficients exceed tolerable values -- beyond " << limit << G4endl;
              if( std::fabs(cf1) > limit )
              {
                ermsg << " |cf1| = " << cf1 << " is > " << limit << " (limit)";
              }
              if( bad_cf2fac)
              {
                ermsg << " bad cf2-factor:  cf2 = " << cf2
                      << " product is > " << 2*limit << " (limit)";
              }
              G4Exception("QSS2::recompute_next_times",
                          "Field/Qss2-", FatalException, ermsg ); 
            }
#endif
            G4double cf1_2 = cf1 * cf1;
            G4double cf2_4 = 4 * cf2;
            G4double disc1 = cf1_2 - cf2_4 * cf0;
            G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
            G4double cf2_d2 = 2 * cf2;
 
            if (unlikely(disc1 < 0 && disc2 < 0))  // no real roots
            {
              mpr = Qss_misc::INF;
            }
            else if (disc2 < 0)
            {
              G4double sd, r1;
              sd = std::sqrt(disc1);
              r1 = (-cf1 + sd) / cf2_d2;
              if (r1 > 0) {
                mpr = r1;
              } else {
                mpr = Qss_misc::INF;
              }
              r1 = (-cf1 - sd) / cf2_d2;
              if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
            }
            else if (disc1 < 0)
            {
              G4double sd, r1;
              sd = std::sqrt(disc2);
              r1 = (-cf1 + sd) / cf2_d2;
              if (r1 > 0) {
                mpr = r1;
              } else {
                mpr = Qss_misc::INF;
              }
              r1 = (-cf1 - sd) / cf2_d2;
              if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
            }
            else
            {
              G4double sd1, r1, sd2, r2;
              sd1 = std::sqrt(disc1);
              sd2 = std::sqrt(disc2);
              r1 = (-cf1 + sd1) / cf2_d2;
              r2 = (-cf1 + sd2) / cf2_d2;
              if (r1 > 0) { mpr = r1; }
              else { mpr = Qss_misc::INF; }
              r1 = (-cf1 - sd1) / cf2_d2;
              if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
              if (r2 > 0 && r2 < mpr) { mpr = r2; }
              r2 = (-cf1 - sd2) / cf2_d2;
              if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
            }
          }
          time[var] += mpr;
        }
      }
    }

reset_state()

void G4QSS2::reset_state ( G4int i, G4double value )

inline

Definition at line 360 of file G4QSS2.hh.

    {
      G4double* const x = simulator->x;
      G4double* const q = simulator->q;
      G4double* const tq = simulator->tq;
      G4double* const tx = simulator->tx;
      const G4int idx = 3 * i;
 
      x[idx] = value;
 
      simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
      if (simulator->lqu[i] < simulator->dQMin[i])
      {
        simulator->lqu[i] = simulator->dQMin[i];
      }
 
      q[idx] = value;
      q[idx + 1] = tq[i] = tx[i] = 0;
    }

update_quantized_state()

void G4QSS2::update_quantized_state ( G4int i )

inline

Definition at line 350 of file G4QSS2.hh.

    {
      const G4int cf0 = i * 3, cf1 = cf0 + 1;
      G4double* const q = simulator->q;
      G4double* const x = simulator->x;
 
      q[cf0] = x[cf0];
      q[cf1] = x[cf1];
    }

---

The documentation for this class was generated from the following file:

• G4QSS2.hh