Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QSS2.hh
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4QSS2
27//
28// G4QSS2 simulator
29
30// Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
31// --------------------------------------------------------------------
32#ifndef G4QSS2_HH
33#define G4QSS2_HH
34
35#include "G4Types.hh"
36#include "G4qss_misc.hh"
37
38#include <cmath>
39#include <cassert>
40
41#define REPORT_CRITICAL_PROBLEM 1
42
43#ifdef REPORT_CRITICAL_PROBLEM
44#include <cassert>
45#include "G4Log.hh"
46#endif
47
48/**
49 * @brief G4QSS2 defines the QSS2 simulator engine used in QSS field stepper.
50 */
51
52class G4QSS2
53{
54 public:
55
56 inline G4QSS2(QSS_simulator sim) : simulator(sim) {}
57
58 inline QSS_simulator getSimulator() const { return this->simulator; }
59
60 inline G4int order() const { return 2; }
61
62 inline void full_definition(G4double coeff)
63 {
64 G4double* const x = simulator->q;
65 G4double* const dx = simulator->x;
66 G4double* const alg = simulator->alg;
67
68 dx[1] = x[9];
69 dx[2] = 0;
70
71 dx[4] = x[12];
72 dx[5] = 0;
73
74 dx[7] = x[15];
75 dx[8] = 0;
76
77 dx[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
78 dx[11] = 0;
79
80 dx[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
81 dx[14] = 0;
82
83 dx[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
84 dx[17] = 0;
85 }
86
87 inline void dependencies(G4int i, G4double coeff)
88 {
89 G4double* const x = simulator->q;
90 G4double* const der = simulator->x;
91 G4double* const alg = simulator->alg;
92
93 switch (i)
94 {
95 case 0:
96 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
97 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
98
99 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
100 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
101
102 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
103 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
104 return;
105 case 1:
106 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
107 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
108
109 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
110 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
111
112 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
113 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
114 return;
115 case 2:
116 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
117 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
118
119 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
120 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
121
122 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
123 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
124 return;
125 case 3:
126 der[1] = x[9];
127 der[2] = (x[10]) / 2;
128
129 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
130 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
131
132 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
133 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
134 return;
135 case 4:
136 der[4] = x[12];
137 der[5] = (x[13]) / 2;
138
139 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
140 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
141
142 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
143 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
144 return;
145 case 5:
146 der[7] = x[15];
147 der[8] = (x[16]) / 2;
148
149 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
150 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
151
152 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
153 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
154 return;
155 }
156 }
157
158 inline void recompute_next_times(G4int* inf, G4double t)
159 {
160 G4int i;
161 G4double* x = simulator->x;
162 G4double* q = simulator->q;
163 G4double* lqu = simulator->lqu;
164 G4double* time = simulator->nextStateTime;
165
166 for (i = 0; i < 3; i++)
167 {
168 const G4int var = inf[i];
169 const G4int icf0 = 3 * var;
170 const G4int icf1 = icf0 + 1;
171 const G4int icf2 = icf1 + 1;
172
173 time[var] = t;
174
175 if (std::fabs(q[icf0] - x[icf0]) < lqu[var])
176 {
177 G4double mpr = -1, mpr2;
178 G4double cf0 = q[icf0] + lqu[var] - x[icf0];
179 G4double cf1 = q[icf1] - x[icf1];
180 G4double cf2 = -x[icf2];
181 G4double cf0Alt = q[icf0] - lqu[var] - x[icf0];
182
183 if (unlikely(cf2 == 0 || (1000 * std::fabs(cf2)) < std::fabs(cf1)))
184 {
185 if (cf1 == 0) {
186 mpr = Qss_misc::INF;
187 } else
188 {
189 mpr = -cf0 / cf1;
190 mpr2 = -cf0Alt / cf1;
191 if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
192 }
193
194 if (mpr < 0) { mpr = Qss_misc::INF; }
195 }
196 else
197 {
198 static G4ThreadLocal unsigned long long okCalls=0LL, badCalls= 0LL;
199 constexpr G4double dangerZone = 1.0e+30;
200 static G4ThreadLocal G4double bigCf1_pr = dangerZone,
201 bigCf2_pr = dangerZone;
202 static G4ThreadLocal G4double bigCf1 = 0.0, bigCf2 = 0.0;
203 if( std::abs(cf1) > dangerZone || std::fabs(cf2) > dangerZone )
204 {
205 badCalls++;
206 if( badCalls == 1
207 || ( badCalls < 1000 && badCalls % 20 == 0 )
208 || ( 1000 < badCalls && badCalls < 10000 && badCalls % 100 == 0 )
209 || ( 10000 < badCalls && badCalls < 100000 && badCalls % 1000 == 0 )
210 || ( 100000 < badCalls && badCalls % 10000 == 0 )
211 || ( std::fabs(cf1) > 1.5 * bigCf1_pr || std::fabs(cf2) > 1.5 * bigCf2_pr )
212 )
213 {
214 std::cout << " cf1 = " << std::setw(15) << cf1 << " cf2= " << std::setw(15) << cf2
215 << " badCall # " << badCalls << " of " << badCalls + okCalls
216 << " fraction = " << double(badCalls) / double(badCalls+okCalls);
217
218 if( std::fabs(cf1) > 1.5 * bigCf1_pr ) { bigCf1_pr = std::fabs(cf1); std::cout << " Bigger cf1 "; }
219 if( std::fabs(cf2) > 1.5 * bigCf2_pr ) { bigCf2_pr = std::fabs(cf2); std::cout << " Bigger cf2 "; }
220 std::cout << std::endl;
221 }
222 if( std::fabs(cf1) > 1.5 * bigCf1 ) { bigCf1 = std::fabs(cf1); }
223 if( std::fabs(cf2) > 1.5 * bigCf2 ) { bigCf2 = std::fabs(cf2); }
224 }
225 else
226 {
227 okCalls++;
228 }
229
230#ifdef REPORT_CRITICAL_PROBLEM
231 constexpr unsigned int exp_limit= 140;
232 constexpr G4double limit= 1.0e+140; // std::pow(10,exp_limit));
233 assert( std::fabs( std::pow(10, exp_limit) - limit ) < 1.0e-14*limit );
234 G4bool bad_cf2fac= G4Log(std::fabs(cf2))
235 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*limit;
236 if( std::fabs(cf1) > limit
237 || G4Log(std::fabs(cf2))
238 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*exp_limit )
239 {
240 G4ExceptionDescription ermsg;
241 ermsg << "QSS2: Coefficients exceed tolerable values -- beyond " << limit << G4endl;
242 if( std::fabs(cf1) > limit )
243 {
244 ermsg << " |cf1| = " << cf1 << " is > " << limit << " (limit)";
245 }
246 if( bad_cf2fac)
247 {
248 ermsg << " bad cf2-factor: cf2 = " << cf2
249 << " product is > " << 2*limit << " (limit)";
250 }
251 G4Exception("QSS2::recompute_next_times",
252 "Field/Qss2-", FatalException, ermsg );
253 }
254#endif
255 G4double cf1_2 = cf1 * cf1;
256 G4double cf2_4 = 4 * cf2;
257 G4double disc1 = cf1_2 - cf2_4 * cf0;
258 G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
259 G4double cf2_d2 = 2 * cf2;
260
261 if (unlikely(disc1 < 0 && disc2 < 0)) // no real roots
262 {
263 mpr = Qss_misc::INF;
264 }
265 else if (disc2 < 0)
266 {
267 G4double sd, r1;
268 sd = std::sqrt(disc1);
269 r1 = (-cf1 + sd) / cf2_d2;
270 if (r1 > 0) {
271 mpr = r1;
272 } else {
273 mpr = Qss_misc::INF;
274 }
275 r1 = (-cf1 - sd) / cf2_d2;
276 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
277 }
278 else if (disc1 < 0)
279 {
280 G4double sd, r1;
281 sd = std::sqrt(disc2);
282 r1 = (-cf1 + sd) / cf2_d2;
283 if (r1 > 0) {
284 mpr = r1;
285 } else {
286 mpr = Qss_misc::INF;
287 }
288 r1 = (-cf1 - sd) / cf2_d2;
289 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
290 }
291 else
292 {
293 G4double sd1, r1, sd2, r2;
294 sd1 = std::sqrt(disc1);
295 sd2 = std::sqrt(disc2);
296 r1 = (-cf1 + sd1) / cf2_d2;
297 r2 = (-cf1 + sd2) / cf2_d2;
298 if (r1 > 0) { mpr = r1; }
299 else { mpr = Qss_misc::INF; }
300 r1 = (-cf1 - sd1) / cf2_d2;
301 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
302 if (r2 > 0 && r2 < mpr) { mpr = r2; }
303 r2 = (-cf1 - sd2) / cf2_d2;
304 if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
305 }
306 }
307 time[var] += mpr;
308 }
309 }
310 }
311
312 inline void recompute_all_state_times(G4double t)
313 {
314 G4double mpr;
315 G4double* const x = simulator->x;
316 G4double* const lqu = simulator->lqu;
317 G4double* const time = simulator->nextStateTime;
318
319 for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 3)
320 {
321 const G4int icf1 = icf0 + 1;
322
323 if (x[icf1] == 0)
324 {
325 time[var] = Qss_misc::INF;
326 }
327 else
328 {
329 mpr = lqu[var] / x[icf1];
330 if (mpr < 0) { mpr *= -1; }
331 time[var] = t + mpr;
332 }
333 }
334 }
335
336 inline void next_time(G4int var, G4double t)
337 {
338 const G4int cf2 = var * 3 + 2;
339 G4double* const x = simulator->x;
340 G4double* const lqu = simulator->lqu;
341 G4double* const time = simulator->nextStateTime;
342
343 if (x[cf2] != 0.0) {
344 time[var] = t + std::sqrt(lqu[var] / std::fabs(x[cf2]));
345 } else {
346 time[var] = Qss_misc::INF;
347 }
348 }
349
350 inline void update_quantized_state(G4int i)
351 {
352 const G4int cf0 = i * 3, cf1 = cf0 + 1;
353 G4double* const q = simulator->q;
354 G4double* const x = simulator->x;
355
356 q[cf0] = x[cf0];
357 q[cf1] = x[cf1];
358 }
359
360 inline void reset_state(G4int i, G4double value)
361 {
362 G4double* const x = simulator->x;
363 G4double* const q = simulator->q;
364 G4double* const tq = simulator->tq;
365 G4double* const tx = simulator->tx;
366 const G4int idx = 3 * i;
367
368 x[idx] = value;
369
370 simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
371 if (simulator->lqu[i] < simulator->dQMin[i])
372 {
373 simulator->lqu[i] = simulator->dQMin[i];
374 }
375
376 q[idx] = value;
377 q[idx + 1] = tq[i] = tx[i] = 0;
378 }
379
380 inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
381 {
382 return (p[i + 2] * dt + p[i + 1]) * dt + p[i];
383 }
384
385 inline void advance_time_q(G4int i, G4double dt) // __attribute__((hot))
386 {
387 G4double* const p = simulator->q;
388 p[i] = p[i] + dt * p[i + 1];
389 }
390
391 inline void advance_time_x(G4int i, G4double dt) // __attribute__((hot))
392 {
393 G4double* const p = simulator->x;
394 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
395 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
396 p[i1] = p[i1] + 2 * dt * p[i2];
397 }
398
399 private:
400
401 QSS_simulator simulator;
402};
403
404#endif
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:169
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4qss_misc.hh
QSS_simulator
struct QSS_simulator_ * QSS_simulator
Definition G4qss_misc.hh:37
unlikely
#define unlikely(x)
Definition G4qss_misc.hh:67
G4QSS2::dependencies
void dependencies(G4int i, G4double coeff)
Definition G4QSS2.hh:87
G4QSS2::full_definition
void full_definition(G4double coeff)
Definition G4QSS2.hh:62
G4QSS2::advance_time_q
void advance_time_q(G4int i, G4double dt)
Definition G4QSS2.hh:385
G4QSS2::evaluate_x_poly
G4double evaluate_x_poly(G4int i, G4double dt, G4double *p)
Definition G4QSS2.hh:380
G4QSS2::getSimulator
QSS_simulator getSimulator() const
Definition G4QSS2.hh:58
G4QSS2::G4QSS2
G4QSS2(QSS_simulator sim)
Definition G4QSS2.hh:56
G4QSS2::order
G4int order() const
Definition G4QSS2.hh:60
G4QSS2::recompute_all_state_times
void recompute_all_state_times(G4double t)
Definition G4QSS2.hh:312
G4QSS2::recompute_next_times
void recompute_next_times(G4int *inf, G4double t)
Definition G4QSS2.hh:158
G4QSS2::advance_time_x
void advance_time_x(G4int i, G4double dt)
Definition G4QSS2.hh:391
G4QSS2::reset_state
void reset_state(G4int i, G4double value)
Definition G4QSS2.hh:360
G4QSS2::update_quantized_state
void update_quantized_state(G4int i)
Definition G4QSS2.hh:350
G4QSS2::next_time
void next_time(G4int var, G4double t)
Definition G4QSS2.hh:336
Qss_misc::INF
constexpr G4double INF
Definition G4qss_misc.hh:60
G4ThreadLocal
#define G4ThreadLocal
Definition tls.hh:77