diff --git a/src/core/magnetostatics/dlc.cpp b/src/core/magnetostatics/dlc.cpp index 15d820d3c8..183564e818 100644 --- a/src/core/magnetostatics/dlc.cpp +++ b/src/core/magnetostatics/dlc.cpp @@ -50,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -430,15 +431,26 @@ double DipolarLayerCorrection::tune_far_cut() const { } auto constexpr limitkc = 200; + auto constexpr exp_max = +709.8; // for IEEE-compatible double + auto constexpr exp_min = -708.4; // for IEEE-compatible double + auto const log_max = std::log(std::numeric_limits::max()); auto const piarea = std::numbers::pi / (lx * ly); auto const nmp = static_cast(count_magnetic_particles(system)); auto const h = dlc.box_h; auto far_cut = -1.; for (int kc = 1; kc < limitkc; kc++) { auto const gc = kc * 2. * std::numbers::pi / lx; - auto const fa0 = sqrt(9. * exp(+2. * gc * h) * g1_DLC_dip(gc, lz - h) + - 9. * exp(-2. * gc * h) * g1_DLC_dip(gc, lz + h) + - 22. * g1_DLC_dip(gc, lz)); + auto const pref1 = g1_DLC_dip(gc, lz - h); + auto const pref2 = g1_DLC_dip(gc, lz + h); + auto const pref0 = g1_DLC_dip(gc, lz); + auto const exp_term = 2. * gc * h; + auto const log_term1 = exp_term + std::log(9. * pref1); + if (exp_term >= exp_max or log_term1 >= log_max or -exp_term <= exp_min) { + // no solution can be found at larger k values due to overflow/underflow + break; + } + auto const fa0 = std::sqrt(9. * std::exp(+exp_term) * pref1 + + 9. * std::exp(-exp_term) * pref2 + 22. * pref0); auto const fa1 = sqrt(0.125 * piarea) * fa0; auto const fa2 = g2_DLC_dip(gc, lz); auto const de = nmp * mu_max_sq / (4. * (exp(gc * lz) - 1.)) * (fa1 + fa2); diff --git a/src/script_interface/electrostatics/CoulombP3M.hpp b/src/script_interface/electrostatics/CoulombP3M.hpp index a641a6f625..5560fe7bff 100644 --- a/src/script_interface/electrostatics/CoulombP3M.hpp +++ b/src/script_interface/electrostatics/CoulombP3M.hpp @@ -116,28 +116,35 @@ class CoulombP3M : public Actor, ::CoulombP3M> { m_tune_verbose = get_value(params, "verbose"); m_tune_limits = {std::nullopt, std::nullopt}; if (params.contains("tune_limits")) { - std::vector range; - try { - auto const val = get_value>(params, "tune_limits"); - assert(val.size() == 2u); - range.emplace_back(val[0u]); - range.emplace_back(val[1u]); - } catch (...) { - range = get_value>(params, "tune_limits"); - assert(range.size() == 2u); - } - if (not is_none(range[0u])) { - m_tune_limits.first = get_value(range[0u]); - } - if (not is_none(range[1u])) { - m_tune_limits.second = get_value(range[1u]); + auto const &variant = params.at("tune_limits"); + std::size_t range_length = 0u; + if (is_type>(variant)) { + auto const range = get_value>(variant); + range_length = range.size(); + if (range_length == 2u) { + m_tune_limits = {range[0u], range[1u]}; + } + } else { + auto const range = get_value>(variant); + range_length = range.size(); + if (range_length == 2u) { + if (not is_none(range[0u])) { + m_tune_limits.first = get_value(range[0u]); + } + if (not is_none(range[1u])) { + m_tune_limits.second = get_value(range[1u]); + } + } } context()->parallel_try_catch([&]() { + if (range_length != 2u) { + throw std::invalid_argument("Parameter 'tune_limits' needs 2 values"); + } if (m_tune_limits.first and *m_tune_limits.first <= 0) { - throw std::domain_error("P3M mesh tuning limits: mesh must be > 0"); + throw std::domain_error("Parameter 'tune_limits' must be > 0"); } if (m_tune_limits.second and *m_tune_limits.second <= 0) { - throw std::domain_error("P3M mesh tuning limits: mesh must be > 0"); + throw std::domain_error("Parameter 'tune_limits' must be > 0"); } }); } diff --git a/src/script_interface/magnetostatics/DipolarP3M.hpp b/src/script_interface/magnetostatics/DipolarP3M.hpp index 8559171a27..4436bd58c5 100644 --- a/src/script_interface/magnetostatics/DipolarP3M.hpp +++ b/src/script_interface/magnetostatics/DipolarP3M.hpp @@ -111,28 +111,35 @@ class DipolarP3M : public Actor, ::DipolarP3M> { m_tune_verbose = get_value(params, "verbose"); m_tune_limits = {std::nullopt, std::nullopt}; if (params.contains("tune_limits")) { - std::vector range; - try { - auto const val = get_value>(params, "tune_limits"); - assert(val.size() == 2u); - range.emplace_back(val[0u]); - range.emplace_back(val[1u]); - } catch (...) { - range = get_value>(params, "tune_limits"); - assert(range.size() == 2u); - } - if (not is_none(range[0u])) { - m_tune_limits.first = get_value(range[0u]); - } - if (not is_none(range[1u])) { - m_tune_limits.second = get_value(range[1u]); + auto const &variant = params.at("tune_limits"); + std::size_t range_length = 0u; + if (is_type>(variant)) { + auto const range = get_value>(variant); + range_length = range.size(); + if (range_length == 2u) { + m_tune_limits = {range[0u], range[1u]}; + } + } else { + auto const range = get_value>(variant); + range_length = range.size(); + if (range_length == 2u) { + if (not is_none(range[0u])) { + m_tune_limits.first = get_value(range[0u]); + } + if (not is_none(range[1u])) { + m_tune_limits.second = get_value(range[1u]); + } + } } context()->parallel_try_catch([&]() { + if (range_length != 2u) { + throw std::invalid_argument("Parameter 'tune_limits' needs 2 values"); + } if (m_tune_limits.first and *m_tune_limits.first <= 0) { - throw std::domain_error("P3M mesh tuning limits: mesh must be > 0"); + throw std::domain_error("Parameter 'tune_limits' must be > 0"); } if (m_tune_limits.second and *m_tune_limits.second <= 0) { - throw std::domain_error("P3M mesh tuning limits: mesh must be > 0"); + throw std::domain_error("Parameter 'tune_limits' must be > 0"); } }); } diff --git a/testsuite/python/p3m_tuning_exceptions.py b/testsuite/python/p3m_tuning_exceptions.py index 8ccfa6f569..e7aa4e82dc 100644 --- a/testsuite/python/p3m_tuning_exceptions.py +++ b/testsuite/python/p3m_tuning_exceptions.py @@ -205,8 +205,9 @@ def check_invalid_params(self, container, class_solver, **custom_params): ('alpha', -2.0, "Parameter 'alpha' must be > 0"), ('accuracy', -2.0, "Parameter 'accuracy' must be > 0"), ('mesh', (-1, -1, -1), "Parameter 'mesh' must be > 0"), - ('tune_limits', (-1, 1), "P3M mesh tuning limits: mesh must be > 0"), - ('tune_limits', (1, 0), "P3M mesh tuning limits: mesh must be > 0"), + ('tune_limits', (-1,), "Parameter 'tune_limits' needs 2 values"), + ('tune_limits', (-1, 1), "Parameter 'tune_limits' must be > 0"), + ('tune_limits', (1, 0), "Parameter 'tune_limits' must be > 0"), ('mesh', (2, 2, 2), "Parameter 'cao' cannot be larger than 'mesh'"), ('mesh_off', (-2, 1, 1), "Parameter 'mesh_off' must be >= 0 and <= 1"), ]