Commit 22359d1f by Timo Koch Committed by Dennis Gläser

### [test][math] Lower precision for polynomial test

```When the roots are very close together the inversion becomes ill-conditioned.
Exclude these cases from the random generation, test some corner cases
with close roots and decrease the tolerance for root comparison.```
parent 2cfc0539
Pipeline #5003 waiting for manual action with stages
 ... @@ -265,8 +265,8 @@ int main() ... @@ -265,8 +265,8 @@ int main() std::uniform_real_distribution<> dis(-10.0, 10.0); std::uniform_real_distribution<> dis(-10.0, 10.0); auto threeRandomNumbers = [&]() mutable { auto threeRandomNumbers = [&]() mutable { const auto x = dis(gen); const auto x = dis(gen); auto y = dis(gen); while (Dune::FloatCmp::eq(x, y, 1e-10)) { y = dis(gen); } auto y = dis(gen); while (Dune::FloatCmp::eq(x, y, 1e-3)) { y = dis(gen); } auto z = dis(gen); while (Dune::FloatCmp::eq(x, z, 1e-10) || Dune::FloatCmp::eq(y, z, 1e-10)) { z = dis(gen); } auto z = dis(gen); while (Dune::FloatCmp::eq(x, z, 1e-3) || Dune::FloatCmp::eq(y, z, 1e-3)) { z = dis(gen); } return std::make_tuple(x, y, z); return std::make_tuple(x, y, z); }; }; ... @@ -287,7 +287,7 @@ int main() ... @@ -287,7 +287,7 @@ int main() using std::isfinite; using std::isfinite; if (!isfinite(roots[i])) if (!isfinite(roots[i])) DUNE_THROW(Dune::Exception, "Expecting finite root"); DUNE_THROW(Dune::Exception, "Expecting finite root"); if (throwing && std::none_of(ref.begin(), ref.end(), [&](auto r){ return Dune::FloatCmp::eq(roots[i], r, 1e-8); })) if (throwing && std::none_of(ref.begin(), ref.end(), [&](auto r){ return Dune::FloatCmp::eq(roots[i], r, 1e-5); })) DUNE_THROW(Dune::Exception, "[invertCubicPolynomial] Root " << Dumux::Fmt::format( DUNE_THROW(Dune::Exception, "[invertCubicPolynomial] Root " << Dumux::Fmt::format( "{} of {}: {:.14e}, does not match reference roots = [{}].\n", i+1, n, roots[i], refToString()) "{} of {}: {:.14e}, does not match reference roots = [{}].\n", i+1, n, roots[i], refToString()) << Dumux::Fmt::format("Polynomial: {:+.5e}x³ {:+.5e}x² {:+.5e}x {:+.5e}\n", coeff[0], coeff[1], coeff[2], coeff[3])); << Dumux::Fmt::format("Polynomial: {:+.5e}x³ {:+.5e}x² {:+.5e}x {:+.5e}\n", coeff[0], coeff[1], coeff[2], coeff[3])); ... @@ -301,6 +301,7 @@ int main() ... @@ -301,6 +301,7 @@ int main() std::array roots{}; std::array roots{}; for (double scaling : {1.0, -1.0, 5.0}) for (double scaling : {1.0, -1.0, 5.0}) { { // random numbers for (int i = 0; i < 100000; ++i) for (int i = 0; i < 100000; ++i) { { const auto [x, y, z] = threeRandomNumbers(); const auto [x, y, z] = threeRandomNumbers(); ... @@ -347,6 +348,33 @@ int main() ... @@ -347,6 +348,33 @@ int main() } } } } } } // some corner cases with close roots { const auto [x, y, z] = std::make_tuple(10.0, 10.0+1e-8, 0.0); const auto [a, b, c, d] = cubicCoefficientsFromRoots(x, y, z, scaling); int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d); if (numRoots != 3) DUNE_THROW(Dune::Exception, "Expecting three roots! Got " << numRoots); if (!checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots, false)) { // Try to increase Newton iterations for increased precision and try again int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d, 10); checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots); } }{ const auto [x, y, z] = std::make_tuple(10.0, 10.0+1e-4, 10.0+2e-4); const auto [a, b, c, d] = cubicCoefficientsFromRoots(x, y, z, scaling); int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d); if (numRoots != 3) DUNE_THROW(Dune::Exception, "Expecting three roots! Got " << numRoots); if (!checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots, false)) { // Try to increase Newton iterations for increased precision and try again int numRoots = Dumux::invertCubicPolynomial(roots.data(), a, b, c, d, 10); checkRoots(roots, {x, y, z}, {a, b, c, d}, numRoots); } } } } // Test one corner case (see MR!2594) // Test one corner case (see MR!2594) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!