| Geant4 Cross Reference |
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 // G4QSStepper
27 //
28 // QSS Integrator Stepper
29
30 // Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires) - 2018-2021
31 // --------------------------------------------------------------------
32 #ifndef QSS_Stepper_HH
33 #define QSS_Stepper_HH 1
34
35 #include "G4FieldTrack.hh"
36 #include "G4FieldUtils.hh"
37 #include "G4LineSection.hh"
38 #include "G4MagIntegratorStepper.hh"
39 #include "G4QSS2.hh"
40 #include "G4QSS3.hh"
41 #include "G4QSSDriver.hh"
42 #include "G4QSSMessenger.hh"
43 #include "G4VIntegrationDriver.hh"
44 #include "G4qss_misc.hh"
45
46 #include <cmath>
47 #include <cassert>
48
49 // Maximum allowed number of QSS substeps per integration step
50 #define QSS_MAX_SUBSTEPS 1000
51
52 template <class QSS>
53 class G4QSStepper : public G4MagIntegratorStepper
54 {
55 public:
56
57 G4QSStepper(G4EquationOfMotion* EqRhs,
58 G4int numberOfVariables = 6,
59 G4bool primary = true);
60 ~G4QSStepper() override;
61
62 void Stepper(const G4double y[],
63 const G4double dydx[],
64 G4double h,
65 G4double yout[],
66 G4double yerr[]) override;
67
68 void Stepper(const G4double yInput[],
69 const G4double dydx[],
70 G4double hstep,
71 G4double yOutput[],
72 G4double yError[],
73 G4double dydxOutput[]);
74
75 // For calculating the output at the tau fraction of Step
76 //
77 inline void SetupInterpolation() {}
78 inline void Interpolate(G4double tau, G4double yOut[]);
79
80 G4double DistChord() const override;
81
82 G4int IntegratorOrder() const override { return method->order(); }
83
84 void reset(const G4FieldTrack* track);
85
86 void SetPrecision(G4double dq_rel, G4double dq_min);
87 // precision parameters for QSS method
88
89 static G4QSStepper<G4QSS2>* build_QSS2(G4EquationOfMotion* EqRhs,
90 G4int numberOfVariables = 6,
91 G4bool primary = true);
92
93 static G4QSStepper<G4QSS3>* build_QSS3(G4EquationOfMotion* EqRhs,
94 G4int numberOfVariables = 6,
95 G4bool primary = true);
96
97 inline G4EquationOfMotion* GetSpecificEquation() { return GetEquationOfMotion(); }
98
99 inline const field_utils::State& GetYOut() const { return fyOut; }
100
101 inline G4double GetLastStepLength() { return fLastStepLength; }
102
103 private:
104
105 G4QSStepper(QSS* method,
106 G4EquationOfMotion* EqRhs,
107 G4int numberOfVariables = 6,
108 G4bool primary = true);
109
110 void initialize_data_structs();
111 static QSS_simulator build_simulator();
112
113 inline void update_field();
114 inline void save_substep(G4double time, G4double length);
115
116 inline void realloc_substeps();
117 inline void get_state_from_poly(G4double* x, G4double* tx,
118 G4double time, G4double* state);
119
120 inline void recompute_derivatives(int index);
121 inline void update_time();
122
123 inline G4double get_coeff() { return fCoeff_local; }
124
125 inline void set_coeff(G4double coeff) { fCoeff_local = coeff; }
126
127 inline void set_charge(G4double q)
128 {
129 f_charge_c2 = q * cLight_local * cLight_local; // 89875.5178737;
130 }
131
132 inline G4double get_qc2() { return f_charge_c2; }
133
134 inline void set_mg() { fMassGamma = f_mass * fGamma2; }
135
136 inline void set_gamma2(G4double gamma2) { fGamma2 = gamma2; }
137 inline void set_velocity(G4double v) { fVelocity = v; }
138
139 inline void velocity_to_momentum(G4double* state);
140
141 inline void set_gamma(G4double p_sq)
142 {
143 set_gamma2(std::sqrt(p_sq / (f_mass * f_mass) + 1));
144 set_mg();
145 set_coeff(get_qc2() / fMassGamma);
146 }
147
148 private:
149
150 QSS_simulator simulator;
151 QSS* method;
152
153 // State
154 //
155 G4double fLastStepLength;
156 field_utils::State fyIn, fyOut;
157
158 G4double f_mass;
159 static constexpr G4double cLight_local = 299.792458; // should use CLHEP
160 G4double f_charge_c2;
161 G4double fMassGamma;
162 G4double fGamma2;
163 G4double fCoeff_local;
164 G4double fVelocity;
165 };
166
167 using G4QSStepper_QSS2 = G4QSStepper<G4QSS2>;
168 using G4QSStepper_QSS3 = G4QSStepper<G4QSS3>;
169
170 template <class QSS>
171 inline G4QSStepper<QSS>::G4QSStepper(QSS* qss, G4EquationOfMotion* EqRhs,
172 G4int noIntegrationVariables, G4bool)
173 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
174 simulator(qss->getSimulator()),
175 method(qss)
176 {
177 SetIsQSS(true); // Replaces virtual method IsQSS
178 fLastStepLength = -1.0;
179
180 f_mass = 0;
181 f_charge_c2 = 0;
182 fMassGamma = 0;
183 fGamma2 = 0;
184 fCoeff_local = 0;
185 fVelocity = 0;
186
187 this->initialize_data_structs();
188 this->SetPrecision(1e-4, 1e-7); // Default values
189 }
190
191 template <class QSS>
192 inline G4QSStepper<QSS>::~G4QSStepper()
193 {
194 for (auto & i : simulator->SD) { free(i); }
195
196 free(SUBSTEPS(this->simulator));
197 free(this->simulator);
198 }
199
200 template <class QSS>
201 inline void G4QSStepper<QSS>::Stepper(const G4double yInput[],
202 const G4double dydx[],
203 G4double hstep,
204 G4double yOutput[],
205 G4double yError[],
206 G4double /*dydxOutput*/[])
207 {
208 Stepper(yInput, dydx, hstep, yOutput, yError);
209 }
210
211 template <class QSS>
212 inline void G4QSStepper<QSS>::update_time()
213 {
214 auto* const sim = this->simulator;
215
216 sim->time = sim->nextStateTime[0];
217 sim->minIndex = 0;
218
219 if (sim->nextStateTime[1] < sim->time) {
220 sim->time = sim->nextStateTime[1];
221 sim->minIndex = 1;
222 }
223 if (sim->nextStateTime[2] < sim->time) {
224 sim->time = sim->nextStateTime[2];
225 sim->minIndex = 2;
226 }
227 if (sim->nextStateTime[3] < sim->time) {
228 sim->time = sim->nextStateTime[3];
229 sim->minIndex = 3;
230 }
231 if (sim->nextStateTime[4] < sim->time) {
232 sim->time = sim->nextStateTime[4];
233 sim->minIndex = 4;
234 }
235 if (sim->nextStateTime[5] < sim->time) {
236 sim->time = sim->nextStateTime[5];
237 sim->minIndex = 5;
238 }
239 }
240
241 template <class QSS>
242 inline void G4QSStepper<QSS>::Stepper(const G4double yInput[],
243 const G4double /*DyDx*/[],
244 G4double max_length,
245 G4double yOut[],
246 G4double[] /*yErr[]*/)
247 {
248 G4double elapsed;
249 G4double t, prev_time = 0;
250 G4double length = 0.;
251 G4int index;
252
253 const G4int coeffs = method->order() + 1;
254 G4double* tq = simulator->tq;
255 G4double* tx = simulator->tx;
256 G4double* dQRel = simulator->dQRel;
257 G4double* dQMin = simulator->dQMin;
258 G4double* lqu = simulator->lqu;
259 G4double* x = simulator->x;
260 G4int** SD = simulator->SD;
261 G4int cf0, infCf0;
262
263 CUR_SUBSTEP(simulator) = 0;
264
265 this->save_substep(0, length);
266
267 this->update_time();
268 t = simulator->time;
269 index = simulator->minIndex;
270
271 while (length < max_length && t < Qss_misc::INF && CUR_SUBSTEP(simulator) < QSS_MAX_SUBSTEPS) {
272 cf0 = index * coeffs;
273 elapsed = t - tx[index];
274 method->advance_time_x(cf0, elapsed);
275 tx[index] = t;
276 lqu[index] = dQRel[index] * std::fabs(x[cf0]);
277 if (lqu[index] < dQMin[index]) {
278 lqu[index] = dQMin[index];
279 }
280 method->update_quantized_state(index);
281 tq[index] = t;
282 method->next_time(index, t);
283 for (G4int i = 0; i < 3; i++) {
284 G4int j = SD[index][i];
285 elapsed = t - tx[j];
286 infCf0 = j * coeffs;
287 if (elapsed > 0) {
288 x[infCf0] = method->evaluate_x_poly(infCf0, elapsed, x);
289 tx[j] = t;
290 }
291 }
292
293 this->update_field();
294 this->recompute_derivatives(index);
295 method->recompute_next_times(SD[index], t);
296
297 if (t > prev_time) {
298 length += fVelocity * (t - prev_time);
299 if (length <= max_length) { this->save_substep(t, length); }
300 else { break; }
301 }
302
303 this->update_time();
304 prev_time = t;
305 t = simulator->time;
306 index = simulator->minIndex;
307 }
308
309 if(CUR_SUBSTEP(simulator) >= QSS_MAX_SUBSTEPS) {
310 max_length = length;
311 }
312
313 auto* const substep = &LAST_SUBSTEP_STRUCT(simulator);
314 t = substep->start_time + (max_length - substep->len) / fVelocity;
315
316 this->get_state_from_poly(substep->x, substep->tx, t, yOut);
317
318 velocity_to_momentum(yOut);
319
320 const G4int numberOfVariables = GetNumberOfVariables();
321 for (G4int i = 0; i < numberOfVariables; ++i) {
322 // Store Input and Final values, for possible use in calculating chord
323 fyIn[i] = yInput[i];
324 fyOut[i] = yOut[i];
325 }
326
327 fLastStepLength = max_length;
328 }
329
330 template<class QSS>
331 inline G4double G4QSStepper<QSS>::DistChord() const
332 {
333 G4double yMid[6];
334 const_cast<G4QSStepper<QSS>*>(this)->Interpolate(0.5, yMid);
335
336 const G4ThreeVector begin = makeVector(fyIn, field_utils::Value3D::Position);
337 const G4ThreeVector end = makeVector(fyOut, field_utils::Value3D::Position);
338 const G4ThreeVector mid = makeVector(yMid, field_utils::Value3D::Position);
339
340 return G4LineSection::Distline(mid, begin, end);
341 }
342
343 template <class QSS>
344 inline void G4QSStepper<QSS>::Interpolate(G4double tau, G4double yOut[])
345 {
346 G4double length = tau * fLastStepLength;
347 G4int idx = 0, j = LAST_SUBSTEP(simulator);
348 G4double end_time;
349
350 if (j >= 15) {
351 G4int i = 0, k = j;
352 idx = j >> 1;
353 while (idx < k && i < j - 1) {
354 if (length < SUBSTEP_LEN(simulator, idx)) {
355 j = idx;
356 } else if (length >= SUBSTEP_LEN(simulator, idx + 1)) {
357 i = idx;
358 } else {
359 break;
360 }
361
362 idx = (i + j) >> 1;
363 }
364 }
365 else {
366 for (; idx < j && length >= SUBSTEP_LEN(simulator, idx + 1); idx++) {;}
367 }
368
369 auto* const substep = &SUBSTEP_STRUCT(simulator, idx);
370 end_time = substep->start_time + (length - substep->len) / fVelocity;
371
372 this->get_state_from_poly(substep->x, substep->tx, end_time, yOut);
373
374 velocity_to_momentum(yOut);
375 }
376
377 template <class QSS>
378 inline void G4QSStepper<QSS>::reset(const G4FieldTrack* track)
379 {
380 using Qss_misc::PXidx;
381 using Qss_misc::PYidx;
382 using Qss_misc::PZidx;
383 using Qss_misc::VXidx;
384 using Qss_misc::VYidx;
385 using Qss_misc::VZidx;
386
387 G4ThreeVector pos = track->GetPosition();
388 G4ThreeVector momentum = track->GetMomentum();
389
390 f_mass = track->GetRestMass();
391 set_charge(track->GetCharge());
392 set_gamma(momentum.mag2());
393 G4double c_mg = cLight_local / fMassGamma;
394 set_velocity(momentum.mag() * c_mg);
395
396 method->reset_state(PXidx, pos.getX());
397 method->reset_state(PYidx, pos.getY());
398 method->reset_state(PZidx, pos.getZ());
399
400 method->reset_state(VXidx, momentum.getX() * c_mg);
401 method->reset_state(VYidx, momentum.getY() * c_mg);
402 method->reset_state(VZidx, momentum.getZ() * c_mg);
403
404 this->update_field();
405 method->full_definition(get_coeff());
406
407 method->recompute_all_state_times(0);
408
409 simulator->time = 0;
410 }
411
412 template <class QSS>
413 inline void G4QSStepper<QSS>::SetPrecision(G4double dq_rel, G4double dq_min)
414 {
415 G4double* dQMin = simulator->dQMin;
416 G4double* dQRel = simulator->dQRel;
417 G4int n_vars = simulator->states;
418
419 if (dq_min <= 0) { dq_min = dq_rel * 1e-3; }
420
421 for (G4int i = 0; i < n_vars; ++i) {
422 dQRel[i] = dq_rel;
423 dQMin[i] = dq_min;
424 }
425 }
426
427 template <class QSS>
428 inline void G4QSStepper<QSS>::initialize_data_structs()
429 {
430 auto sim = this->simulator;
431 auto states = (G4int*)calloc(Qss_misc::VAR_IDX_END, sizeof(G4int));
432
433 sim->states = Qss_misc::VAR_IDX_END;
434 sim->it = 0.;
435
436 for (unsigned int i = 0; i < Qss_misc::VAR_IDX_END; i++) {
437 sim->SD[i] = (G4int*)malloc(3 * sizeof(G4int));
438 }
439
440 sim->SD[0][states[0]++] = 3;
441 sim->SD[0][states[0]++] = 4;
442 sim->SD[0][states[0]++] = 5;
443
444 sim->SD[1][states[1]++] = 3;
445 sim->SD[1][states[1]++] = 4;
446 sim->SD[1][states[1]++] = 5;
447
448 sim->SD[2][states[2]++] = 3;
449 sim->SD[2][states[2]++] = 4;
450 sim->SD[2][states[2]++] = 5;
451
452 sim->SD[3][states[3]++] = 0;
453 sim->SD[3][states[3]++] = 4;
454 sim->SD[3][states[3]++] = 5;
455
456 sim->SD[4][states[4]++] = 1;
457 sim->SD[4][states[4]++] = 3;
458 sim->SD[4][states[4]++] = 5;
459
460 sim->SD[5][states[5]++] = 2;
461 sim->SD[5][states[5]++] = 3;
462 sim->SD[5][states[5]++] = 4;
463
464 free(states);
465 }
466
467 template <class QSS>
468 inline QSS_simulator G4QSStepper<QSS>::build_simulator()
469 {
470 QSS_simulator sim = (QSS_simulator)malloc(sizeof(*sim));
471 MAX_SUBSTEP(sim) = Qss_misc::MIN_SUBSTEPS;
472 SUBSTEPS(sim) = (QSSSubstep)malloc(Qss_misc::MIN_SUBSTEPS * sizeof(*SUBSTEPS(sim)));
473 return sim;
474 }
475
476 template <class QSS>
477 inline void G4QSStepper<QSS>::recompute_derivatives(G4int index)
478 {
479 const G4int coeffs = method->order() + 1;
480 G4double e;
481 G4int idx = 0;
482
483 e = simulator->time - simulator->tq[0];
484 if (likely(e > 0)) { method->advance_time_q(idx, e); }
485 simulator->tq[0] = simulator->time;
486
487 idx += coeffs;
488 e = simulator->time - simulator->tq[1];
489 if (likely(e > 0)) { method->advance_time_q(idx, e); }
490 simulator->tq[1] = simulator->time;
491
492 idx += coeffs;
493 e = simulator->time - simulator->tq[2];
494 if (likely(e > 0)) { method->advance_time_q(idx, e); }
495 simulator->tq[2] = simulator->time;
496
497 idx += coeffs;
498 e = simulator->time - simulator->tq[3];
499 if (likely(e > 0)) { method->advance_time_q(idx, e); }
500 simulator->tq[3] = simulator->time;
501
502 idx += coeffs;
503 e = simulator->time - simulator->tq[4];
504 if (likely(e > 0)) { method->advance_time_q(idx, e); }
505 simulator->tq[4] = simulator->time;
506
507 idx += coeffs;
508 e = simulator->time - simulator->tq[5];
509 if (likely(e > 0)) { method->advance_time_q(idx, e); }
510 simulator->tq[5] = simulator->time;
511
512 method->dependencies(index, get_coeff());
513 }
514
515 template <class QSS>
516 inline void G4QSStepper<QSS>::update_field()
517 {
518 using Qss_misc::PXidx;
519 using Qss_misc::PYidx;
520 using Qss_misc::PZidx;
521
522 const G4int order1 = method->order() + 1;
523 G4double* const _field = simulator->alg;
524 G4double* const _point = _field + order1;
525
526 _point[PXidx] = simulator->x[PXidx];
527 _point[PYidx] = simulator->x[PYidx * order1];
528 _point[PZidx] = simulator->x[PZidx * order1];
529
530 this->GetEquationOfMotion()->GetFieldValue(_point, _field);
531 }
532
533 template <class QSS>
534 inline void G4QSStepper<QSS>::save_substep(G4double time, G4double length)
535 {
536 memcpy(CUR_SUBSTEP_X(simulator), simulator->x,
537 (Qss_misc::VAR_IDX_END * (Qss_misc::MAX_QSS_STEPPER_ORDER + 2)) * sizeof(G4double));
538
539 CUR_SUBSTEP_START(simulator) = time;
540 CUR_SUBSTEP_LEN(simulator) = length;
541 CUR_SUBSTEP(simulator)++;
542
543 if (unlikely(CUR_SUBSTEP(simulator) == MAX_SUBSTEP(simulator))) {
544 this->realloc_substeps();
545 }
546 }
547
548 template <class QSS>
549 inline void G4QSStepper<QSS>::realloc_substeps()
550 {
551 const G4int prev_index = MAX_SUBSTEP(simulator), new_index = 2 * prev_index;
552
553 MAX_SUBSTEP(simulator) = new_index;
554 SUBSTEPS(simulator) =
555 (QSSSubstep)realloc(SUBSTEPS(simulator), new_index * sizeof(*SUBSTEPS(simulator)));
556 }
557
558 template <class QSS>
559 inline void G4QSStepper<QSS>::get_state_from_poly(
560 G4double* x, G4double* tx, G4double time, G4double* state)
561 {
562 unsigned int coeff_index = 0, i;
563 const unsigned int x_order = method->order(), x_order1 = x_order + 1;
564
565 for (i = 0; i < Qss_misc::VAR_IDX_END; ++i) {
566 assert(tx[i] <= time);
567 state[i] = method->evaluate_x_poly(coeff_index, time - tx[i], x);
568 coeff_index += x_order1;
569 }
570 }
571
572 template <class QSS>
573 inline void G4QSStepper<QSS>::velocity_to_momentum(G4double* state)
574 {
575 using Qss_misc::VXidx;
576 using Qss_misc::VYidx;
577 using Qss_misc::VZidx;
578 G4double coeff = fMassGamma / cLight_local;
579
580 state[VXidx] *= coeff;
581 state[VYidx] *= coeff;
582 state[VZidx] *= coeff;
583 }
584
585 #endif
586