Skip to content

Commit

Permalink
Skip infinite values in DLC parameter tuning
Browse files Browse the repository at this point in the history
  • Loading branch information
@jngrad
jngrad committed Dec 12, 2024
1 parent f7f164e commit 75f0812
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 39 deletions.
18 changes: 15 additions & 3 deletions src/core/magnetostatics/dlc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@
#include <cmath>
#include <cstdio>
#include <functional>
#include <limits>
#include <numbers>
#include <numeric>
#include <stdexcept>
Expand Down Expand Up @@ -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<double>::max());
auto const piarea = std::numbers::pi / (lx * ly);
auto const nmp = static_cast<double>(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);
Expand Down
41 changes: 24 additions & 17 deletions src/script_interface/electrostatics/CoulombP3M.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,28 +116,35 @@ class CoulombP3M : public Actor<CoulombP3M<Architecture>, ::CoulombP3M> {
m_tune_verbose = get_value<bool>(params, "verbose");
m_tune_limits = {std::nullopt, std::nullopt};
if (params.contains("tune_limits")) {
std::vector<Variant> range;
try {
auto const val = get_value<std::vector<int>>(params, "tune_limits");
assert(val.size() == 2u);
range.emplace_back(val[0u]);
range.emplace_back(val[1u]);
} catch (...) {
range = get_value<std::vector<Variant>>(params, "tune_limits");
assert(range.size() == 2u);
}
if (not is_none(range[0u])) {
m_tune_limits.first = get_value<int>(range[0u]);
}
if (not is_none(range[1u])) {
m_tune_limits.second = get_value<int>(range[1u]);
auto const &variant = params.at("tune_limits");
std::size_t range_length = 0u;
if (is_type<std::vector<int>>(variant)) {
auto const range = get_value<std::vector<int>>(variant);
range_length = range.size();
if (range_length == 2u) {
m_tune_limits = {range[0u], range[1u]};
}
} else {
auto const range = get_value<std::vector<Variant>>(variant);
range_length = range.size();
if (range_length == 2u) {
if (not is_none(range[0u])) {
m_tune_limits.first = get_value<int>(range[0u]);
}
if (not is_none(range[1u])) {
m_tune_limits.second = get_value<int>(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");
}
});
}
Expand Down
41 changes: 24 additions & 17 deletions src/script_interface/magnetostatics/DipolarP3M.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,28 +111,35 @@ class DipolarP3M : public Actor<DipolarP3M<Architecture>, ::DipolarP3M> {
m_tune_verbose = get_value<bool>(params, "verbose");
m_tune_limits = {std::nullopt, std::nullopt};
if (params.contains("tune_limits")) {
std::vector<Variant> range;
try {
auto const val = get_value<std::vector<int>>(params, "tune_limits");
assert(val.size() == 2u);
range.emplace_back(val[0u]);
range.emplace_back(val[1u]);
} catch (...) {
range = get_value<std::vector<Variant>>(params, "tune_limits");
assert(range.size() == 2u);
}
if (not is_none(range[0u])) {
m_tune_limits.first = get_value<int>(range[0u]);
}
if (not is_none(range[1u])) {
m_tune_limits.second = get_value<int>(range[1u]);
auto const &variant = params.at("tune_limits");
std::size_t range_length = 0u;
if (is_type<std::vector<int>>(variant)) {
auto const range = get_value<std::vector<int>>(variant);
range_length = range.size();
if (range_length == 2u) {
m_tune_limits = {range[0u], range[1u]};
}
} else {
auto const range = get_value<std::vector<Variant>>(variant);
range_length = range.size();
if (range_length == 2u) {
if (not is_none(range[0u])) {
m_tune_limits.first = get_value<int>(range[0u]);
}
if (not is_none(range[1u])) {
m_tune_limits.second = get_value<int>(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");
}
});
}
Expand Down
5 changes: 3 additions & 2 deletions testsuite/python/p3m_tuning_exceptions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"),
]
Expand Down

0 comments on commit 75f0812

Please sign in to comment.