| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4QSS2
27 //
28 // G4QSS2 simulator
29
30 // Authors: Lucio Santi, Rodrigo Castro (Univ.
31 // -------------------------------------------
32 #ifndef _G4QSS2_H_
33 #define _G4QSS2_H_ 1
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 class G4QSS2
49 {
50 public:
51
52 G4QSS2(QSS_simulator sim) : simulator(sim)
53
54 inline QSS_simulator getSimulator() const
55
56 inline G4int order() const { return 2; }
57
58 inline void full_definition(G4double coeff
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
74 dx[11] = 0;
75
76 dx[13] = coeff * (alg[0] * x[15] - alg[2
77 dx[14] = 0;
78
79 dx[16] = coeff * (alg[1] * x[9] - alg[0]
80 dx[17] = 0;
81 }
82
83 inline void dependencies(G4int i, G4double
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] -
93 der[11] = ((alg[2] * x[13] - x[16] *
94
95 der[13] = coeff * (alg[0] * x[15] -
96 der[14] = ((alg[0] * x[16] - alg[2]
97
98 der[16] = coeff * (alg[1] * x[9] - a
99 der[17] = (-coeff * (alg[0] * x[13]
100 return;
101 case 1:
102 der[10] = coeff * (alg[2] * x[12] -
103 der[11] = ((alg[2] * x[13] - x[16] *
104
105 der[13] = coeff * (alg[0] * x[15] -
106 der[14] = ((alg[0] * x[16] - alg[2]
107
108 der[16] = coeff * (alg[1] * x[9] - a
109 der[17] = (-coeff * (alg[0] * x[13]
110 return;
111 case 2:
112 der[10] = coeff * (alg[2] * x[12] -
113 der[11] = ((alg[2] * x[13] - x[16] *
114
115 der[13] = coeff * (alg[0] * x[15] -
116 der[14] = ((alg[0] * x[16] - alg[2]
117
118 der[16] = coeff * (alg[1] * x[9] - a
119 der[17] = (-coeff * (alg[0] * x[13]
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] -
126 der[14] = ((alg[0] * x[16] - alg[2]
127
128 der[16] = coeff * (alg[1] * x[9] - a
129 der[17] = (-coeff * (alg[0] * x[13]
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] -
136 der[11] = ((alg[2] * x[13] - x[16] *
137
138 der[16] = coeff * (alg[1] * x[9] - a
139 der[17] = (-coeff * (alg[0] * x[13]
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] -
146 der[11] = ((alg[2] * x[13] - x[16] *
147
148 der[13] = coeff * (alg[0] * x[15] -
149 der[14] = ((alg[0] * x[16] - alg[2]
150 return;
151 }
152 }
153
154 inline void recompute_next_times(G4int* in
155 {
156 G4int i;
157 G4double* x = simulator->x;
158 G4double* q = simulator->q;
159 G4double* lqu = simulator->lqu;
160 G4double* time = simulator->nextStateTim
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
172 {
173 G4double mpr = -1, mpr2;
174 G4double cf0 = q[icf0] + lqu[var] -
175 G4double cf1 = q[icf1] - x[icf1];
176 G4double cf2 = -x[icf2];
177 G4double cf0Alt = q[icf0] - lqu[var]
178
179 if (unlikely(cf2 == 0 || (1000 * std
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
188 }
189
190 if (mpr < 0) { mpr = Qss_misc::INF
191 }
192 else
193 {
194 static G4ThreadLocal unsigned long
195 constexpr G4double dangerZone = 1.
196 static G4ThreadLocal G4double bigC
197 bigC
198 static G4ThreadLocal G4double bigC
199 if( std::abs(cf1) > dangerZone ||
200 {
201 badCalls++;
202 if( badCalls == 1
203 || ( badCalls < 1000 && badCa
204 || ( 1000 < badCalls && bad
205 || ( 10000 < badCalls && bad
206 || ( 100000 < badCalls &&
207 || ( std::fabs(cf1) > 1.5 * b
208 )
209 {
210 std::cout << " cf1 = " << std:
211 << " badCall # " <<
212 << " fraction = " <
213
214 if( std::fabs(cf1) > 1.5 * big
215 if( std::fabs(cf2) > 1.5 * big
216 std::cout << std::endl;
217 }
218 if( std::fabs(cf1) > 1.5 * bigCf
219 if( std::fabs(cf2) > 1.5 * bigCf
220 }
221 else
222 {
223 okCalls++;
224 }
225
226 #ifdef REPORT_CRITICAL_PROBLEM
227 constexpr unsigned int exp_limit=
228 constexpr G4double limit= 1.0e+140
229 assert( std::fabs( std::pow(10, ex
230 G4bool bad_cf2fac= G4Log(std::fabs
231 + G4Log(std::max(
232 if( std::fabs(cf1) > limit
233 || G4Log(std::fabs(cf2))
234 + G4Log(std::max( std::fabs(cf
235 {
236 G4ExceptionDescription ermsg;
237 ermsg << "QSS2: Coefficients exc
238 if( std::fabs(cf1) > limit )
239 {
240 ermsg << " |cf1| = " << cf1 <<
241 }
242 if( bad_cf2fac)
243 {
244 ermsg << " bad cf2-factor: cf
245 << " product is > " << 2
246 }
247 G4Exception("QSS2::recompute_nex
248 "Field/Qss2-", Fatal
249 }
250 #endif
251 G4double cf1_2 = cf1 * cf1;
252 G4double cf2_4 = 4 * cf2;
253 G4double disc1 = cf1_2 - cf2_4 * c
254 G4double disc2 = cf1_2 - cf2_4 * c
255 G4double cf2_d2 = 2 * cf2;
256
257 if (unlikely(disc1 < 0 && disc2 <
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)) { mp
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)) { mp
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)) { mp
298 if (r2 > 0 && r2 < mpr) { mpr =
299 r2 = (-cf1 - sd2) / cf2_d2;
300 if ((r2 > 0) && (r2 < mpr)) { mp
301 }
302 }
303 time[var] += mpr;
304 }
305 }
306 }
307
308 inline void recompute_all_state_times(G4do
309 {
310 G4double mpr;
311 G4double* const x = simulator->x;
312 G4double* const lqu = simulator->lqu;
313 G4double* const time = simulator->nextSt
314
315 for (G4int var = 0, icf0 = 0; var < 6; v
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 }
331
332 inline void next_time(G4int var, G4double
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->nextSt
338
339 if (x[cf2] != 0.0) {
340 time[var] = t + std::sqrt(lqu[var] / s
341 } else {
342 time[var] = Qss_misc::INF;
343 }
344 }
345
346 inline void update_quantized_state(G4int i
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 }
355
356 inline void reset_state(G4int i, G4double
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]
367 if (simulator->lqu[i] < simulator->dQMin
368 {
369 simulator->lqu[i] = simulator->dQMin[i
370 }
371
372 q[idx] = value;
373 q[idx + 1] = tq[i] = tx[i] = 0;
374 }
375
376 inline G4double evaluate_x_poly(G4int i, G
377 {
378 return (p[i + 2] * dt + p[i + 1]) * dt +
379 }
380
381 inline void advance_time_q(G4int i, G4doub
382 {
383 G4double* const p = simulator->q;
384 p[i] = p[i] + dt * p[i + 1];
385 }
386
387 inline void advance_time_x(G4int i, G4doub
388 {
389 G4double* const p = simulator->x;
390 const G4int i0 = i, i1 = i0 + 1, i2 = i1
391 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0
392 p[i1] = p[i1] + 2 * dt * p[i2];
393 }
394
395 private:
396
397 QSS_simulator simulator;
398 };
399
400 #endif
401