Geant4 Cross Reference |
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 47 G4BulirschStoer::G4BulirschStoer(G4EquationOfMotion* equation, 48 G4int nvar, G4double eps_rel, G4double max_dt) 49 : fnvar(nvar), m_eps_rel(eps_rel), m_midpoint(equation,nvar), m_max_dt(max_dt) 50 { 51 /* initialize sequence of stage numbers and work */ 52 53 for(G4int i = 0; i < m_k_max + 1; ++i) 54 { 55 m_interval_sequence[i] = 2 * (i + 1); 56 if (i == 0) 57 { 58 m_cost[i] = m_interval_sequence[i]; 59 } 60 else 61 { 62 m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; 63 } 64 for(G4int k = 0; k < i; ++k) 65 { 66 const G4double r = static_cast<G4double>(m_interval_sequence[i]) 67 / static_cast<G4double>(m_interval_sequence[k]); 68 m_coeff[i][k] = 1.0 / (r * r - 1.0); // coefficients for extrapolation 69 } 70 71 // crude estimate of optimal order 72 m_current_k_opt = 4; 73 74 // no calculation because log10 might not exist for value_type! 75 76 //const G4double logfact = -log10(std::max(eps_rel, 1.0e-12)) * 0.6 + 0.5; 77 //m_current_k_opt = std::max(1., 78 // std::min(static_cast<G4double>(m_k_max-1), logfact)); 79 } 80 } 81 82 G4BulirschStoer::step_result 83 G4BulirschStoer::try_step( const G4double in[], const G4double dxdt[], 84 G4double& t, G4double out[], G4double& dt) 85 { 86 if(m_max_dt < dt) 87 { 88 // given step size is bigger then max_dt set limit and return fail 89 // 90 dt = m_max_dt; 91 return step_result::fail; 92 } 93 94 if (dt != m_dt_last) 95 { 96 reset(); // step size changed from outside -> reset 97 } 98 99 G4bool reject = true; 100 101 G4double new_h = dt; 102 103 /* m_current_k_opt is the estimated current optimal stage number */ 104 105 for(G4int k = 0; k <= m_current_k_opt+1; ++k) 106 { 107 // the stage counts are stored in m_interval_sequence 108 // 109 m_midpoint.SetSteps(m_interval_sequence[k]); 110 if(k == 0) 111 { 112 m_midpoint.DoStep(in, dxdt, out, dt); 113 /* the first step, nothing more to do */ 114 } 115 else 116 { 117 m_midpoint.DoStep(in, dxdt, m_table[k-1], dt); 118 extrapolate(k, out); 119 // get error estimate 120 for (G4int i = 0; i < fnvar; ++i) 121 { 122 m_err[i] = out[i] - m_table[0][i]; 123 } 124 const G4double error = 125 field_utils::relativeError(out, m_err, dt, m_eps_rel); 126 h_opt[k] = calc_h_opt(dt, error, k); 127 work[k] = static_cast<G4double>(m_cost[k]) / h_opt[k]; 128 129 if( (k == m_current_k_opt-1) || m_first) // convergence before k_opt ? 130 { 131 if(error < 1.0) 132 { 133 // convergence 134 reject = false; 135 if( (work[k] < KFAC2 * work[k-1]) || (m_current_k_opt <= 2) ) 136 { 137 // leave order as is (except we were in first round) 138 m_current_k_opt = std::min(m_k_max - 1 , std::max(2 , k + 1)); 139 new_h = h_opt[k]; 140 new_h *= static_cast<G4double>(m_cost[k + 1]) 141 / static_cast<G4double>(m_cost[k]); 142 } 143 else 144 { 145 m_current_k_opt = std::min(m_k_max - 1, std::max(2, k)); 146 new_h = h_opt[k]; 147 } 148 break; 149 } 150 if(should_reject(error , k) && !m_first) 151 { 152 reject = true; 153 new_h = h_opt[k]; 154 break; 155 } 156 } 157 if(k == m_current_k_opt) // convergence at k_opt ? 158 { 159 if(error < 1.0) 160 { 161 // convergence 162 reject = false; 163 if(work[k-1] < KFAC2 * work[k]) 164 { 165 m_current_k_opt = std::max( 2 , m_current_k_opt-1 ); 166 new_h = h_opt[m_current_k_opt]; 167 } 168 else if( (work[k] < KFAC2 * work[k-1]) && !m_last_step_rejected ) 169 { 170 m_current_k_opt = std::min(m_k_max - 1, m_current_k_opt + 1); 171 new_h = h_opt[k]; 172 new_h *= static_cast<G4double>(m_cost[m_current_k_opt]) 173 / static_cast<G4double>(m_cost[k]); 174 } 175 else 176 { 177 new_h = h_opt[m_current_k_opt]; 178 } 179 break; 180 } 181 if(should_reject(error, k)) 182 { 183 reject = true; 184 new_h = h_opt[m_current_k_opt]; 185 break; 186 } 187 } 188 if(k == m_current_k_opt + 1) // convergence at k_opt+1 ? 189 { 190 if(error < 1.0) // convergence 191 { 192 reject = false; 193 if(work[k-2] < KFAC2 * work[k-1]) 194 { 195 m_current_k_opt = std::max(2, m_current_k_opt - 1); 196 } 197 if((work[k] < KFAC2 * work[m_current_k_opt]) && !m_last_step_rejected) 198 { 199 m_current_k_opt = std::min(m_k_max - 1 , k); 200 } 201 new_h = h_opt[m_current_k_opt]; 202 } 203 else 204 { 205 reject = true; 206 new_h = h_opt[m_current_k_opt]; 207 } 208 break; 209 } 210 } 211 } 212 213 if(!reject) 214 { 215 t += dt; 216 } 217 218 if(!m_last_step_rejected || new_h < dt) 219 { 220 // limit step size 221 new_h = std::min(m_max_dt, new_h); 222 m_dt_last = new_h; 223 dt = new_h; 224 } 225 226 m_last_step_rejected = reject; 227 m_first = false; 228 229 return reject ? step_result::fail : step_result::success; 230 } 231 232 void G4BulirschStoer::reset() 233 { 234 m_first = true; 235 m_last_step_rejected = false; 236 } 237 238 void G4BulirschStoer::extrapolate(std::size_t k , G4double xest[]) 239 { 240 /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf 241 * uses the obtained intermediate results to extrapolate to dt->0 */ 242 243 for(std::size_t j = k - 1 ; j > 0; --j) 244 { 245 for (G4int i = 0; i < fnvar; ++i) 246 { 247 m_table[j-1][i] = m_table[j][i] * (1. + m_coeff[k][j]) 248 - m_table[j-1][i] * m_coeff[k][j]; 249 } 250 } 251 for (G4int i = 0; i < fnvar; ++i) 252 { 253 xest[i] = m_table[0][i] * (1. + m_coeff[k][0]) - xest[i] * m_coeff[k][0]; 254 } 255 } 256 257 G4double 258 G4BulirschStoer::calc_h_opt(G4double h , G4double error , std::size_t k) const 259 { 260 /* calculates the optimal step size for a given error and stage number */ 261 262 const G4double expo = 1.0 / (2 * k + 1); 263 const G4double facmin = std::pow(STEPFAC3, expo); 264 G4double fac; 265 266 G4double facminInv= 1.0 / facmin; 267 if (error == 0.0) 268 { 269 fac = facminInv; 270 } 271 else 272 { 273 fac = STEPFAC2 * std::pow(error * inv_STEPFAC1 , -expo); 274 fac = std::max(facmin * inv_STEPFAC4, std::min( facminInv, fac)); 275 } 276 277 return h * fac; 278 } 279 280 //why is not used!!?? 281 G4bool G4BulirschStoer::set_k_opt(std::size_t k, G4double& dt) 282 { 283 /* calculates the optimal stage number */ 284 285 if(k == 1) 286 { 287 m_current_k_opt = 2; 288 return true; 289 } 290 if( (work[k-1] < KFAC1 * work[k]) || (k == m_k_max) ) // order decrease 291 { 292 m_current_k_opt = (G4int)k - 1; 293 dt = h_opt[ m_current_k_opt ]; 294 return true; 295 } 296 if( (work[k] < KFAC2 * work[k-1]) 297 || m_last_step_rejected || (k == m_k_max-1) ) 298 { // same order - also do this if last step got rejected 299 m_current_k_opt = (G4int)k; 300 dt = h_opt[m_current_k_opt]; 301 return true; 302 } 303 else { // order increase - only if last step was not rejected 304 m_current_k_opt = (G4int)k + 1; 305 dt = h_opt[m_current_k_opt - 1] * m_cost[m_current_k_opt] 306 / m_cost[m_current_k_opt - 1]; 307 return true; 308 } 309 } 310 311 G4bool G4BulirschStoer::in_convergence_window(G4int k) const 312 { 313 if( (k == m_current_k_opt - 1) && !m_last_step_rejected ) 314 { 315 return true; // decrease stepsize only if last step was not rejected 316 } 317 return (k == m_current_k_opt) || (k == m_current_k_opt + 1); 318 } 319 320 321 G4bool G4BulirschStoer::should_reject(G4double error, G4int k) const 322 { 323 if(k == m_current_k_opt - 1) 324 { 325 const auto d = G4double(m_interval_sequence[m_current_k_opt] 326 * m_interval_sequence[m_current_k_opt+1]); 327 const auto e = G4double(m_interval_sequence[0]); 328 const G4double e2 = e*e; 329 // step will fail, criterion 17.3.17 in NR 330 return error * e2 * e2 > d * d; // was return error > dOld * dOld; (where dOld= d/e; ) 331 } 332 if(k == m_current_k_opt) 333 { 334 const auto d = G4double(m_interval_sequence[m_current_k_opt]); 335 const auto e = G4double(m_interval_sequence[0]); 336 return error * e * e > d * d; // was return error > dOld * dOld; (where dOld= d/e; ) 337 } 338 else 339 { 340 return error > 1.0; 341 } 342 } 343