ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BulirschStoer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BulirschStoer.cc
1 // ********************************************************************
2 // * License and Disclaimer *
3 // * *
4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitation of liability. *
16 // * *
17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Software license. *
23 // ********************************************************************
24 //
25 // G4BulirschStoer class implementation
26 // Based on bulirsch_stoer.hpp from boost
27 //
28 // Author: Dmitry Sorokin, Google Summer of Code 2016
29 // --------------------------------------------------------------------
30 
31 #include "G4BulirschStoer.hh"
32 
33 #include "G4FieldUtils.hh"
34 
35 namespace
36 {
37  constexpr G4double STEPFAC1 = 0.65;
38  constexpr G4double STEPFAC2 = 0.94;
39  constexpr G4double STEPFAC3 = 0.02;
40  constexpr G4double STEPFAC4 = 4.0;
41  constexpr G4double KFAC1 = 0.8;
42  constexpr G4double KFAC2 = 0.9;
43  constexpr G4double inv_STEPFAC1 = 1.0 / STEPFAC1;
44  constexpr G4double inv_STEPFAC4 = 1.0 / STEPFAC4;
45 } // namespace
46 
48  G4int nvar, G4double eps_rel, G4double max_dt)
49  : fnvar(nvar), m_eps_rel(eps_rel), m_midpoint(equation,nvar),
50  m_last_step_rejected(false), m_first(true), m_dt_last(0.0), m_max_dt(max_dt)
51 {
52  /* initialize sequence of stage numbers and work */
53 
54  for(G4int i = 0; i < m_k_max + 1; ++i)
55  {
56  m_interval_sequence[i] = 2 * (i + 1);
57  if (i == 0)
58  {
60  }
61  else
62  {
63  m_cost[i] = m_cost[i-1] + m_interval_sequence[i];
64  }
65  for(G4int k = 0; k < i; ++k)
66  {
67  const G4double r = static_cast<G4double>(m_interval_sequence[i])
68  / static_cast<G4double>(m_interval_sequence[k]);
69  m_coeff[i][k] = 1.0 / (r * r - 1.0); // coefficients for extrapolation
70  }
71 
72  // crude estimate of optimal order
73  m_current_k_opt = 4;
74 
75  // no calculation because log10 might not exist for value_type!
76 
77  //const G4double logfact = -log10(std::max(eps_rel, 1.0e-12)) * 0.6 + 0.5;
78  //m_current_k_opt = std::max(1.,
79  // std::min(static_cast<G4double>(m_k_max-1), logfact));
80  }
81 }
82 
85  G4double& t, G4double out[], G4double& dt)
86 {
87  if(m_max_dt < dt)
88  {
89  // given step size is bigger then max_dt set limit and return fail
90  //
91  dt = m_max_dt;
92  return step_result::fail;
93  }
94 
95  if (dt != m_dt_last)
96  {
97  reset(); // step size changed from outside -> reset
98  }
99 
100  G4bool reject = true;
101 
102  G4double new_h = dt;
103 
104  /* m_current_k_opt is the estimated current optimal stage number */
105 
106  for(G4int k = 0; k <= m_current_k_opt+1; ++k)
107  {
108  // the stage counts are stored in m_interval_sequence
109  //
111  if(k == 0)
112  {
113  m_midpoint.DoStep(in, dxdt, out, dt);
114  /* the first step, nothing more to do */
115  }
116  else
117  {
118  m_midpoint.DoStep(in, dxdt, m_table[k-1], dt);
119  extrapolate(k, out);
120  // get error estimate
121  for (G4int i = 0; i < fnvar; ++i)
122  {
123  m_err[i] = out[i] - m_table[0][i];
124  }
125  const G4double error =
127  h_opt[k] = calc_h_opt(dt, error, k);
128  work[k] = static_cast<G4double>(m_cost[k]) / h_opt[k];
129 
130  if( (k == m_current_k_opt-1) || m_first) // convergence before k_opt ?
131  {
132  if(error < 1.0)
133  {
134  // convergence
135  reject = false;
136  if( (work[k] < KFAC2 * work[k-1]) || (m_current_k_opt <= 2) )
137  {
138  // leave order as is (except we were in first round)
139  m_current_k_opt = std::min(m_k_max - 1 , std::max(2 , k + 1));
140  new_h = h_opt[k];
141  new_h *= static_cast<G4double>(m_cost[k + 1])
142  / static_cast<G4double>(m_cost[k]);
143  }
144  else
145  {
147  new_h = h_opt[k];
148  }
149  break;
150  }
151  else if(should_reject(error , k) && !m_first)
152  {
153  reject = true;
154  new_h = h_opt[k];
155  break;
156  }
157  }
158  if(k == m_current_k_opt) // convergence at k_opt ?
159  {
160  if(error < 1.0)
161  {
162  // convergence
163  reject = false;
164  if(work[k-1] < KFAC2 * work[k])
165  {
167  new_h = h_opt[m_current_k_opt];
168  }
169  else if( (work[k] < KFAC2 * work[k-1]) && !m_last_step_rejected )
170  {
172  new_h = h_opt[k];
173  new_h *= static_cast<G4double>(m_cost[m_current_k_opt])
174  / static_cast<G4double>(m_cost[k]);
175  }
176  else
177  {
178  new_h = h_opt[m_current_k_opt];
179  }
180  break;
181  }
182  else if(should_reject(error, k))
183  {
184  reject = true;
185  new_h = h_opt[m_current_k_opt];
186  break;
187  }
188  }
189  if(k == m_current_k_opt + 1) // convergence at k_opt+1 ?
190  {
191  if(error < 1.0) // convergence
192  {
193  reject = false;
194  if(work[k-2] < KFAC2 * work[k-1])
195  {
197  }
198  if((work[k] < KFAC2 * work[m_current_k_opt]) && !m_last_step_rejected)
199  {
200  m_current_k_opt = std::min(m_k_max - 1 , k);
201  }
202  new_h = h_opt[m_current_k_opt];
203  }
204  else
205  {
206  reject = true;
207  new_h = h_opt[m_current_k_opt];
208  }
209  break;
210  }
211  }
212  }
213 
214  if(!reject)
215  {
216  t += dt;
217  }
218 
219  if(!m_last_step_rejected || new_h < dt)
220  {
221  // limit step size
222  new_h = std::min(m_max_dt, new_h);
223  m_dt_last = new_h;
224  dt = new_h;
225  }
226 
227  m_last_step_rejected = reject;
228  m_first = false;
229 
230  return reject ? step_result::fail : step_result::success;
231 }
232 
234 {
235  m_first = true;
236  m_last_step_rejected = false;
237 }
238 
240 {
241  /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
242  * uses the obtained intermediate results to extrapolate to dt->0 */
243 
244  for(G4int j = k - 1 ; j > 0; --j)
245  {
246  for (G4int i = 0; i < fnvar; ++i)
247  {
248  m_table[j-1][i] = m_table[j][i] * (1. + m_coeff[k][j])
249  - m_table[j-1][i] * m_coeff[k][j];
250  }
251  }
252  for (G4int i = 0; i < fnvar; ++i)
253  {
254  xest[i] = m_table[0][i] * (1. + m_coeff[k][0]) - xest[i] * m_coeff[k][0];
255  }
256 }
257 
258 G4double
260 {
261  /* calculates the optimal step size for a given error and stage number */
262 
263  const G4double expo = 1.0 / (2 * k + 1);
264  const G4double facmin = std::pow(STEPFAC3, expo);
265  G4double fac;
266 
267  G4double facminInv= 1.0 / facmin;
268  if (error == 0.0)
269  {
270  fac = facminInv;
271  }
272  else
273  {
274  fac = STEPFAC2 * std::pow(error * inv_STEPFAC1 , -expo);
275  fac = std::max(facmin * inv_STEPFAC4, std::min( facminInv, fac));
276  }
277 
278  return h * fac;
279 }
280 
281 //why is not used!!??
283 {
284  /* calculates the optimal stage number */
285 
286  if(k == 1)
287  {
288  m_current_k_opt = 2;
289  return true;
290  }
291  if( (work[k-1] < KFAC1 * work[k]) || (k == m_k_max) ) // order decrease
292  {
293  m_current_k_opt = k - 1;
294  dt = h_opt[ m_current_k_opt ];
295  return true;
296  }
297  else if( (work[k] < KFAC2 * work[k-1])
298  || m_last_step_rejected || (k == m_k_max-1) )
299  { // same order - also do this if last step got rejected
300  m_current_k_opt = k;
301  dt = h_opt[m_current_k_opt];
302  return true;
303  }
304  else { // order increase - only if last step was not rejected
305  m_current_k_opt = k + 1;
307  / m_cost[m_current_k_opt - 1];
308  return true;
309  }
310 }
311 
313 {
314  if( (k == m_current_k_opt - 1) && !m_last_step_rejected )
315  {
316  return true; // decrease stepsize only if last step was not rejected
317  }
318  return (k == m_current_k_opt) || (k == m_current_k_opt + 1);
319 }
320 
321 
323 {
324  if(k == m_current_k_opt - 1)
325  {
329  const G4double e2 = e*e;
330  // step will fail, criterion 17.3.17 in NR
331  return error * e2 * e2 > d * d; // was return error > dOld * dOld; (where dOld= d/e; )
332  }
333  else if(k == m_current_k_opt)
334  {
337  return error * e * e > d * d; // was return error > dOld * dOld; (where dOld= d/e; )
338  }
339  else
340  {
341  return error > 1.0;
342  }
343 }