Geant4 Cross Reference |
1 // ******************************************* 1 2 // * License and Disclaimer 3 // * 4 // * The Geant4 software is copyright of th 5 // * the Geant4 Collaboration. It is provided 6 // * conditions of the Geant4 Software License 7 // * LICENSE and available at http://cern.ch/ 8 // * include a list of copyright holders. 9 // * 10 // * Neither the authors of this software syst 11 // * institutes,nor the agencies providing fin 12 // * work make any representation or warran 13 // * regarding this software system or assum 14 // * use. Please see the license in the file 15 // * for the full disclaimer and the limitatio 16 // * 17 // * This code implementation is the result 18 // * technical work of the GEANT4 collaboratio 19 // * By using, copying, modifying or distri 20 // * any work based on the software) you ag 21 // * use in resulting scientific publicati 22 // * acceptance of all terms of the Geant4 Sof 23 // ******************************************* 24 // 25 // G4BulirschStoer class implementation 26 // Based on bulirsch_stoer.hpp from boost 27 // 28 // Author: Dmitry Sorokin, Google Summer of Co 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 / STEP 44 constexpr G4double inv_STEPFAC4 = 1.0 / STEP 45 } // namespace 46 47 G4BulirschStoer::G4BulirschStoer(G4EquationOfM 48 G4int nvar, G 49 : fnvar(nvar), m_eps_rel(eps_rel), m_midpoin 50 { 51 /* initialize sequence of stage numbers and 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_seq 63 } 64 for(G4int k = 0; k < i; ++k) 65 { 66 const G4double r = static_cast<G4double> 67 / static_cast<G4double> 68 m_coeff[i][k] = 1.0 / (r * r - 1.0); // 69 } 70 71 // crude estimate of optimal order 72 m_current_k_opt = 4; 73 74 // no calculation because log10 might not 75 76 //const G4double logfact = -log10(std::max 77 //m_current_k_opt = std::max(1., 78 // std::min(static_cast<G 79 } 80 } 81 82 G4BulirschStoer::step_result 83 G4BulirschStoer::try_step( const G4double in[] 84 G4double& t, G4doub 85 { 86 if(m_max_dt < dt) 87 { 88 // given step size is bigger then max_dt s 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 97 } 98 99 G4bool reject = true; 100 101 G4double new_h = dt; 102 103 /* m_current_k_opt is the estimated current 104 105 for(G4int k = 0; k <= m_current_k_opt+1; ++k 106 { 107 // the stage counts are stored in m_interv 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] 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_ 126 h_opt[k] = calc_h_opt(dt, error, k); 127 work[k] = static_cast<G4double>(m_cost[k 128 129 if( (k == m_current_k_opt-1) || m_first) 130 { 131 if(error < 1.0) 132 { 133 // convergence 134 reject = false; 135 if( (work[k] < KFAC2 * work[k-1]) || 136 { 137 // leave order as is (except we we 138 m_current_k_opt = std::min(m_k_max 139 new_h = h_opt[k]; 140 new_h *= static_cast<G4double>(m_c 141 / static_cast<G4double>(m_c 142 } 143 else 144 { 145 m_current_k_opt = std::min(m_k_max 146 new_h = h_opt[k]; 147 } 148 break; 149 } 150 if(should_reject(error , k) && !m_firs 151 { 152 reject = true; 153 new_h = h_opt[k]; 154 break; 155 } 156 } 157 if(k == m_current_k_opt) // convergence 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_ 166 new_h = h_opt[m_current_k_opt]; 167 } 168 else if( (work[k] < KFAC2 * work[k-1 169 { 170 m_current_k_opt = std::min(m_k_max 171 new_h = h_opt[k]; 172 new_h *= static_cast<G4double>(m_c 173 / static_cast<G4double>(m_c 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) // converg 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_cu 196 } 197 if((work[k] < KFAC2 * work[m_current 198 { 199 m_current_k_opt = std::min(m_k_max 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_res 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 239 { 240 /* polynomial extrapolation, see http://www. 241 * uses the obtained intermediate results to 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. + 248 - m_table[j-1][i] * m_co 249 } 250 } 251 for (G4int i = 0; i < fnvar; ++i) 252 { 253 xest[i] = m_table[0][i] * (1. + m_coeff[k] 254 } 255 } 256 257 G4double 258 G4BulirschStoer::calc_h_opt(G4double h , G4dou 259 { 260 /* calculates the optimal step size for a gi 261 262 const G4double expo = 1.0 / (2 * k + 1); 263 const G4double facmin = std::pow(STEPFAC3, e 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_STEP 274 fac = std::max(facmin * inv_STEPFAC4, std: 275 } 276 277 return h * fac; 278 } 279 280 //why is not used!!?? 281 G4bool G4BulirschStoer::set_k_opt(std::size_t 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 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 298 { // same order - also do this if last step 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 st 304 m_current_k_opt = (G4int)k + 1; 305 dt = h_opt[m_current_k_opt - 1] * m_cost[m 306 / m_cost[m_current_k_opt - 1]; 307 return true; 308 } 309 } 310 311 G4bool G4BulirschStoer::in_convergence_window( 312 { 313 if( (k == m_current_k_opt - 1) && !m_last_st 314 { 315 return true; // decrease stepsize only if 316 } 317 return (k == m_current_k_opt) || (k == m_cur 318 } 319 320 321 G4bool G4BulirschStoer::should_reject(G4double 322 { 323 if(k == m_current_k_opt - 1) 324 { 325 const auto d = G4double(m_interval_sequen 326 * m_interval_sequence[m_curr 327 const auto e = G4double(m_interval_sequen 328 const G4double e2 = e*e; 329 // step will fail, criterion 17.3.17 in NR 330 return error * e2 * e2 > d * d; // was r 331 } 332 if(k == m_current_k_opt) 333 { 334 const auto d = G4double(m_interval_sequen 335 const auto e = G4double(m_interval_sequen 336 return error * e * e > d * d; // was retu 337 } 338 else 339 { 340 return error > 1.0; 341 } 342 } 343