intrp::callbacks::FindApparentHorizon< InterpolationTargetTag, Frame > Struct Template Reference

post interpolation callback (see intrp::protocols::PostInterpolationCallback) that does a FastFlow iteration and triggers another one until convergence. More...

#include <FindApparentHorizon.hpp>

Public Types
using	const_global_cache_tags = Parallel::get_const_global_cache_tags_from_actions< typename InterpolationTargetTag::post_horizon_find_callbacks >

Static Public Member Functions
template<typename DbTags , typename Metavariables , typename TemporalId >
static bool	apply (const gsl::not_null< db::DataBox< DbTags > * > box, const gsl::not_null< Parallel::GlobalCache< Metavariables > * > cache, const TemporalId &temporal_id)

Detailed Description

template<typename InterpolationTargetTag, typename Frame>
struct intrp::callbacks::FindApparentHorizon< InterpolationTargetTag, Frame >

post interpolation callback (see intrp::protocols::PostInterpolationCallback) that does a FastFlow iteration and triggers another one until convergence.

Assumes that InterpolationTargetTag contains an additional type alias called post_horizon_find_callbacks, which is a list of structs, each of which has a function

    if constexpr (MakeHorizonFinderFailOnPurpose) {
      // Make the analytic solution off-center on purpose, so that
      // the domain only partially contains the horizon and therefore
      // the interpolation fails.
      return {0.5, 0.0, 0.0};
    }
    return {0.0, 0.0, 0.0};
  }();
 
  // Create volume data and send it to the interpolator, for each temporal_id.
  for (const auto& temporal_id : temporal_ids) {
    for (const auto& element_id : element_ids) {
      const auto& block = domain.blocks()[element_id.block_id()];
      ::Mesh<3> mesh{domain_creator->initial_extents()[element_id.block_id()],
                     Spectral::Basis::Legendre,
                     Spectral::Quadrature::GaussLobatto};
 
      // If the map is time-independent, we always compute
      // analytic_solution_coords in the inertial frame.
      tnsr::I<
          DataVector, 3,
          std::conditional_t<IsTimeDependent::value, Frame, ::Frame::Inertial>>
          analytic_solution_coords{};
      if constexpr (std::is_same_v<Frame, ::Frame::Grid> and
                    IsTimeDependent::value) {
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        analytic_solution_coords =
            map_logical_to_grid(logical_coordinates(mesh));
      } else if constexpr (IsTimeDependent::value) {
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        // We don't have an Element ParallelComponent in this test, so
        // get the cache from the target component.
        const auto& cache =
            ActionTesting::cache<target_component>(runner, 0_st);
        const auto& functions_of_time =
            get<domain::Tags::FunctionsOfTime>(cache);
        if constexpr (std::is_same_v<Frame, ::Frame::Distorted>) {
          analytic_solution_coords = block.moving_mesh_grid_to_distorted_map()(
              map_logical_to_grid(logical_coordinates(mesh)), temporal_id,
              functions_of_time);
        } else {
          static_assert(std::is_same_v<Frame, ::Frame::Inertial>);
          analytic_solution_coords = block.moving_mesh_grid_to_inertial_map()(
              map_logical_to_grid(logical_coordinates(mesh)), temporal_id,
              functions_of_time);
        }
      } else {
        // Time-independent
        ElementMap<3, ::Frame::Inertial> map{
            element_id, block.stationary_map().get_clone()};
        analytic_solution_coords = map(logical_coordinates(mesh));
      }
 
      // Compute psi, pi, phi for KerrSchild.
      // Horizon is always at 0,0,0 in analytic_solution_coordinates.
      // Note that we always use Inertial frame if the map is time-independent.
      gr::Solutions::KerrSchild solution(mass, dimensionless_spin,
                                         analytic_solution_center);
      const auto solution_vars = solution.variables(
          analytic_solution_coords, 0.0,
          typename gr::Solutions::KerrSchild::tags<
              DataVector, std::conditional_t<IsTimeDependent::value, Frame,
                                             ::Frame::Inertial>>{});
 
      // Fill output variables with solution.
      typename ::Tags::Variables<typename metavars::interpolator_source_vars>::
          type output_vars(mesh.number_of_grid_points());
 
      if constexpr (std::is_same_v<Frame, ::Frame::Inertial> or
                    not IsTimeDependent::value) {
        // Easy case: Grid and Inertial frame are the same
 
        const auto& lapse = get<gr::Tags::Lapse<DataVector>>(solution_vars);
        const auto& dt_lapse =
            get<Tags::dt<gr::Tags::Lapse<DataVector>>>(solution_vars);
        const auto& d_lapse =
            get<typename gr::Solutions::KerrSchild ::DerivLapse<
                DataVector, ::Frame::Inertial>>(solution_vars);
        const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(solution_vars);
        const auto& d_shift =
            get<typename gr::Solutions::KerrSchild ::DerivShift<
                DataVector, ::Frame::Inertial>>(solution_vars);
        const auto& dt_shift =
            get<Tags::dt<gr::Tags::Shift<DataVector, 3>>>(solution_vars);
        const auto& g =
            get<gr::Tags::SpatialMetric<DataVector, 3>>(solution_vars);
        const auto& dt_g =
            get<Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>>(
                solution_vars);
        const auto& d_g =
            get<typename gr::Solutions::KerrSchild ::DerivSpatialMetric<
                DataVector, ::Frame::Inertial>>(solution_vars);
 
        get<::gr::Tags::SpacetimeMetric<DataVector, 3>>(output_vars) =
            gr::spacetime_metric(lapse, shift, g);
        get<::gh::Tags::Phi<DataVector, 3>>(output_vars) =
            gh::phi(lapse, d_lapse, shift, d_shift, g, d_g);
        get<::gh::Tags::Pi<DataVector, 3>>(output_vars) =
            gh::pi(lapse, dt_lapse, shift, dt_shift, g, dt_g,
                   get<::gh::Tags::Phi<DataVector, 3>>(output_vars));
      } else {
        // Frame is not Inertial, and we are time-dependent,
        // so need to transform tensors to
        // Inertial frame, since InterpolatorReceiveVolumeData always gets
        // its volume data in the Inertial frame.
 
        // The difficult parts are Pi and Phi, which are not tensors,
        // so they are not so easy to transform because we do not have
        // Hessians.
        //
        // So what we do is compute only the 3-metric, lapse, shift,
        // and extrinsic curvature, transform them (because they are
        // 3-tensors), and compute the other components by numerical
        // differentiation.
        const auto& lapse = get<gr::Tags::Lapse<DataVector>>(solution_vars);
        const auto& shift =
            get<gr::Tags::Shift<DataVector, 3, Frame>>(solution_vars);
        const auto& g =
            get<gr::Tags::SpatialMetric<DataVector, 3, Frame>>(solution_vars);
        const auto& K = get<gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>>(
            solution_vars);
 
        auto& cache = ActionTesting::cache<target_component>(runner, 0_st);
        const auto& functions_of_time =
            get<domain::Tags::FunctionsOfTime>(cache);
        const auto coords_frame_velocity_jacobians = [&block,
                                                      &analytic_solution_coords,
                                                      &temporal_id,
                                                      &functions_of_time]() {
          if constexpr (std::is_same_v<Frame, ::Frame::Grid>) {
            return block.moving_mesh_grid_to_inertial_map()
                .coords_frame_velocity_jacobians(
                    analytic_solution_coords, temporal_id, functions_of_time);
          } else {
            static_assert(std::is_same_v<Frame, ::Frame::Distorted>);
            return block.moving_mesh_distorted_to_inertial_map()
                .coords_frame_velocity_jacobians(
                    analytic_solution_coords, temporal_id, functions_of_time);
          }
        }();
        const auto& inv_jacobian = std::get<1>(coords_frame_velocity_jacobians);
        const auto& jacobian = std::get<2>(coords_frame_velocity_jacobians);
        const auto& frame_velocity =
            std::get<3>(coords_frame_velocity_jacobians);
 
        using inertial_metric_vars_tags =
            tmpl::list<gr::Tags::Lapse<DataVector>,
                       gr::Tags::Shift<DataVector, 3>,
                       gr::Tags::SpatialMetric<DataVector, 3>>;
        Variables<inertial_metric_vars_tags> inertial_metric_vars(
            mesh.number_of_grid_points());
 
        auto& shift_inertial =
            get<::gr::Tags::Shift<DataVector, 3>>(inertial_metric_vars);
        auto& lower_metric_inertial =
            get<::gr::Tags::SpatialMetric<DataVector, 3>>(inertial_metric_vars);
        // Just copy lapse, since it doesn't transform. Need it for derivs.
        get<gr::Tags::Lapse<DataVector>>(inertial_metric_vars) = lapse;
 
        for (size_t k = 0; k < 3; ++k) {
          shift_inertial.get(k) = -frame_velocity.get(k);
          for (size_t j = 0; j < 3; ++j) {
            shift_inertial.get(k) += jacobian.get(k, j) * shift.get(j);
          }
        }
 
        auto transform =
            [&inv_jacobian](
                const gsl::not_null<tnsr::ii<DataVector, 3, ::Frame::Inertial>*>
                    u_inertial,
                const tnsr::ii<DataVector, 3, Frame>& u_frame) {
              for (size_t i = 0; i < 3; ++i) {
                for (size_t j = i; j < 3; ++j) {  // symmetry
                  u_inertial->get(i, j) = 0.0;
                  for (size_t k = 0; k < 3; ++k) {
                    for (size_t p = 0; p < 3; ++p) {
                      u_inertial->get(i, j) += inv_jacobian.get(k, i) *
                                               inv_jacobian.get(p, j) *
                                               u_frame.get(k, p);
                    }
                  }
                }
              }
            };
 
        transform(make_not_null(&lower_metric_inertial), g);
        tnsr::ii<DataVector, 3, ::Frame::Inertial> extrinsic_curvature_inertial(
            mesh.number_of_grid_points());
        transform(make_not_null(&extrinsic_curvature_inertial), K);
 
        // Take spatial derivatives
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        const auto grid_deriv_inertial_metric_vars =
            partial_derivatives<inertial_metric_vars_tags>(
                inertial_metric_vars, mesh,
                map_logical_to_grid.inv_jacobian(logical_coordinates(mesh)));
        tnsr::ijj<DataVector, 3, ::Frame::Inertial> d_g(
            mesh.number_of_grid_points());
        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = i; j < 3; ++j) {  // symmetry
            for (size_t k = 0; k < 3; ++k) {
              d_g.get(k, i, j) = 0.0;
              for (size_t l = 0; l < 3; ++l) {
                d_g.get(k, i, j) +=
                    inv_jacobian.get(l, k) *
                    get<Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>,
                                    tmpl::size_t<3>, ::Frame::Grid>>(
                        grid_deriv_inertial_metric_vars)
                        .get(l, i, j);
              }
            }
          }
        }
        tnsr::iJ<DataVector, 3, ::Frame::Inertial> d_shift(
            mesh.number_of_grid_points());
        for (size_t i = 0; i < 3; ++i) {
          for (size_t k = 0; k < 3; ++k) {
            d_shift.get(k, i) = 0.0;
            for (size_t l = 0; l < 3; ++l) {
              d_shift.get(k, i) +=
                  inv_jacobian.get(l, k) *
                  get<Tags::deriv<gr::Tags::Shift<DataVector, 3>,
                                  tmpl::size_t<3>, ::Frame::Grid>>(
                      grid_deriv_inertial_metric_vars)
                      .get(l, i);
            }
          }
        }
        tnsr::i<DataVector, 3, ::Frame::Inertial> d_lapse(
            mesh.number_of_grid_points());
        for (size_t k = 0; k < 3; ++k) {
          d_lapse.get(k) = 0.0;
          for (size_t l = 0; l < 3; ++l) {
            d_lapse.get(k) +=
                inv_jacobian.get(l, k) *
                get<Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
                                ::Frame::Grid>>(grid_deriv_inertial_metric_vars)
                    .get(l);
          }
        }
 
        get<::gr::Tags::SpacetimeMetric<DataVector, 3>>(output_vars) =
            gr::spacetime_metric(lapse, shift_inertial, lower_metric_inertial);
        get<::gh::Tags::Phi<DataVector, 3>>(output_vars) =
            gh::phi(lapse, d_lapse, shift_inertial, d_shift,
                    lower_metric_inertial, d_g);
        // Compute Pi from extrinsic curvature and Phi.  Fill in zero
        // for zero components of Pi, since they won't be used at all
        // (can't fill in NaNs, because they will still be interpolated).
        const auto spacetime_normal_vector =
            gr::spacetime_normal_vector(lapse, shift_inertial);
        auto& Pi = get<::gh::Tags::Pi<DataVector, 3>>(output_vars);
        const auto& Phi = get<::gh::Tags::Phi<DataVector, 3>>(output_vars);
        for (size_t i = 0; i < 3; ++i) {
          Pi.get(i + 1, 0) = 0.0;
          for (size_t j = i; j < 3; ++j) {  // symmetry
            Pi.get(i + 1, j + 1) = 2.0 * extrinsic_curvature_inertial.get(i, j);
            for (size_t c = 0; c < 4; ++c) {
              Pi.get(i + 1, j + 1) -=
                  spacetime_normal_vector.get(c) *
                  (Phi.get(i, j + 1, c) + Phi.get(j, i + 1, c));
            }
          }
        }
        Pi.get(0, 0) = 0.0;
      }
 
      // Call the InterpolatorReceiveVolumeData action on each element_id.
      ActionTesting::simple_action<
          interp_component,
          intrp::Actions::InterpolatorReceiveVolumeData<::Tags::Time>>(
          make_not_null(&runner), mock_core_for_each_element.at(element_id),
          temporal_id, element_id, mesh, output_vars);
    }
  }
 
  // Invoke remaining actions in random order.
  MAKE_GENERATOR(generator);
  auto array_indices_with_queued_simple_actions =
      ActionTesting::array_indices_with_queued_simple_actions<
          typename metavars::component_list>(make_not_null(&runner));
  while (ActionTesting::number_of_elements_with_queued_simple_actions<
             typename metavars::component_list>(
             array_indices_with_queued_simple_actions) > 0) {
    ActionTesting::invoke_random_queued_simple_action<
        typename metavars::component_list>(
        make_not_null(&runner), make_not_null(&generator),
        array_indices_with_queued_simple_actions);
    array_indices_with_queued_simple_actions =
        ActionTesting::array_indices_with_queued_simple_actions<
            typename metavars::component_list>(make_not_null(&runner));
  }
 
  // Make sure function was called three times per
  // post_horizon_find_callback.
  CHECK(*test_horizon_called == 3 * tmpl::size<PostHorizonFindCallbacks>{});
}
 
// This tests the entire AH finder including numerical interpolation.
// We increase the timeout because we want to test time-independent
// and time-dependent cases, and we want to test more than one frame.
// Already the resolution used for the tests is very low
// (lmax=3, num_pts_per_dim=3 to 7) and the error
// tolerance Approx::custom().epsilon() is pretty large (1e-2 and 1e-3).
// [[TimeOut, 60]]
SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Interpolator.ApparentHorizonFinder",
                  "[Unit]") {
  domain::creators::register_derived_with_charm();
  domain::creators::time_dependence::register_derived_with_charm();
  domain::FunctionsOfTime::register_derived_with_charm();
 
  // Time-independent tests.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  // Note we have 2 post_horizon_find_callbacks as the first template
  // argument in the line below, just to test that everything works
  // with 2 post_horizon_find_callbacks.
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Inertial>,
                                   TestSchwarzschildHorizon<Frame::Inertial>>,
                        std::false_type>(&test_schwarzschild_horizon_called, 3,
                                         3, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Inertial>>,
                        std::false_type>(&test_kerr_horizon_called, 3, 5, 1.1,
                                         {{0.12, 0.23, 0.45}});
 
  // Time-independent tests with different frame tags.
  test_schwarzschild_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Grid>>,
                        std::false_type, Frame::Grid>(
      &test_schwarzschild_horizon_called, 3, 3, 1.0, {{0.0, 0.0, 0.0}});
 
  // Time-dependent tests.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Inertial>>,
                        std::true_type>(&test_schwarzschild_horizon_called, 3,
                                        4, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Inertial>>,
                        std::true_type>(&test_kerr_horizon_called, 3, 5, 1.1,
                                        {{0.12, 0.23, 0.45}});
 
  // Time-dependent tests in grid frame.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Grid>>,
                        std::true_type, Frame::Grid>(
      &test_schwarzschild_horizon_called, 3, 6, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Grid>>,
                        std::true_type, Frame::Grid>(
      &test_kerr_horizon_called, 3, 7, 1.1, {{0.12, 0.23, 0.45}});
 
  // Time-dependent tests in distorted frame using both translation and shape
  // maps.
  tmpl::for_each<tmpl::list<std::true_type, std::false_type>>(

post_horizon_find_callback_example

that is called if the FastFlow iteration has converged. InterpolationTargetTag also is assumed to contain an additional struct called horizon_find_failure_callback, which has a function

    if constexpr (MakeHorizonFinderFailOnPurpose) {
      // Make the analytic solution off-center on purpose, so that
      // the domain only partially contains the horizon and therefore
      // the interpolation fails.
      return {0.5, 0.0, 0.0};
    }
    return {0.0, 0.0, 0.0};
  }();
 
  // Create volume data and send it to the interpolator, for each temporal_id.
  for (const auto& temporal_id : temporal_ids) {
    for (const auto& element_id : element_ids) {
      const auto& block = domain.blocks()[element_id.block_id()];
      ::Mesh<3> mesh{domain_creator->initial_extents()[element_id.block_id()],
                     Spectral::Basis::Legendre,
                     Spectral::Quadrature::GaussLobatto};
 
      // If the map is time-independent, we always compute
      // analytic_solution_coords in the inertial frame.
      tnsr::I<
          DataVector, 3,
          std::conditional_t<IsTimeDependent::value, Frame, ::Frame::Inertial>>
          analytic_solution_coords{};
      if constexpr (std::is_same_v<Frame, ::Frame::Grid> and
                    IsTimeDependent::value) {
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        analytic_solution_coords =
            map_logical_to_grid(logical_coordinates(mesh));
      } else if constexpr (IsTimeDependent::value) {
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        // We don't have an Element ParallelComponent in this test, so
        // get the cache from the target component.
        const auto& cache =
            ActionTesting::cache<target_component>(runner, 0_st);
        const auto& functions_of_time =
            get<domain::Tags::FunctionsOfTime>(cache);
        if constexpr (std::is_same_v<Frame, ::Frame::Distorted>) {
          analytic_solution_coords = block.moving_mesh_grid_to_distorted_map()(
              map_logical_to_grid(logical_coordinates(mesh)), temporal_id,
              functions_of_time);
        } else {
          static_assert(std::is_same_v<Frame, ::Frame::Inertial>);
          analytic_solution_coords = block.moving_mesh_grid_to_inertial_map()(
              map_logical_to_grid(logical_coordinates(mesh)), temporal_id,
              functions_of_time);
        }
      } else {
        // Time-independent
        ElementMap<3, ::Frame::Inertial> map{
            element_id, block.stationary_map().get_clone()};
        analytic_solution_coords = map(logical_coordinates(mesh));
      }
 
      // Compute psi, pi, phi for KerrSchild.
      // Horizon is always at 0,0,0 in analytic_solution_coordinates.
      // Note that we always use Inertial frame if the map is time-independent.
      gr::Solutions::KerrSchild solution(mass, dimensionless_spin,
                                         analytic_solution_center);
      const auto solution_vars = solution.variables(
          analytic_solution_coords, 0.0,
          typename gr::Solutions::KerrSchild::tags<
              DataVector, std::conditional_t<IsTimeDependent::value, Frame,
                                             ::Frame::Inertial>>{});
 
      // Fill output variables with solution.
      typename ::Tags::Variables<typename metavars::interpolator_source_vars>::
          type output_vars(mesh.number_of_grid_points());
 
      if constexpr (std::is_same_v<Frame, ::Frame::Inertial> or
                    not IsTimeDependent::value) {
        // Easy case: Grid and Inertial frame are the same
 
        const auto& lapse = get<gr::Tags::Lapse<DataVector>>(solution_vars);
        const auto& dt_lapse =
            get<Tags::dt<gr::Tags::Lapse<DataVector>>>(solution_vars);
        const auto& d_lapse =
            get<typename gr::Solutions::KerrSchild ::DerivLapse<
                DataVector, ::Frame::Inertial>>(solution_vars);
        const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(solution_vars);
        const auto& d_shift =
            get<typename gr::Solutions::KerrSchild ::DerivShift<
                DataVector, ::Frame::Inertial>>(solution_vars);
        const auto& dt_shift =
            get<Tags::dt<gr::Tags::Shift<DataVector, 3>>>(solution_vars);
        const auto& g =
            get<gr::Tags::SpatialMetric<DataVector, 3>>(solution_vars);
        const auto& dt_g =
            get<Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>>(
                solution_vars);
        const auto& d_g =
            get<typename gr::Solutions::KerrSchild ::DerivSpatialMetric<
                DataVector, ::Frame::Inertial>>(solution_vars);
 
        get<::gr::Tags::SpacetimeMetric<DataVector, 3>>(output_vars) =
            gr::spacetime_metric(lapse, shift, g);
        get<::gh::Tags::Phi<DataVector, 3>>(output_vars) =
            gh::phi(lapse, d_lapse, shift, d_shift, g, d_g);
        get<::gh::Tags::Pi<DataVector, 3>>(output_vars) =
            gh::pi(lapse, dt_lapse, shift, dt_shift, g, dt_g,
                   get<::gh::Tags::Phi<DataVector, 3>>(output_vars));
      } else {
        // Frame is not Inertial, and we are time-dependent,
        // so need to transform tensors to
        // Inertial frame, since InterpolatorReceiveVolumeData always gets
        // its volume data in the Inertial frame.
 
        // The difficult parts are Pi and Phi, which are not tensors,
        // so they are not so easy to transform because we do not have
        // Hessians.
        //
        // So what we do is compute only the 3-metric, lapse, shift,
        // and extrinsic curvature, transform them (because they are
        // 3-tensors), and compute the other components by numerical
        // differentiation.
        const auto& lapse = get<gr::Tags::Lapse<DataVector>>(solution_vars);
        const auto& shift =
            get<gr::Tags::Shift<DataVector, 3, Frame>>(solution_vars);
        const auto& g =
            get<gr::Tags::SpatialMetric<DataVector, 3, Frame>>(solution_vars);
        const auto& K = get<gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>>(
            solution_vars);
 
        auto& cache = ActionTesting::cache<target_component>(runner, 0_st);
        const auto& functions_of_time =
            get<domain::Tags::FunctionsOfTime>(cache);
        const auto coords_frame_velocity_jacobians = [&block,
                                                      &analytic_solution_coords,
                                                      &temporal_id,
                                                      &functions_of_time]() {
          if constexpr (std::is_same_v<Frame, ::Frame::Grid>) {
            return block.moving_mesh_grid_to_inertial_map()
                .coords_frame_velocity_jacobians(
                    analytic_solution_coords, temporal_id, functions_of_time);
          } else {
            static_assert(std::is_same_v<Frame, ::Frame::Distorted>);
            return block.moving_mesh_distorted_to_inertial_map()
                .coords_frame_velocity_jacobians(
                    analytic_solution_coords, temporal_id, functions_of_time);
          }
        }();
        const auto& inv_jacobian = std::get<1>(coords_frame_velocity_jacobians);
        const auto& jacobian = std::get<2>(coords_frame_velocity_jacobians);
        const auto& frame_velocity =
            std::get<3>(coords_frame_velocity_jacobians);
 
        using inertial_metric_vars_tags =
            tmpl::list<gr::Tags::Lapse<DataVector>,
                       gr::Tags::Shift<DataVector, 3>,
                       gr::Tags::SpatialMetric<DataVector, 3>>;
        Variables<inertial_metric_vars_tags> inertial_metric_vars(
            mesh.number_of_grid_points());
 
        auto& shift_inertial =
            get<::gr::Tags::Shift<DataVector, 3>>(inertial_metric_vars);
        auto& lower_metric_inertial =
            get<::gr::Tags::SpatialMetric<DataVector, 3>>(inertial_metric_vars);
        // Just copy lapse, since it doesn't transform. Need it for derivs.
        get<gr::Tags::Lapse<DataVector>>(inertial_metric_vars) = lapse;
 
        for (size_t k = 0; k < 3; ++k) {
          shift_inertial.get(k) = -frame_velocity.get(k);
          for (size_t j = 0; j < 3; ++j) {
            shift_inertial.get(k) += jacobian.get(k, j) * shift.get(j);
          }
        }
 
        auto transform =
            [&inv_jacobian](
                const gsl::not_null<tnsr::ii<DataVector, 3, ::Frame::Inertial>*>
                    u_inertial,
                const tnsr::ii<DataVector, 3, Frame>& u_frame) {
              for (size_t i = 0; i < 3; ++i) {
                for (size_t j = i; j < 3; ++j) {  // symmetry
                  u_inertial->get(i, j) = 0.0;
                  for (size_t k = 0; k < 3; ++k) {
                    for (size_t p = 0; p < 3; ++p) {
                      u_inertial->get(i, j) += inv_jacobian.get(k, i) *
                                               inv_jacobian.get(p, j) *
                                               u_frame.get(k, p);
                    }
                  }
                }
              }
            };
 
        transform(make_not_null(&lower_metric_inertial), g);
        tnsr::ii<DataVector, 3, ::Frame::Inertial> extrinsic_curvature_inertial(
            mesh.number_of_grid_points());
        transform(make_not_null(&extrinsic_curvature_inertial), K);
 
        // Take spatial derivatives
        ElementMap<3, ::Frame::Grid> map_logical_to_grid{
            element_id, block.moving_mesh_logical_to_grid_map().get_clone()};
        const auto grid_deriv_inertial_metric_vars =
            partial_derivatives<inertial_metric_vars_tags>(
                inertial_metric_vars, mesh,
                map_logical_to_grid.inv_jacobian(logical_coordinates(mesh)));
        tnsr::ijj<DataVector, 3, ::Frame::Inertial> d_g(
            mesh.number_of_grid_points());
        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = i; j < 3; ++j) {  // symmetry
            for (size_t k = 0; k < 3; ++k) {
              d_g.get(k, i, j) = 0.0;
              for (size_t l = 0; l < 3; ++l) {
                d_g.get(k, i, j) +=
                    inv_jacobian.get(l, k) *
                    get<Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>,
                                    tmpl::size_t<3>, ::Frame::Grid>>(
                        grid_deriv_inertial_metric_vars)
                        .get(l, i, j);
              }
            }
          }
        }
        tnsr::iJ<DataVector, 3, ::Frame::Inertial> d_shift(
            mesh.number_of_grid_points());
        for (size_t i = 0; i < 3; ++i) {
          for (size_t k = 0; k < 3; ++k) {
            d_shift.get(k, i) = 0.0;
            for (size_t l = 0; l < 3; ++l) {
              d_shift.get(k, i) +=
                  inv_jacobian.get(l, k) *
                  get<Tags::deriv<gr::Tags::Shift<DataVector, 3>,
                                  tmpl::size_t<3>, ::Frame::Grid>>(
                      grid_deriv_inertial_metric_vars)
                      .get(l, i);
            }
          }
        }
        tnsr::i<DataVector, 3, ::Frame::Inertial> d_lapse(
            mesh.number_of_grid_points());
        for (size_t k = 0; k < 3; ++k) {
          d_lapse.get(k) = 0.0;
          for (size_t l = 0; l < 3; ++l) {
            d_lapse.get(k) +=
                inv_jacobian.get(l, k) *
                get<Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
                                ::Frame::Grid>>(grid_deriv_inertial_metric_vars)
                    .get(l);
          }
        }
 
        get<::gr::Tags::SpacetimeMetric<DataVector, 3>>(output_vars) =
            gr::spacetime_metric(lapse, shift_inertial, lower_metric_inertial);
        get<::gh::Tags::Phi<DataVector, 3>>(output_vars) =
            gh::phi(lapse, d_lapse, shift_inertial, d_shift,
                    lower_metric_inertial, d_g);
        // Compute Pi from extrinsic curvature and Phi.  Fill in zero
        // for zero components of Pi, since they won't be used at all
        // (can't fill in NaNs, because they will still be interpolated).
        const auto spacetime_normal_vector =
            gr::spacetime_normal_vector(lapse, shift_inertial);
        auto& Pi = get<::gh::Tags::Pi<DataVector, 3>>(output_vars);
        const auto& Phi = get<::gh::Tags::Phi<DataVector, 3>>(output_vars);
        for (size_t i = 0; i < 3; ++i) {
          Pi.get(i + 1, 0) = 0.0;
          for (size_t j = i; j < 3; ++j) {  // symmetry
            Pi.get(i + 1, j + 1) = 2.0 * extrinsic_curvature_inertial.get(i, j);
            for (size_t c = 0; c < 4; ++c) {
              Pi.get(i + 1, j + 1) -=
                  spacetime_normal_vector.get(c) *
                  (Phi.get(i, j + 1, c) + Phi.get(j, i + 1, c));
            }
          }
        }
        Pi.get(0, 0) = 0.0;
      }
 
      // Call the InterpolatorReceiveVolumeData action on each element_id.
      ActionTesting::simple_action<
          interp_component,
          intrp::Actions::InterpolatorReceiveVolumeData<::Tags::Time>>(
          make_not_null(&runner), mock_core_for_each_element.at(element_id),
          temporal_id, element_id, mesh, output_vars);
    }
  }
 
  // Invoke remaining actions in random order.
  MAKE_GENERATOR(generator);
  auto array_indices_with_queued_simple_actions =
      ActionTesting::array_indices_with_queued_simple_actions<
          typename metavars::component_list>(make_not_null(&runner));
  while (ActionTesting::number_of_elements_with_queued_simple_actions<
             typename metavars::component_list>(
             array_indices_with_queued_simple_actions) > 0) {
    ActionTesting::invoke_random_queued_simple_action<
        typename metavars::component_list>(
        make_not_null(&runner), make_not_null(&generator),
        array_indices_with_queued_simple_actions);
    array_indices_with_queued_simple_actions =
        ActionTesting::array_indices_with_queued_simple_actions<
            typename metavars::component_list>(make_not_null(&runner));
  }
 
  // Make sure function was called three times per
  // post_horizon_find_callback.
  CHECK(*test_horizon_called == 3 * tmpl::size<PostHorizonFindCallbacks>{});
}
 
// This tests the entire AH finder including numerical interpolation.
// We increase the timeout because we want to test time-independent
// and time-dependent cases, and we want to test more than one frame.
// Already the resolution used for the tests is very low
// (lmax=3, num_pts_per_dim=3 to 7) and the error
// tolerance Approx::custom().epsilon() is pretty large (1e-2 and 1e-3).
// [[TimeOut, 60]]
SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Interpolator.ApparentHorizonFinder",
                  "[Unit]") {
  domain::creators::register_derived_with_charm();
  domain::creators::time_dependence::register_derived_with_charm();
  domain::FunctionsOfTime::register_derived_with_charm();
 
  // Time-independent tests.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  // Note we have 2 post_horizon_find_callbacks as the first template
  // argument in the line below, just to test that everything works
  // with 2 post_horizon_find_callbacks.
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Inertial>,
                                   TestSchwarzschildHorizon<Frame::Inertial>>,
                        std::false_type>(&test_schwarzschild_horizon_called, 3,
                                         3, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Inertial>>,
                        std::false_type>(&test_kerr_horizon_called, 3, 5, 1.1,
                                         {{0.12, 0.23, 0.45}});
 
  // Time-independent tests with different frame tags.
  test_schwarzschild_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Grid>>,
                        std::false_type, Frame::Grid>(
      &test_schwarzschild_horizon_called, 3, 3, 1.0, {{0.0, 0.0, 0.0}});
 
  // Time-dependent tests.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Inertial>>,
                        std::true_type>(&test_schwarzschild_horizon_called, 3,
                                        4, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Inertial>>,
                        std::true_type>(&test_kerr_horizon_called, 3, 5, 1.1,
                                        {{0.12, 0.23, 0.45}});
 
  // Time-dependent tests in grid frame.
  test_schwarzschild_horizon_called = 0;
  test_kerr_horizon_called = 0;
  test_apparent_horizon<tmpl::list<TestSchwarzschildHorizon<Frame::Grid>>,
                        std::true_type, Frame::Grid>(
      &test_schwarzschild_horizon_called, 3, 6, 1.0, {{0.0, 0.0, 0.0}});
  test_apparent_horizon<tmpl::list<TestKerrHorizon<Frame::Grid>>,
                        std::true_type, Frame::Grid>(
      &test_kerr_horizon_called, 3, 7, 1.1, {{0.12, 0.23, 0.45}});
 
  // Time-dependent tests in distorted frame using both translation and shape
  // maps.
  tmpl::for_each<tmpl::list<std::true_type, std::false_type>>(

horizon_find_failure_callback_example

that is called if the FastFlow iteration or the interpolation has failed.

Uses:

• Metavariables:
  • temporal_id
• DataBox:
  • logging::Tags::Verbosity<InterpolationTargetTag>
  • gr::Tags::InverseSpatialMetric<DataVector, 3, Frame>
  • gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame>
  • gr::Tags::SpatialChristoffelSecondKind<DataVector, 3, Frame>
  • ah::Tags::FastFlow
  • ylm::Tags::Strahlkorper<Frame>

Modifies:

• DataBox:
  • ah::Tags::FastFlow
  • ylm::Tags::Strahlkorper<Frame>

This can be used in InterpolationTargetTag::post_interpolation_callbacks; see intrp::protocols::InterpolationTargetTag for details on InterpolationTargetTag.

Output

Optionally, a single line of output is printed to stdout either on each iteration (if verbosity > Verbosity::Quiet) or on convergence (if verbosity > Verbosity::Silent). The output consists of the following labels, and values associated with each label:

• t = time given by value of temporal_id argument to apply
• its = current iteration number.
• R = min and max of residual over all prolonged grid points.
• |R| = L2 norm of residual, counting only L modes solved for.
• |R_mesh| = L2 norm of residual over prolonged grid points.
• r = min and max radius of trial horizon surface.

Difference between |R| and |R_mesh|:

The horizon is represented in a $Y_{lm}$ expansion up to $l=l_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}$; the residual |R| represents the failure of that surface to satisfy the apparent horizon equation.

However, at each iteration we also interpolate the horizon surface to a higher resolution ("prolongation"). The prolonged surface includes $Y_{lm}$ coefficients up to $l=l_{\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}}$, where $l_{\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{h}}>l_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}$. The residual computed on this higher-resolution surface is |R_mesh|.

As iterations proceed, |R| should decrease until it reaches numerical roundoff error, because we are varying all the $Y_{lm}$ coefficients up to $l=l_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}$ to minimize the residual. However, |R_mesh| should eventually stop decreasing as iterations proceed, hitting a floor that represents the numerical truncation error of the solution.

The convergence criterion looks at both |R| and |R_mesh|: Once |R| is small enough and |R_mesh| has stabilized, it is pointless to keep iterating (even though one could iterate until |R| reaches roundoff).

Problems with convergence of the apparent horizon finder can often be diagnosed by looking at the behavior of |R| and |R_mesh| as a function of iteration.

---

The documentation for this struct was generated from the following file:

• src/ParallelAlgorithms/ApparentHorizonFinder/Callbacks/FindApparentHorizon.hpp