Line data Source code
1 : #include "crpropa/module/PropagationBP.h"
2 :
3 : #include <sstream>
4 : #include <stdexcept>
5 : #include <vector>
6 :
7 : namespace crpropa {
8 28 : void PropagationBP::tryStep(const Y &y, Y &out, Y &error, double h,
9 : ParticleState &particle, double z, double q, double m) const {
10 28 : out = dY(y.x, y.u, h, z, q, m); // 1 step with h
11 :
12 28 : Y outHelp = dY(y.x, y.u, h/2, z, q, m); // 2 steps with h/2
13 28 : Y outCompare = dY(outHelp.x, outHelp.u, h/2, z, q, m);
14 :
15 28 : error = errorEstimation(out.x , outCompare.x , h);
16 28 : }
17 :
18 :
19 103 : PropagationBP::Y PropagationBP::dY(Vector3d pos, Vector3d dir, double step,
20 : double z, double q, double m) const {
21 : // half leap frog step in the position
22 : pos += dir * step / 2.;
23 :
24 : // get B field at particle position
25 103 : Vector3d B = getFieldAtPosition(pos, z);
26 :
27 : // Boris help vectors
28 : Vector3d t = B * q / 2 / m * step / c_light;
29 103 : Vector3d s = t * 2 / (1 + t.dot(t));
30 : Vector3d v_help;
31 :
32 : // Boris push
33 : v_help = dir + dir.cross(t);
34 : dir = dir + v_help.cross(s);
35 :
36 : // the other half leap frog step in the position
37 : pos += dir * step / 2.;
38 103 : return Y(pos, dir);
39 : }
40 :
41 :
42 : // with a fixed step size
43 10 : PropagationBP::PropagationBP(ref_ptr<MagneticField> field, double fixedStep) :
44 10 : minStep(0) {
45 10 : setField(field);
46 10 : setTolerance(0.42);
47 10 : setMaximumStep(fixedStep);
48 10 : setMinimumStep(fixedStep);
49 10 : }
50 :
51 :
52 : // with adaptive step size
53 5 : PropagationBP::PropagationBP(ref_ptr<MagneticField> field, double tolerance, double minStep, double maxStep) :
54 5 : minStep(0) {
55 7 : setField(field);
56 5 : setTolerance(tolerance);
57 4 : setMaximumStep(maxStep);
58 4 : setMinimumStep(minStep);
59 3 : }
60 :
61 :
62 37 : void PropagationBP::process(Candidate *candidate) const {
63 : // save the new previous particle state
64 37 : ParticleState ¤t = candidate->current;
65 : candidate->previous = current;
66 :
67 37 : Y yIn(current.getPosition(), current.getDirection());
68 :
69 : // calculate charge of particle
70 37 : double q = current.getCharge();
71 37 : double step = maxStep;
72 :
73 : // rectilinear propagation for neutral particles
74 37 : if (q == 0) {
75 2 : step = clip(candidate->getNextStep(), minStep, maxStep);
76 1 : current.setPosition(yIn.x + yIn.u * step);
77 1 : candidate->setCurrentStep(step);
78 1 : candidate->setNextStep(maxStep);
79 : return;
80 : }
81 :
82 : Y yOut, yErr;
83 36 : double newStep = step;
84 36 : double z = candidate->getRedshift();
85 36 : double m = current.getEnergy()/(c_light * c_light);
86 :
87 : // if minStep is the same as maxStep the adaptive algorithm with its error
88 : // estimation is not needed and the computation time can be saved:
89 36 : if (minStep == maxStep){
90 19 : yOut = dY(yIn.x, yIn.u, step, z, q, m);
91 : } else {
92 17 : step = clip(candidate->getNextStep(), minStep, maxStep);
93 17 : newStep = step;
94 : double r = 42; // arbitrary value
95 :
96 : // try performing step until the target error (tolerance) or the minimum/maximum step size has been reached
97 : while (true) {
98 28 : tryStep(yIn, yOut, yErr, step, current, z, q, m);
99 28 : r = yErr.u.getR() / tolerance; // ratio of absolute direction error and tolerance
100 28 : if (r > 1) { // large direction error relative to tolerance, try to decrease step size
101 18 : if (step == minStep) // already minimum step size
102 : break;
103 : else {
104 11 : newStep = step * 0.95 * pow(r, -0.2);
105 19 : newStep = std::max(newStep, 0.1 * step); // limit step size decrease
106 11 : newStep = std::max(newStep, minStep); // limit step size to minStep
107 : step = newStep;
108 : }
109 : } else { // small direction error relative to tolerance, try to increase step size
110 10 : if (step != maxStep) { // only update once if maximum step size yet not reached
111 10 : newStep = step * 0.95 * pow(r, -0.2);
112 17 : newStep = std::min(newStep, 5 * step); // limit step size increase
113 10 : newStep = std::min(newStep, maxStep); // limit step size to maxStep
114 : }
115 : break;
116 : }
117 : }
118 : }
119 :
120 36 : current.setPosition(yOut.x);
121 36 : current.setDirection(yOut.u.getUnitVector());
122 36 : candidate->setCurrentStep(step);
123 36 : candidate->setNextStep(newStep);
124 : }
125 :
126 :
127 16 : void PropagationBP::setField(ref_ptr<MagneticField> f) {
128 16 : field = f;
129 16 : }
130 :
131 :
132 2 : ref_ptr<MagneticField> PropagationBP::getField() const {
133 2 : return field;
134 : }
135 :
136 :
137 104 : Vector3d PropagationBP::getFieldAtPosition(Vector3d pos, double z) const {
138 : Vector3d B(0, 0, 0);
139 : try {
140 : // check if field is valid and use the field vector at the
141 : // position pos with the redshift z
142 104 : if (field.valid())
143 104 : B = field->getField(pos, z);
144 0 : } catch (std::exception &e) {
145 0 : KISS_LOG_ERROR << "PropagationBP: Exception in PropagationBP::getFieldAtPosition.\n"
146 0 : << e.what();
147 0 : }
148 104 : return B;
149 : }
150 :
151 :
152 28 : double PropagationBP::errorEstimation(const Vector3d x1, const Vector3d x2, double step) const {
153 : // compare the position after one step with the position after two steps with step/2.
154 : Vector3d diff = (x1 - x2);
155 :
156 28 : double S = diff.getR() / (step * (1 - 1/4.) ); // 1/4 = (1/2)² number of steps for x1 divided by number of steps for x2 to the power of p (order)
157 :
158 28 : return S;
159 : }
160 :
161 :
162 22 : void PropagationBP::setTolerance(double tol) {
163 22 : if ((tol > 1) or (tol < 0))
164 : throw std::runtime_error(
165 2 : "PropagationBP: target error not in range 0-1");
166 20 : tolerance = tol;
167 20 : }
168 :
169 :
170 24 : void PropagationBP::setMinimumStep(double min) {
171 24 : if (min < 0)
172 1 : throw std::runtime_error("PropagationBP: minStep < 0 ");
173 23 : if (min > maxStep)
174 2 : throw std::runtime_error("PropagationBP: minStep > maxStep");
175 21 : minStep = min;
176 21 : }
177 :
178 :
179 22 : void PropagationBP::setMaximumStep(double max) {
180 22 : if (max < minStep)
181 1 : throw std::runtime_error("PropagationBP: maxStep < minStep");
182 21 : maxStep = max;
183 21 : }
184 :
185 :
186 2 : double PropagationBP::getTolerance() const {
187 2 : return tolerance;
188 : }
189 :
190 :
191 2 : double PropagationBP::getMinimumStep() const {
192 2 : return minStep;
193 : }
194 :
195 :
196 5 : double PropagationBP::getMaximumStep() const {
197 5 : return maxStep;
198 : }
199 :
200 :
201 0 : std::string PropagationBP::getDescription() const {
202 0 : std::stringstream s;
203 0 : s << "Propagation in magnetic fields using the adaptive Boris push method.";
204 0 : s << " Target error: " << tolerance;
205 0 : s << ", Minimum Step: " << minStep / kpc << " kpc";
206 0 : s << ", Maximum Step: " << maxStep / kpc << " kpc";
207 0 : return s.str();
208 0 : }
209 : } // namespace crpropa
|
