crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Bug fix in quartic solver when d2 is close to zero

Browse source 
This bug was resulted in lost particles for torus models with R ~ r. Bug was
reported to author of quartic solver and a fix was made in the quarticpp
project: https://github.com/cridemichel/quarticpp/issues/1
This commit is contained in:
Paul Romano 2021-12-29 09:55:54 -05:00
parent d497d54966
commit 9a0f56a891
1 changed files with 1 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/external/quartic_solver.cpp vendored
View file

@ -11,11 +11,6 @@ constexpr double CUBIC_RESCAL_FACT = 3.488062113727083e+102;
constexpr double QUART_RESCAL_FACT = 7.156344627944542e+76;
constexpr double MACHEPS = std::numeric_limits<double>::epsilon();
double max3(double a, double b, double c)
{
return std::max(std::max(a, b), c);
}
double solve_cubic_analytic_depressed_handle_inf(double b, double c)
{
/* find analytically the dominant root of a depressed cubic x^3+b*x+c
@ -479,7 +474,7 @@ void quartic_solver(double coeff[5], std::complex<double> roots[4])
* better) */
if (realcase[0] == -1 ||
(std::abs(d2) <=
MACHEPS * oqs::max3(std::abs(2. * b / 3.), std::abs(phi0), l1 * l1))) {
MACHEPS * (std::abs(2. * b / 3.) + std::abs(phi0) + l1 * l1))) {
double d3 = d - l3 * l3;
double err0 = 0.0;
if (realcase[0] == 1)