Fix truncation error estimate for functions with symmetry

Improves p-AMR considerably.

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<double>::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<double>(last_index) + 0.5));
+    weight_value = exp(-square(0.5 - static_cast<double>(last_index - index)));
     // Add weighted power monitor
     weighted_average += weight_value * log10(mode);
     // Add term to weighted sum

src/NumericalAlgorithms/LinearOperators/PowerMonitors.hpp

@@ -70,7 +70,12 @@ std::array<DataVector, Dim> 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()`).

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<size_t, 2> coeff = {x_mode, y_mode};
+  const size_t x_mode = 0;
+  const size_t y_mode = 1;
+  const std::array<size_t, 2> 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<DataVector, 2> expected_power_monitors{
-      check_data_vector_x, check_data_vector_y};
+  const std::array<DataVector, 2> expected_power_monitors{check_data_vector_x,
+                                                          check_data_vector_y};
 
   CHECK_ITERABLE_APPROX(test_power_monitors, expected_power_monitors);
 }
@@ -156,11 +156,81 @@ 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<double>::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<double>::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",
                   "[NumericalAlgorithms][LinearOperators][Unit]") {
   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();
 }

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<Dim> mesh{4, Spectral::Basis::Legendre,
-                       Spectral::Quadrature::GaussLobatto};
-  const auto logical_coords = logical_coordinates(mesh);
-  // Manufacture some test data
-  tnsr::I<DataVector, Dim> 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<Dim> mesh{num_points, Spectral::Basis::Legendre,
+                         Spectral::Quadrature::GaussLobatto};
+    const auto logical_coords = logical_coordinates(mesh);
+    // Manufacture some test data
+    tnsr::I<DataVector, Dim> 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<Metavariables<Dim>> empty_cache{};
-  auto databox =
-      db::create<tmpl::list<::domain::Tags::Mesh<Dim>, TestVector<Dim>>>(
-          mesh, std::move(test_data));
-  ObservationBox<
-      tmpl::list<>,
-      db::DataBox<tmpl::list<::domain::Tags::Mesh<Dim>, TestVector<Dim>>>>
-      box{make_not_null(&databox)};
+    Parallel::GlobalCache<Metavariables<Dim>> empty_cache{};
+    auto databox =
+        db::create<tmpl::list<::domain::Tags::Mesh<Dim>, TestVector<Dim>>>(
+            mesh, std::move(test_data));
+    ObservationBox<
+        tmpl::list<>,
+        db::DataBox<tmpl::list<::domain::Tags::Mesh<Dim>, TestVector<Dim>>>>
+        box{make_not_null(&databox)};
 
-  const auto flags = criterion.evaluate(box, empty_cache, ElementId<Dim>{0});
-  CHECK(flags[0] == amr::Flag::IncreaseResolution);
-  CHECK(flags[1] == amr::Flag::DecreaseResolution);
+    return criterion.evaluate(box, empty_cache, ElementId<Dim>{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