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