diff --git a/src/NumericalAlgorithms/LinearOperators/PowerMonitors.cpp b/src/NumericalAlgorithms/LinearOperators/PowerMonitors.cpp index bca5c3683ec2e..b66bdfdb96f8b 100644 --- a/src/NumericalAlgorithms/LinearOperators/PowerMonitors.cpp +++ b/src/NumericalAlgorithms/LinearOperators/PowerMonitors.cpp @@ -63,19 +63,31 @@ double relative_truncation_error(const DataVector& power_monitor, "Number of modes needs less or equal than the number of power monitors"); ASSERT(2_st <= num_modes_to_use, "Number of modes needs to be larger or equal than 2."); + const size_t last_index = num_modes_to_use - 1; + const double max_mode = blaze::max(power_monitor); + const double cutoff = + 100. * std::numeric_limits::epsilon() * max_mode; + // If the last two or more modes are zero, assume that the function is + // represented exactly and return a relative truncation error of zero. + // Just one zero mode is not enough to make this assumption, as the function + // could have zero modes by symmetry. + if (num_modes_to_use >= 2 and power_monitor[last_index] < cutoff and + power_monitor[last_index - 1] < cutoff) { + return -log10(cutoff) + 2.; + } // Compute weighted average and total sum in the current dimension double weighted_average = 0.0; double weight_sum = 0.0; double weight_value = 0.0; - const size_t last_index = num_modes_to_use - 1; for (size_t index = 0; index <= last_index; ++index) { const double mode = power_monitor[index]; - if (mode == 0.) { + if (mode < cutoff) { + // Ignore modes below this cutoff, so modes that are zero (e.g. by + // symmetry) don't make us underestimate the truncation error. continue; } // Compute current weight - weight_value = - exp(-square(index - static_cast(last_index) + 0.5)); + weight_value = exp(-square(0.5 - static_cast(last_index - index))); // Add weighted power monitor weighted_average += weight_value * log10(mode); // Add term to weighted sum diff --git a/src/NumericalAlgorithms/LinearOperators/PowerMonitors.hpp b/src/NumericalAlgorithms/LinearOperators/PowerMonitors.hpp index 8bb84e58fc690..fb4e19dd62ab4 100644 --- a/src/NumericalAlgorithms/LinearOperators/PowerMonitors.hpp +++ b/src/NumericalAlgorithms/LinearOperators/PowerMonitors.hpp @@ -70,7 +70,12 @@ std::array power_monitors(const DataVector& u, * average with larger weights toward the highest modes. This number should * correspond to the number of digits resolved by the spectral expansion. * - * \note Modes that are identically zero are ignored in the weighted average. + * \note Modes below a cutoff of $100 \epsilon \mathrm{max}_k(P_k)$ are ignored + * in the weighted average, where $\epsilon$ is the machine epsilon. This + * ensures that we don't underestimate the truncation error if some modes are + * zero (e.g. by symmetry). Furthermore, if the last two or more modes are zero, + * we assume that the function is represented exactly and return a relative + * truncation error of zero. * * \details The number of modes (`num_modes_to_use`) argument needs to be less * or equal than the total number of power monitors (`power_monitor.size()`). diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PowerMonitors.cpp b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PowerMonitors.cpp index 4668297c542e9..dded77888e648 100644 --- a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PowerMonitors.cpp +++ b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_PowerMonitors.cpp @@ -24,8 +24,8 @@ namespace { void test_power_monitors_impl() { - size_t number_of_points_per_dimension = 4; - size_t number_of_points = pow<2>(number_of_points_per_dimension); + const size_t number_of_points_per_dimension = 4; + const size_t number_of_points = pow<2>(number_of_points_per_dimension); // Test a constant function const DataVector test_data_vector = DataVector{number_of_points, 1.0}; @@ -34,7 +34,7 @@ void test_power_monitors_impl() { Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto}; - auto test_power_monitors = + const auto test_power_monitors = PowerMonitors::power_monitors<2_st>(test_data_vector, mesh); // The only non-zero modal coefficient of a constant is the one corresponding @@ -50,7 +50,7 @@ void test_power_monitors_impl() { } void test_power_monitors_second_impl() { - size_t number_of_points_per_dimension = 4; + const size_t number_of_points_per_dimension = 4; const Mesh<2_st> mesh{number_of_points_per_dimension, Spectral::Basis::Legendre, @@ -60,9 +60,9 @@ void test_power_monitors_second_impl() { // Build a test function containing only one Legendre basis function // per dimension - size_t x_mode = 0; - size_t y_mode = 1; - std::array coeff = {x_mode, y_mode}; + const size_t x_mode = 0; + const size_t y_mode = 1; + const std::array coeff = {x_mode, y_mode}; DataVector u_nodal(mesh.number_of_grid_points(), 1.0); for (size_t dim = 0; dim < 2; ++dim) { @@ -71,7 +71,7 @@ void test_power_monitors_second_impl() { gsl::at(coeff, dim), logical_coords.get(dim)); } - auto test_power_monitors = + const auto test_power_monitors = PowerMonitors::power_monitors<2_st>(u_nodal, mesh); // The only non-zero modal coefficient of a constant is the one corresponding @@ -88,8 +88,8 @@ void test_power_monitors_second_impl() { check_data_vector_y[y_mode] = 1.0 / sqrt(number_of_points_per_dimension); // We compare against the expected array - const std::array expected_power_monitors{ - check_data_vector_x, check_data_vector_y}; + const std::array expected_power_monitors{check_data_vector_x, + check_data_vector_y}; CHECK_ITERABLE_APPROX(test_power_monitors, expected_power_monitors); } @@ -156,6 +156,74 @@ void test_relative_truncation_error_impl() { CHECK_ITERABLE_APPROX(test_truncation_error, expected_truncation_error_x); } +void test_relative_truncation_error_with_symmetry() { + // Try to resolve half a period of a sinusoid + const size_t num_modes = 12; + const Mesh<1> mesh{num_modes, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const auto xi = Spectral::collocation_points(mesh); + const double wave_number = 0.5; + const DataVector u_nodal = sin((xi + 1.) * M_PI * wave_number); + CAPTURE(u_nodal); + auto modes = PowerMonitors::power_monitors(u_nodal, mesh)[0]; + // Add some more noise to the modes + modes += 10. * std::numeric_limits::epsilon(); + CAPTURE(modes); + const double relative_truncation_error = + PowerMonitors::relative_truncation_error(modes, num_modes); + // Last mode should be zero by symmetry + REQUIRE(modes[num_modes - 1] == approx(0.)); + // Expect the relative truncation error to be the ratio of the first and last + // nonzero modes + const double expected_relative_truncation_error = + log10(modes[0] / modes[num_modes - 2]); + Approx custom_approx = Approx::custom().epsilon(1e-2).scale(1.); + CHECK(relative_truncation_error == + custom_approx(expected_relative_truncation_error)); +} + +void test_relative_truncation_error_linear_function() { + // Resolve a linear function with a few modes. We technically need only 2. + const auto get_modes = [](const size_t num_modes) { + const Mesh<1> mesh{num_modes, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const auto xi = Spectral::collocation_points(mesh); + const DataVector u_nodal = (xi + 1.) * 0.5; + auto modes = PowerMonitors::power_monitors(u_nodal, mesh)[0]; + // Add some noise to the modes + modes += 10. * std::numeric_limits::epsilon(); + const double relative_truncation_error = + PowerMonitors::relative_truncation_error(modes, num_modes); + return std::make_pair(modes, relative_truncation_error); + }; + { + INFO("2 modes"); + const auto [modes, rel_error_digits] = get_modes(2); + CAPTURE(modes); + CHECK_ITERABLE_APPROX(modes, (DataVector{0.5, 0.5})); + // We don't know for sure that we have resolved the function exactly, + // because we have two nonzero modes and nothing else. + CHECK(rel_error_digits == approx(0.)); + } + { + INFO("3 modes"); + const auto [modes, rel_error_digits] = get_modes(3); + CAPTURE(modes); + CHECK_ITERABLE_APPROX(modes, (DataVector{0.5, 0.5, 0.})); + // The last mode is zero, but we still don't know if we have resolved the + // function because the last mode could be zero by symmetry. + CHECK(rel_error_digits == approx(0.)); + } + { + INFO("4 modes"); + const auto [modes, rel_error_digits] = get_modes(4); + CAPTURE(modes); + CHECK_ITERABLE_APPROX(modes, (DataVector{0.5, 0.5, 0., 0.})); + // We have two zero modes, so we know we have resolved the function exactly. + CHECK(rel_error_digits > 14.); + } +} + } // namespace SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.PowerMonitors", @@ -163,4 +231,6 @@ SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.PowerMonitors", test_power_monitors_impl(); test_power_monitors_second_impl(); test_relative_truncation_error_impl(); + test_relative_truncation_error_with_symmetry(); + test_relative_truncation_error_linear_function(); } diff --git a/tests/Unit/ParallelAlgorithms/Amr/Criteria/Test_TruncationError.cpp b/tests/Unit/ParallelAlgorithms/Amr/Criteria/Test_TruncationError.cpp index f27be4e383e34..2c2498dc59e70 100644 --- a/tests/Unit/ParallelAlgorithms/Amr/Criteria/Test_TruncationError.cpp +++ b/tests/Unit/ParallelAlgorithms/Amr/Criteria/Test_TruncationError.cpp @@ -59,29 +59,57 @@ SPECTRE_TEST_CASE("Unit.Amr.Criteria.TruncationError", " AbsoluteTarget: 1.e-3\n" " RelativeTarget: 1.e-3\n"); - const Mesh mesh{4, Spectral::Basis::Legendre, - Spectral::Quadrature::GaussLobatto}; - const auto logical_coords = logical_coordinates(mesh); - // Manufacture some test data - tnsr::I test_data{}; - // X-component is linear in x and y, so it is exactly represented on the mesh - get<0>(test_data) = get<0>(logical_coords) + 2. * get<1>(logical_coords); - // Y-component is nonlinear in one dimension and linear in the other - get<1>(test_data) = - exp(sin(M_PI * get<0>(logical_coords))) + 2. * get<1>(logical_coords); + const auto evaluate_criterion = [&criterion](const size_t num_points) { + const Mesh mesh{num_points, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const auto logical_coords = logical_coordinates(mesh); + // Manufacture some test data + tnsr::I test_data{}; + // X-component is linear in x and y, so it is exactly represented on the + // mesh + get<0>(test_data) = get<0>(logical_coords) + 2. * get<1>(logical_coords); + // Y-component is nonlinear in one dimension and linear in the other + get<1>(test_data) = + exp(sin(M_PI * get<0>(logical_coords))) + 2. * get<1>(logical_coords); - Parallel::GlobalCache> empty_cache{}; - auto databox = - db::create, TestVector>>( - mesh, std::move(test_data)); - ObservationBox< - tmpl::list<>, - db::DataBox, TestVector>>> - box{make_not_null(&databox)}; + Parallel::GlobalCache> empty_cache{}; + auto databox = + db::create, TestVector>>( + mesh, std::move(test_data)); + ObservationBox< + tmpl::list<>, + db::DataBox, TestVector>>> + box{make_not_null(&databox)}; - const auto flags = criterion.evaluate(box, empty_cache, ElementId{0}); - CHECK(flags[0] == amr::Flag::IncreaseResolution); - CHECK(flags[1] == amr::Flag::DecreaseResolution); + return criterion.evaluate(box, empty_cache, ElementId{0}); + }; + + // Expectation: + // - In the first dimension, one of the components is nonlinear, so we need + // lots of resolution there. + // - In the second dimension, both components are linear. Technically, we need + // only 2 modes to represent the functions exactly. However, if we have only + // 2 modes we can't be sure numerically that we're resolving the function. + // Also with 3 modes we can't be sure, because the third mode might be zero + // by symmetry. So we need 4 modes in this dimension. + { + INFO("3 modes"); + const auto flags = evaluate_criterion(3); + CHECK(flags[0] == amr::Flag::IncreaseResolution); + CHECK(flags[1] == amr::Flag::IncreaseResolution); + } + { + INFO("4 modes"); + const auto flags = evaluate_criterion(4); + CHECK(flags[0] == amr::Flag::IncreaseResolution); + CHECK(flags[1] == amr::Flag::DoNothing); + } + { + INFO("5 modes"); + const auto flags = evaluate_criterion(5); + CHECK(flags[0] == amr::Flag::IncreaseResolution); + CHECK(flags[1] == amr::Flag::DecreaseResolution); + } } } // namespace amr::Criteria