Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4QSS2 Class Reference

#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

Definition at line 48 of file G4QSS2.hh.

Constructor & Destructor Documentation

◆ G4QSS2()

G4QSS2::G4QSS2 ( QSS_simulator sim)
inline

Definition at line 52 of file G4QSS2.hh.

52: simulator(sim) {}

Member Function Documentation

◆ advance_time_q()

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

Definition at line 381 of file G4QSS2.hh.

382 {
383 G4double* const p = simulator->q;
384 p[i] = p[i] + dt * p[i + 1];
385 }
G4double
double G4double
Definition G4Types.hh:83
QSS_simulator_::q
double q[Qss_misc::VAR_IDX_END *(Qss_misc::MAX_QSS_STEPPER_ORDER+1)]
Definition G4qss_misc.hh:100

◆ advance_time_x()

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

Definition at line 387 of file G4QSS2.hh.

388 {
389 G4double* const p = simulator->x;
390 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
391 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
392 p[i1] = p[i1] + 2 * dt * p[i2];
393 }
G4int
int G4int
Definition G4Types.hh:85
QSS_simulator_::x
double x[Qss_misc::VAR_IDX_END *(Qss_misc::MAX_QSS_STEPPER_ORDER+1)]
Definition G4qss_misc.hh:97

◆ dependencies()

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

Definition at line 83 of file G4QSS2.hh.

84 {
85 G4double* const x = simulator->q;
86 G4double* const der = simulator->x;
87 G4double* const alg = simulator->alg;
88
89 switch (i)
90 {
91 case 0:
92 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
93 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
94
95 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
96 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
97
98 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
99 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
100 return;
101 case 1:
102 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
103 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
104
105 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
106 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
107
108 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
109 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
110 return;
111 case 2:
112 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
113 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
114
115 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
116 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
117
118 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
119 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
120 return;
121 case 3:
122 der[1] = x[9];
123 der[2] = (x[10]) / 2;
124
125 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
126 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
127
128 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
129 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
130 return;
131 case 4:
132 der[4] = x[12];
133 der[5] = (x[13]) / 2;
134
135 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
136 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
137
138 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
139 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
140 return;
141 case 5:
142 der[7] = x[15];
143 der[8] = (x[16]) / 2;
144
145 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
146 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
147
148 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
149 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
150 return;
151 }
152 }
QSS_simulator_::alg
double alg[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:111

◆ evaluate_x_poly()

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

Definition at line 376 of file G4QSS2.hh.

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

◆ full_definition()

void G4QSS2::full_definition ( G4double coeff)
inline

Definition at line 58 of file G4QSS2.hh.

59 {
60 G4double* const x = simulator->q;
61 G4double* const dx = simulator->x;
62 G4double* const alg = simulator->alg;
63
64 dx[1] = x[9];
65 dx[2] = 0;
66
67 dx[4] = x[12];
68 dx[5] = 0;
69
70 dx[7] = x[15];
71 dx[8] = 0;
72
73 dx[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
74 dx[11] = 0;
75
76 dx[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
77 dx[14] = 0;
78
79 dx[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
80 dx[17] = 0;
81 }

◆ getSimulator()

QSS_simulator G4QSS2::getSimulator ( ) const
inline

Definition at line 54 of file G4QSS2.hh.

54{ return this->simulator; }

◆ next_time()

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

Definition at line 332 of file G4QSS2.hh.

333 {
334 const G4int cf2 = var * 3 + 2;
335 G4double* const x = simulator->x;
336 G4double* const lqu = simulator->lqu;
337 G4double* const time = simulator->nextStateTime;
338
339 if (x[cf2] != 0.0) {
340 time[var] = t + std::sqrt(lqu[var] / std::fabs(x[cf2]));
341 } else {
342 time[var] = Qss_misc::INF;
343 }
344 }
Qss_misc::INF
constexpr G4double INF
Definition G4qss_misc.hh:49
QSS_simulator_::nextStateTime
double nextStateTime[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:103
QSS_simulator_::lqu
double lqu[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:109

◆ order()

G4int G4QSS2::order ( ) const
inline

Definition at line 56 of file G4QSS2.hh.

56{ return 2; }

◆ recompute_all_state_times()

void G4QSS2::recompute_all_state_times ( G4double t)
inline

Definition at line 308 of file G4QSS2.hh.

309 {
310 G4double mpr;
311 G4double* const x = simulator->x;
312 G4double* const lqu = simulator->lqu;
313 G4double* const time = simulator->nextStateTime;
314
315 for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 3)
316 {
317 const G4int icf1 = icf0 + 1;
318
319 if (x[icf1] == 0)
320 {
321 time[var] = Qss_misc::INF;
322 }
323 else
324 {
325 mpr = lqu[var] / x[icf1];
326 if (mpr < 0) { mpr *= -1; }
327 time[var] = t + mpr;
328 }
329 }
330 }

◆ recompute_next_times()

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

Definition at line 154 of file G4QSS2.hh.

155 {
156 G4int i;
157 G4double* x = simulator->x;
158 G4double* q = simulator->q;
159 G4double* lqu = simulator->lqu;
160 G4double* time = simulator->nextStateTime;
161
162 for (i = 0; i < 3; i++)
163 {
164 const G4int var = inf[i];
165 const G4int icf0 = 3 * var;
166 const G4int icf1 = icf0 + 1;
167 const G4int icf2 = icf1 + 1;
168
169 time[var] = t;
170
171 if (std::fabs(q[icf0] - x[icf0]) < lqu[var])
172 {
173 G4double mpr = -1, mpr2;
174 G4double cf0 = q[icf0] + lqu[var] - x[icf0];
175 G4double cf1 = q[icf1] - x[icf1];
176 G4double cf2 = -x[icf2];
177 G4double cf0Alt = q[icf0] - lqu[var] - x[icf0];
178
179 if (unlikely(cf2 == 0 || (1000 * std::fabs(cf2)) < std::fabs(cf1)))
180 {
181 if (cf1 == 0) {
182 mpr = Qss_misc::INF;
183 } else
184 {
185 mpr = -cf0 / cf1;
186 mpr2 = -cf0Alt / cf1;
187 if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
188 }
189
190 if (mpr < 0) { mpr = Qss_misc::INF; }
191 }
192 else
193 {
194 static G4ThreadLocal unsigned long long okCalls=0LL, badCalls= 0LL;
195 constexpr G4double dangerZone = 1.0e+30;
196 static G4ThreadLocal G4double bigCf1_pr = dangerZone,
197 bigCf2_pr = dangerZone;
198 static G4ThreadLocal G4double bigCf1 = 0.0, bigCf2 = 0.0;
199 if( std::abs(cf1) > dangerZone || std::fabs(cf2) > dangerZone )
200 {
201 badCalls++;
202 if( badCalls == 1
203 || ( badCalls < 1000 && badCalls % 20 == 0 )
204 || ( 1000 < badCalls && badCalls < 10000 && badCalls % 100 == 0 )
205 || ( 10000 < badCalls && badCalls < 100000 && badCalls % 1000 == 0 )
206 || ( 100000 < badCalls && badCalls % 10000 == 0 )
207 || ( std::fabs(cf1) > 1.5 * bigCf1_pr || std::fabs(cf2) > 1.5 * bigCf2_pr )
208 )
209 {
210 std::cout << " cf1 = " << std::setw(15) << cf1 << " cf2= " << std::setw(15) << cf2
211 << " badCall # " << badCalls << " of " << badCalls + okCalls
212 << " fraction = " << double(badCalls) / double(badCalls+okCalls);
213
214 if( std::fabs(cf1) > 1.5 * bigCf1_pr ) { bigCf1_pr = std::fabs(cf1); std::cout << " Bigger cf1 "; }
215 if( std::fabs(cf2) > 1.5 * bigCf2_pr ) { bigCf2_pr = std::fabs(cf2); std::cout << " Bigger cf2 "; }
216 std::cout << std::endl;
217 }
218 if( std::fabs(cf1) > 1.5 * bigCf1 ) { bigCf1 = std::fabs(cf1); }
219 if( std::fabs(cf2) > 1.5 * bigCf2 ) { bigCf2 = std::fabs(cf2); }
220 }
221 else
222 {
223 okCalls++;
224 }
225
226#ifdef REPORT_CRITICAL_PROBLEM
227 constexpr unsigned int exp_limit= 140;
228 constexpr G4double limit= 1.0e+140; // std::pow(10,exp_limit));
229 assert( std::fabs( std::pow(10, exp_limit) - limit ) < 1.0e-14*limit );
230 G4bool bad_cf2fac= G4Log(std::fabs(cf2))
231 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*limit;
232 if( std::fabs(cf1) > limit
233 || G4Log(std::fabs(cf2))
234 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*exp_limit )
235 {
236 G4ExceptionDescription ermsg;
237 ermsg << "QSS2: Coefficients exceed tolerable values -- beyond " << limit << G4endl;
238 if( std::fabs(cf1) > limit )
239 {
240 ermsg << " |cf1| = " << cf1 << " is > " << limit << " (limit)";
241 }
242 if( bad_cf2fac)
243 {
244 ermsg << " bad cf2-factor: cf2 = " << cf2
245 << " product is > " << 2*limit << " (limit)";
246 }
247 G4Exception("QSS2::recompute_next_times",
248 "Field/Qss2-", FatalException, ermsg );
249 }
250#endif
251 G4double cf1_2 = cf1 * cf1;
252 G4double cf2_4 = 4 * cf2;
253 G4double disc1 = cf1_2 - cf2_4 * cf0;
254 G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
255 G4double cf2_d2 = 2 * cf2;
256
257 if (unlikely(disc1 < 0 && disc2 < 0)) // no real roots
258 {
259 mpr = Qss_misc::INF;
260 }
261 else if (disc2 < 0)
262 {
263 G4double sd, r1;
264 sd = std::sqrt(disc1);
265 r1 = (-cf1 + sd) / cf2_d2;
266 if (r1 > 0) {
267 mpr = r1;
268 } else {
269 mpr = Qss_misc::INF;
270 }
271 r1 = (-cf1 - sd) / cf2_d2;
272 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
273 }
274 else if (disc1 < 0)
275 {
276 G4double sd, r1;
277 sd = std::sqrt(disc2);
278 r1 = (-cf1 + sd) / cf2_d2;
279 if (r1 > 0) {
280 mpr = r1;
281 } else {
282 mpr = Qss_misc::INF;
283 }
284 r1 = (-cf1 - sd) / cf2_d2;
285 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
286 }
287 else
288 {
289 G4double sd1, r1, sd2, r2;
290 sd1 = std::sqrt(disc1);
291 sd2 = std::sqrt(disc2);
292 r1 = (-cf1 + sd1) / cf2_d2;
293 r2 = (-cf1 + sd2) / cf2_d2;
294 if (r1 > 0) { mpr = r1; }
295 else { mpr = Qss_misc::INF; }
296 r1 = (-cf1 - sd1) / cf2_d2;
297 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
298 if (r2 > 0 && r2 < mpr) { mpr = r2; }
299 r2 = (-cf1 - sd2) / cf2_d2;
300 if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
301 }
302 }
303 time[var] += mpr;
304 }
305 }
306 }
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4Log
G4double G4Log(G4double x)
Definition G4Log.hh:227
G4bool
bool G4bool
Definition G4Types.hh:86
G4endl
#define G4endl
Definition G4ios.hh:67
unlikely
#define unlikely(x)
Definition G4qss_misc.hh:56
G4ThreadLocal
#define G4ThreadLocal
Definition tls.hh:77

◆ reset_state()

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

Definition at line 356 of file G4QSS2.hh.

357 {
358 G4double* const x = simulator->x;
359 G4double* const q = simulator->q;
360 G4double* const tq = simulator->tq;
361 G4double* const tx = simulator->tx;
362 const G4int idx = 3 * i;
363
364 x[idx] = value;
365
366 simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
367 if (simulator->lqu[i] < simulator->dQMin[i])
368 {
369 simulator->lqu[i] = simulator->dQMin[i];
370 }
371
372 q[idx] = value;
373 q[idx + 1] = tq[i] = tx[i] = 0;
374 }
QSS_simulator_::tx
double tx[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:98
QSS_simulator_::dQRel
double dQRel[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:108
QSS_simulator_::tq
double tq[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:101
QSS_simulator_::dQMin
double dQMin[Qss_misc::VAR_IDX_END]
Definition G4qss_misc.hh:107

◆ update_quantized_state()

void G4QSS2::update_quantized_state ( G4int i)
inline

Definition at line 346 of file G4QSS2.hh.

347 {
348 const G4int cf0 = i * 3, cf1 = cf0 + 1;
349 G4double* const q = simulator->q;
350 G4double* const x = simulator->x;
351
352 q[cf0] = x[cf0];
353 q[cf1] = x[cf1];
354 }
The documentation for this class was generated from the following file: