SpECTRE  v2024.12.16
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
control_system::size Namespace Reference

Classes and functions used in implementation of size control. More...

Classes

struct  ControlErrorArgs
 Packages some of the inputs to the State::control_error, so that State::control_error doesn't need a large number of arguments. More...
 
struct  CrossingTimeInfo
 Holds information about crossing times, as computed by ZeroCrossingPredictors. More...
 
struct  ErrorDiagnostics
 A simple struct to hold diagnostic information about computing the size control error. More...
 
struct  Info
 Holds information that is saved between calls of SizeControl. More...
 
struct  is_size
 
struct  is_size< Systems::Size< Label, DerivOrder > >
 
class  State
 Represents a 'state' of the size control system. More...
 
struct  StateHistory
 A struct for holding a history of control errors for each state in the control_system::Systems::Size control system. More...
 
struct  StateUpdateArgs
 Packages some of the inputs to the State::update, so that State::update doesn't need a large number of arguments. More...
 

Functions

double control_error_delta_r (const double horizon_00, const double dt_horizon_00, const double lambda_00, const double dt_lambda_00, const double grid_frame_excision_sphere_radius)
 Function that computes the control error for control_system::size::States::DeltaR. More...
 
void comoving_char_speed_derivative (gsl::not_null< Scalar< DataVector > * > result, double lambda_00, double dt_lambda_00, double horizon_00, double dt_horizon_00, double grid_frame_excision_sphere_radius, const tnsr::i< DataVector, 3, Frame::Distorted > &excision_rhat, const tnsr::i< DataVector, 3, Frame::Distorted > &excision_normal_one_form, const Scalar< DataVector > &excision_normal_one_form_norm, const tnsr::I< DataVector, 3, Frame::Distorted > &distorted_components_of_grid_shift, const tnsr::II< DataVector, 3, Frame::Distorted > &inverse_spatial_metric_on_excision_boundary, const tnsr::Ijj< DataVector, 3, Frame::Distorted > &spatial_christoffel_second_kind, const tnsr::i< DataVector, 3, Frame::Distorted > &deriv_lapse, const tnsr::iJ< DataVector, 3, Frame::Distorted > &deriv_of_distorted_shift, const InverseJacobian< DataVector, 3, Frame::Grid, Frame::Distorted > &inverse_jacobian_grid_to_distorted)
 Computes the derivative of the comoving characteristic speed with respect to the size map parameter. More...
 
template<typename Frame >
ErrorDiagnostics control_error (const gsl::not_null< Info * > info, const gsl::not_null< intrp::ZeroCrossingPredictor * > predictor_char_speed, const gsl::not_null< intrp::ZeroCrossingPredictor * > predictor_comoving_char_speed, const gsl::not_null< intrp::ZeroCrossingPredictor * > predictor_delta_radius, double time, double control_error_delta_r, std::optional< double > control_error_delta_r_outward, std::optional< double > max_allowed_radial_distance, double dt_lambda_00, const ylm::Strahlkorper< Frame > &apparent_horizon, const ylm::Strahlkorper< Frame > &excision_boundary, const Scalar< DataVector > &lapse_on_excision_boundary, const tnsr::I< DataVector, 3, Frame > &frame_components_of_grid_shift, const tnsr::ii< DataVector, 3, Frame > &spatial_metric_on_excision_boundary, const tnsr::II< DataVector, 3, Frame > &inverse_spatial_metric_on_excision_boundary)
 Computes the size control error, updating the stored info. More...
 
void register_derived_with_charm ()
 
template<size_t DerivOrder, ::domain::ObjectLabel Horizon, typename Metavariables >
void update_averager (const gsl::not_null< Averager< DerivOrder - 1 > * > averager, const gsl::not_null< ControlErrors::Size< DerivOrder, Horizon > * > control_error, const Parallel::GlobalCache< Metavariables > &cache, const double time, const DataVector &current_timescale, const std::string &function_of_time_name, const int current_measurement)
 Updates the Averager with information from the Size control error . More...
 
template<size_t DerivOrder, ::domain::ObjectLabel Horizon, typename Metavariables >
void update_tuner (const gsl::not_null< TimescaleTuner< false > * > tuner, const gsl::not_null< ControlErrors::Size< DerivOrder, Horizon > * > control_error, const Parallel::GlobalCache< Metavariables > &cache, const double time, const std::string &function_of_time_name)
 Updates the TimescaleTuner with information from the Size control error . More...
 

Variables

template<typename T >
constexpr bool is_size_v = is_size<T>::value
 Check whether a control system is the control_system::Systems::Size system or not.
 

Detailed Description

Classes and functions used in implementation of size control.

Function Documentation

◆ comoving_char_speed_derivative()

void control_system::size::comoving_char_speed_derivative ( gsl::not_null< Scalar< DataVector > * >  result,
double  lambda_00,
double  dt_lambda_00,
double  horizon_00,
double  dt_horizon_00,
double  grid_frame_excision_sphere_radius,
const tnsr::i< DataVector, 3, Frame::Distorted > &  excision_rhat,
const tnsr::i< DataVector, 3, Frame::Distorted > &  excision_normal_one_form,
const Scalar< DataVector > &  excision_normal_one_form_norm,
const tnsr::I< DataVector, 3, Frame::Distorted > &  distorted_components_of_grid_shift,
const tnsr::II< DataVector, 3, Frame::Distorted > &  inverse_spatial_metric_on_excision_boundary,
const tnsr::Ijj< DataVector, 3, Frame::Distorted > &  spatial_christoffel_second_kind,
const tnsr::i< DataVector, 3, Frame::Distorted > &  deriv_lapse,
const tnsr::iJ< DataVector, 3, Frame::Distorted > &  deriv_of_distorted_shift,
const InverseJacobian< DataVector, 3, Frame::Grid, Frame::Distorted > &  inverse_jacobian_grid_to_distorted 
)

Computes the derivative of the comoving characteristic speed with respect to the size map parameter.

Parameters
resultthe derivative of the comoving char speed dvc/dλ00, which is computed here using Eq. ( 29).
lambda_00the map parameter λ00. This is the usual spherical harmonic coefficient, not a Spherepack value.
dt_lambda_00the time derivative of the map parameter
horizon_00the average coefficient of the horizon S^00. This is the usual spherical harmonic coefficient, not a Spherepack value.
dt_horizon_00the time derivative of horizon_00
grid_frame_excision_sphere_radiusradius of the excision boundary in the grid frame, rEB.
excision_rhatthe direction cosine ξi^. Not a spacetime tensor: it is raised/lowered with δij
excision_normal_one_formthe unnormalized one-form s^i^
excision_normal_one_form_normthe norm of the one-form a
distorted_components_of_grid_shiftthe quantity βixi^xi evaluated on the excision boundary. This is not the shift in the distorted frame.
inverse_spatial_metric_on_excision_boundarymetric in the distorted frame.
spatial_christoffel_second_kindthe Christoffel symbols Γi^j^k^
deriv_lapsethe spatial derivative of the lapse i^α
deriv_of_distorted_shiftthe spatial derivative of the shift in the distorted frame j^β^i^. This is not the derivative of distorted_components_of_grid_shift.
inverse_jacobian_grid_to_distortedthe quantity Jk^i=k^xi, where xi are the grid frame coordinates and xk^ are the distorted frame coordinates.

Background

The characteristic speed on the excision boundary is

(1)v=α+niβi

where α is the lapse (invariant under frame transformations), βi is the grid-frame shift, and ni is the metric-normalized outward-pointing (i.e. pointing out of the black hole, toward larger radius) normal one-form to the excision boundary in the grid frame. (Note that the usual expression for the characteristic speed, as in eq. 87 of [91], has a minus sign and defines ni as the inward-pointing (i.e. out of the computational domain) normal; here we have a plus sign and we define ni as outward-pointing because the outward-pointing normal is passed into comoving_char_speed_derivative.)

The size/shape map at the excision boundary is given by Eq. 72 of [91] :

(2)x^i=xirEB(1λ00Y00>0Ymλm),

where x^i are the distorted-frame coordinates and xi are the grid-frame coordinates, and where we have separated the =0 piece from the sum. Here Ym are spherical harmonics, λm are the map parameters, and rEB is the radius of the excision boundary in the grid frame (where the excision boundary is a sphere). The final term with the sum over >0 is independent of λ00, and will not be important because below we will be differentiating the map with respect to λ00.

The comoving characteristic speed is given by rewriting Eq. 98 of [91] in terms of the distorted-frame shift:

(3)vc=α+n^i^β^i^Y00n^i^ξi^[S^˙00(λ00rEB/Y00)/S^00+1Y00>0Ymλ˙m],(4)=α+n^i^βi^Y00n^i^ξi^[S^˙00(λ00rEB/Y00)/S^00λ˙00],

where in the last line we have rewritten β^i^ in terms of βi^ (see Eq. ( 15) below) and we have substituted the time derivative of Eq. ( 2). Here λ˙00 is the time derivative of λ00, and S^00 is the constant spherical-harmonic coefficient of the horizon and S^˙00 is its time derivative. The symbol ξi^ is a direction cosine, i.e. xi/rEB evaluated on the excision boundary, which is the same as x^i/r^EB evaluated on the excision boundary because the size and shape maps preserve angles. Note that rEB is a constant but r^EB is a function of angles. Note also that ξi^ is not a vector; it is a coordinate quantity. In particular, the lower-index ξi^ is δijxj/rEB. The non-vectorness of ξi^ (and of xi itself in Eq. ( 2)) might cause some confusion when using the Einstein summation convention; we attempt to alleviate that confusion by never using the lower-index ξi^ and by keeping δij in formulas below. The normal n^i^ is the same as ni transformed into the distorted frame, that is n^i^=njxj/x^i^. We have put a hat on n^ in addition to putting a hat on its index (despite the usual convention that tensors have decorations on indices and not on the tensors themselves) to reduce later ambiguities in notation that arise because Eq. ( 2) has the same index on both sides of the equation and because ξi^ and xi are not tensors. The quantity βi^ in Eq. ( 4) is the distorted-frame component of the grid-frame shift, defined by

(5)βi^=βix^i^xi.

This is not the shift in the distorted frame β^i^, because the shift does not transform like a spatial tensor under the maps.

If the comoving characteristic speed vc is negative and remains negative forever, then size control will fail. Therefore vc is used in making decisions in size control. We care about dvc/dλ00 because the sign of that quantity tells us whether vc will increase or decrease when we increase or decrease the map parameter λ00; this information will be used to decide whether to transition between different size control states.

Derivation of derivative of comoving characteristic speed

This function computes dvc/dλ00 on the excision boundary, where the total derivative means that all other map parameters (like λm for >0) are held fixed, and the coordinates of the excision boundary (the grid coordinates) are held fixed. We also hold fixed λ˙00 because we are interested in how vc changes from a configuration with a given λ00 and λ˙000 to another configuration with a different nearby λ00 and also with λ˙000.

Here we derive an expression for dvc/dλ00. This expression will be complicated, mostly because of the normals n^i^ that appear in Eq. ( 4) and because of the Jacobians.

Derivative of the Jacobian

First, note that by differentiating Eq. ( 2) we obtain

(6)dx^idλ00=xiY00rEB(7)=ξiY00,

where the last line follows from the definition of the direction cosines.

The Jacobian of the map is

(8)x^ixj=(1λ00Y00/rEB+B)δji+xiBxj,

where B represents the term with the sum over >0 in Eq. ( 2); this term is independent of λ00. Therefore, we have

(9)ddλ00x^ixj=Y00rEBδji.

But we want the derivative of the inverse Jacobian, not the forward Jacobian. By taking the derivative of the identity

(10)x^i^xkxkx^j^=δj^i^

and by using Eq. ( 9) we can derive

(11)ddλ00xix^j=+Y00rEBxix^kxkx^j.

Note that the right-hand side of Eq. ( 11) has two inverse Jacobians contracted with each other, which is not the same as δji.

Derivative of a function of space

Assume we have an arbitrary function of space f(x^i). Here we treat f as a function of the distorted-frame coordinates x^i and not a function of the grid-frame coordinates. This is because we consider the metric functions to be defined in the inertial frame (and equivalently for our purposes the functions are defined in the distorted frame because the distorted-to-inertial map is independent of λ00), and we consider λ00 a parameter in a map that moves the grid with respect to these distorted-frame metric functions. The derivative of f can be written

(12)ddλ00f=fx^idx^idλ00(13)=ξi^Y00fx^i.

This is how we will evaluate derivatives of metric functions like the lapse.

Derivative of the distorted-frame components of the grid-frame shift.

To differentiate the quantity defined by Eq. ( 5) note that

(14)βi^=βix^i^xi(15)=β^i^+x^i^t,

where β^i^α2g0^i^ is the shift in the distorted frame. From the map, Eq. ( 2), we see that

(16)ddλ00x^i^t=0,

because there is no remaining λ00 in x^i^t. So

(17)ddλ00βi^=ddλ00β^i^(18)=ξj^Y00j^β^i^,

where we have used Eq. ( 13) in the last line. Note that we cannot use Eq. ( 13) on βi^ directly, because βi^ depends in a complicated way on the grid-to-distorted map. In particular, we will be evaluating j^β^i^ numerically, and numerical spatial derivatives j^β^i^ are not the same as numerical spatial derivatives j^βi^.

Derivative of the normal one-form

The normal to the surface is the most complicated expression in Eq. ( 4), because of how it depends on the map and on the metric. The grid-frame un-normalized outward-pointing one-form to the excision boundary is

(19)si=ξjδij,

because the excision boundary is a sphere of fixed radius in the grid frame. Therefore si doesn't depend on λ00.

The normalized one-form n^i^ is given by

(20)n^i^=s^i^a,

where

(21)s^i^=sixix^i^,(22)a2=s^i^s^j^γi^j^.

Here γi^j^ is the inverse 3-metric in the distorted frame. Again, to avoid ambiguity later, we have put hats on n and s, despite the usual convention that when transforming tensors one puts hats on the indices and not on the tensors.

Now

(23)ddλ00s^i^=Y00rEBs^kxkx^i^,(24)ddλ00a2=2Y00rEBs^kxkx^i^s^j^γi^j^+s^i^s^j^γi^k^γj^l^ξm^Y00m^γk^l^.

Here we have used Eq. ( 11) to differentiate the Jacobian, and Eq. ( 13) to differentiate the 3-metric. We have also refrained from raising and lowering indices on n^i^, s^i^, and ξi^ to alleviate potential confusion over whether to raise or lower using γi^j^ or using δi^j^. The factor s^kxk/x^i^ is unusal and is not a tensor ( s^k is a tensor but the Jacobian it is being multiplied by is the inverse of the one that would transform it into a different frame); this factor arises because some quantities being differentiated are not tensors.

Given the above, the derivative of the normalized normal one-form is

(25)ddλ00n^i^=1addλ00s^i^s^i^12a3ddλ00a2(26)=s^iY00arEBxix^i^s^i^1a3s^iY00rEBxix^k^s^j^γk^j^s^i^Y002a3s^p^s^j^γp^k^γj^l^ξm^m^γk^l^(27)=n^ixix^k^Y00rEB(δi^k^n^k^n^i^)s^i^Y002a3s^p^s^j^γp^k^γj^l^ξm^m^γk^l^(28)=n^ixix^k^Y00rEB(δi^k^n^k^n^i^)Y00n^i^n^p^n^j^γp^k^ξm^Γk^m^j^,

where we have eliminated s^i^ and a in favor of n^i^ and we have substituted 3-Christoffel symbols for spatial derivatives of the 3-metric (and the factor of 2 on the penultimate line has been absorbed into the 3-Christoffel symbol on the last line). Note that the last term in Eq. ( 28) could also be derived by differentiating n^i^n^j^γi^j^=1. The first term in Eq. ( 28) is strange because the inverse Jacobian (as opposed to the forward Jacobian) is contracted with n^i, so that is not a tensor transformation.

We can now differentiate Eq. ( 4) to obtain

ddλ00vc=ξi^Y00i^α+[βi^Y00ξi^S^˙00(λ00rEB/Y00)/S^00+Y00ξi^λ˙00]ddλ00n^i^(29)n^i^ξj^Y00j^β^i^Y00n^i^ξi^S^˙00/S^00,

where ddλ00n^i^ is given by Eq. ( 28).

◆ control_error()

template<typename Frame >
ErrorDiagnostics control_system::size::control_error ( const gsl::not_null< Info * >  info,
const gsl::not_null< intrp::ZeroCrossingPredictor * >  predictor_char_speed,
const gsl::not_null< intrp::ZeroCrossingPredictor * >  predictor_comoving_char_speed,
const gsl::not_null< intrp::ZeroCrossingPredictor * >  predictor_delta_radius,
double  time,
double  control_error_delta_r,
std::optional< double >  control_error_delta_r_outward,
std::optional< double >  max_allowed_radial_distance,
double  dt_lambda_00,
const ylm::Strahlkorper< Frame > &  apparent_horizon,
const ylm::Strahlkorper< Frame > &  excision_boundary,
const Scalar< DataVector > &  lapse_on_excision_boundary,
const tnsr::I< DataVector, 3, Frame > &  frame_components_of_grid_shift,
const tnsr::ii< DataVector, 3, Frame > &  spatial_metric_on_excision_boundary,
const tnsr::II< DataVector, 3, Frame > &  inverse_spatial_metric_on_excision_boundary 
)

Computes the size control error, updating the stored info.

Template Parameters
Frameshould be Frame::Distorted if Frame::Distorted exists.
Parameters
infostruct containing parameters that will be used/filled. Some of the fields in info will guide the behavior of other control system components like the averager and (possibly) the time step.
predictor_char_speedZeroCrossingPredictor for the characteristic speed.
predictor_comoving_char_speedZeroCrossingPredictor for the comoving characteristic speed.
predictor_delta_radiusZeroCrossingPredictor for the difference in radius between the horizon and the excision boundary.
timethe current time.
control_error_delta_rthe control error for the DeltaR state. This is used in other states as well.
control_error_delta_r_outwardthe control error for the DeltaRDriftOutward state. If std::nullopt, then DeltaRDriftOutward will not be used.
max_allowed_radial_distancethe maximum average radial distance between the horizon and the excision boundary that is allowed without triggering the DeltaRDriftOutward state. If std::nullopt, then DeltaRDriftOutward will not be used.
dt_lambda_00the time derivative of the map parameter lambda_00
apparent_horizonthe current horizon in frame Frame.
excision_boundarya Strahlkorper representing the excision boundary in frame Frame. Note that the excision boundary is assumed to be a sphere in the grid frame.
lapse_on_excision_boundaryLapse on the excision boundary.
frame_components_of_grid_shiftThe quantity βixi^xi (see below) evaluated on the excision boundary. This is a tensor in frame Frame.
spatial_metric_on_excision_boundarymetric in frame Frame.
inverse_spatial_metric_on_excision_boundarymetric in frame Frame.

Returns: Returns an ErrorDiagnostics object which, in addition to the actual control error, holds a lot of diagnostic information about how the control error was calculated. This information could be used to print to a file if desired.

The characteristic speed that is needed here is

(30)v=αniβi(31)v=αni^β^i^ni^xi^t(32)v=αni¯β¯i¯ni¯xi¯t(33)v=αni^xi^xiβi,

where we have written many equivalent forms in terms of quantities defined in different frames.

Here α is the lapse, which is invariant under frame transformations, ni, ni^, and ni¯ are the metric-normalized normal one-form to the Strahlkorper in the grid, distorted, and inertial frames, and βi, β^i^, and β¯i¯ are the shift in the grid, distorted, and inertial frames.

Note that we decorate the shift with hats and bars in addition to decorating its index, because the shift transforms in a non-obvious way under frame transformations so it is easy to make mistakes. To be clear, these different shifts are defined by

(34)βi=α2g0i,(35)β^i^=α2g0^i^,(36)β¯i¯=α2g0¯i¯,

where gab is the spacetime metric, and they transform like

(37)β^i^=βixi^xixi^t.

The quantity we pass as frame_components_of_grid_shift is

(38)βixi^xi=β^i^+xi^t(39)=β¯j¯xi^xixixj¯+xi^xj¯xj¯t,

where we have listed several equivalent formulas that involve quantities in different frames.

◆ control_error_delta_r()

double control_system::size::control_error_delta_r ( const double  horizon_00,
const double  dt_horizon_00,
const double  lambda_00,
const double  dt_lambda_00,
const double  grid_frame_excision_sphere_radius 
)

Function that computes the control error for control_system::size::States::DeltaR.

This is helpful to have calculated separately because other control errors may make use of this quantity. The equation for the control error is given in Eq. 96 in [91].

Parameters
horizon_00The l=0,m=0 coefficient of the apparent horizon in the distorted frame.
dt_horizon_00The l=0,m=0 coefficient of the time derivative of the apparent horizon in the distorted frame, where the derivative is taken in the distorted frame as well.
lambda_00The l=0,m=0 component of the function of time for the time dependent map
dt_lambda_00The l=0,m=0 component of the time derivative of the function of time for the time dependent map
grid_frame_excision_sphere_radiusRadius of the excision sphere in the grid frame

◆ update_averager()

template<size_t DerivOrder, ::domain::ObjectLabel Horizon, typename Metavariables >
void control_system::size::update_averager ( const gsl::not_null< Averager< DerivOrder - 1 > * >  averager,
const gsl::not_null< ControlErrors::Size< DerivOrder, Horizon > * >  control_error,
const Parallel::GlobalCache< Metavariables > &  cache,
const double  time,
const DataVector current_timescale,
const std::string function_of_time_name,
const int  current_measurement 
)

Updates the Averager with information from the Size control error .

Details

After we calculate the control error, we check if a discontinuous change has occurred (the internal control_system::size::State changed). If no discontinuous change has occurred, then we do nothing. If one did occur, we get a history of control errors from the Size control error and use that to repopulate the Averager history.

Note
See control_system::Tags::Averager for why the Averager is templated on DerivOrder - 1.

◆ update_tuner()

template<size_t DerivOrder, ::domain::ObjectLabel Horizon, typename Metavariables >
void control_system::size::update_tuner ( const gsl::not_null< TimescaleTuner< false > * >  tuner,
const gsl::not_null< ControlErrors::Size< DerivOrder, Horizon > * >  control_error,
const Parallel::GlobalCache< Metavariables > &  cache,
const double  time,
const std::string function_of_time_name 
)

Updates the TimescaleTuner with information from the Size control error .

Details

We check for a suggested timescale from the Size control error . If one is suggested and it is smaller than the current damping timescale, we set the timescale in the TimescaleTuner to this suggested value. Otherwise, we let the TimescaleTuner adjust the timescale. Regardless of whether a timescale was suggested or not, we always reset the size control error.