SpECTRE  v2024.04.12
Public Types | Static Public Member Functions | List of all members
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, 40]]
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>>(
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:48
ElementMap
The CoordinateMap for the Element from the Logical frame to the TargetFrame
Definition: ElementMap.hpp:35
Tensor
Definition: ContractFirstNIndices.hpp:16
gr::Solutions::KerrSchild
Kerr black hole in Kerr-Schild coordinates.
Definition: KerrSchild.hpp:219
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
logical_coordinates
void logical_coordinates(gsl::not_null< tnsr::I< DataVector, VolumeDim, Frame::ElementLogical > * > logical_coords, const Mesh< VolumeDim > &mesh)
Compute the logical coordinates in an Element.
db::get
const auto & get(const DataBox< TagList > &box)
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1264
gh::phi
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > * > phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric)
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein's equations...
gr::spacetime_metric
void spacetime_metric(gsl::not_null< tnsr::aa< DataType, Dim, Frame > * > spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, Dim, Frame > &shift, const tnsr::ii< DataType, Dim, Frame > &spatial_metric)
Computes the spacetime metric from the spatial metric, lapse, and shift.
gh::pi
void pi(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > * > pi, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi)
Computes the conjugate momentum of the spacetime metric .
gr::shift
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric)
Compute shift from spacetime metric and inverse spatial metric.
gr::spacetime_normal_vector
tnsr::A< DataType, SpatialDim, Frame > spacetime_normal_vector(const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift)
Computes spacetime normal vector from lapse and shift.
gr::lapse
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric)
Compute lapse from shift and spacetime metric.
MAKE_GENERATOR
#define MAKE_GENERATOR(...)
MAKE_GENERATOR(NAME [, SEED]) declares a variable of name NAME containing a generator of type std::mt...
Definition: TestHelpers.hpp:419
TimeSteppers::transform
void transform(const gsl::not_null< History< DestVars > * > dest, const History< SourceVars > &source, ValueTransformer &&value_transformer, DerivativeTransformer &&derivative_transformer)
Initialize a History object based on the contents of another, applying a transformation to each value...
Definition: History.hpp:1046
value
constexpr T & value(T &t)
Returns t.value() if t is a std::optional otherwise returns t.
Definition: OptionalHelpers.hpp:32
ActionTesting::array_indices_with_queued_simple_actions
auto array_indices_with_queued_simple_actions(const gsl::not_null< MockRuntimeSystem * > runner) -> detail::array_indices_for_each_component< ComponentList >
Returns a vector of array_indices for each Component. The vector is filled with only those array_indi...
Definition: MockRuntimeSystemFreeFunctions.hpp:427
ActionTesting::number_of_elements_with_queued_simple_actions
size_t number_of_elements_with_queued_simple_actions(const detail::array_indices_for_each_component< ComponentList > &array_indices)
Given the output of array_indices_with_queued_simple_actions, returns the total number of array indic...
Definition: MockRuntimeSystemFreeFunctions.hpp:448
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index)
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:379
ActionTesting::invoke_random_queued_simple_action
void invoke_random_queued_simple_action(const gsl::not_null< MockRuntimeSystem * > runner, const gsl::not_null< Generator * > generator, const detail::array_indices_for_each_component< ComponentList > &array_indices)
Invokes the next queued action on a random Component on a random array_index of that component....
Definition: MockRuntimeSystemFreeFunctions.hpp:465
ActionTesting::simple_action
void simple_action(const gsl::not_null< MockRuntimeSystem< Metavariables > * > runner, const typename Component::array_index &array_index, Args &&... args)
Runs the simple action Action on the array_indexth element of the parallel component Component.
Definition: MockRuntimeSystemFreeFunctions.hpp:286
Frame
Indicates the Frame that a TensorIndexType is in.
Definition: IndexType.hpp:36
domain
Holds entities related to the computational domain.
Definition: BoundaryCondition.hpp:14
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh)
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
transform
Holds functions related to transforming between frames.
Definition: FrameTransform.hpp:18
ylm::jacobian
ylm::Tags::aliases::Jacobian< Fr > jacobian(const tnsr::i< DataVector, 2, ::Frame::Spherical< Fr > > &theta_phi)
ylm::inv_jacobian
ylm::Tags::aliases::InvJacobian< Fr > inv_jacobian(const tnsr::i< DataVector, 2, ::Frame::Spherical< Fr > > &theta_phi)
Frame::Grid
Definition: IndexType.hpp:45
Frame::Inertial
Definition: IndexType.hpp:46
gh::Tags::Phi
Auxiliary variable which is analytically the spatial derivative of the spacetime metric.
Definition: Tags.hpp:38
gr::Tags::Shift
Definition: Tags.hpp:62
gr::Tags::SpatialMetric
Definition: Tags.hpp:26
intrp::Actions::InterpolatorReceiveVolumeData
Adds volume data from an Element.
Definition: InterpolatorReceiveVolumeData.hpp:66

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, 40]]
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:

Modifies:

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{surface}}\); 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{mesh}}\), where \(l_{\mathrm{mesh}} > l_{\mathrm{surface}}\). 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{surface}}\) 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: