Numerical Algorithms

Generic numerical algorithms. More...

## Classes

class  intrp::BarycentricRational
A barycentric rational interpolation class. More...

class  intrp::Irregular< Dim >
Interpolates a Variables onto an arbitrary set of points. More...

class  intrp::RegularGrid< Dim >
Interpolate a Variables from a Mesh onto a regular grid of points. More...

## Enumerations

enum  RootFinder::Verbosity { RootFinder::Verbosity::Silent, RootFinder::Verbosity::Quiet, RootFinder::Verbosity::Verbose, RootFinder::Verbosity::Debug }
The different options for the verbosity of gsl_multiroot. More...

enum  RootFinder::Method { RootFinder::Method::Hybrids, RootFinder::Method::Hybrid, RootFinder::Method::Newton }
The different options for the rootfinding method of gsl_multiroot. More...

enum  RootFinder::StoppingCondition { RootFinder::StoppingCondition::AbsoluteAndRelative, RootFinder::StoppingCondition::Absolute }
The different options for the convergence criterion of gsl_multiroot. More...

## Functions

void find_generalized_eigenvalues (gsl::not_null< DataVector *> eigenvalues_real_part, gsl::not_null< DataVector *> eigenvalues_imaginary_part, gsl::not_null< Matrix *> eigenvectors, Matrix matrix_a, Matrix matrix_b) noexcept
Solve the generalized eigenvalue problem for two matrices. More...

template<size_t Dim>
double definite_integral (const DataVector &integrand, const Mesh< Dim > &mesh) noexcept
Compute the definite integral of a grid-function over a manifold. More...

template<typename FluxTags , size_t Dim, typename DerivativeFrame >
auto divergence (const Variables< FluxTags > &F, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept -> Variables< db::wrap_tags_in< Tags::div, FluxTags >>
Compute the (Euclidean) divergence of fluxes.

template<size_t Dim>
double mean_value (const DataVector &f, const Mesh< Dim > &mesh) noexcept
Compute the mean value of a grid-function over a manifold. $mean value = \int f dV / \int dV$. More...

template<typename T >
void raw_transpose (const gsl::not_null< T *> result, const T *const data, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data. More...

template<typename Function >
double RootFinder::newton_raphson (const Function &f, const double initial_guess, const double lower_bound, const double upper_bound, const size_t digits, const size_t max_iterations=50)
Finds the root of the function f with the Newton-Raphson method. More...

template<typename Function >
DataVector RootFinder::newton_raphson (const Function &f, const DataVector &initial_guess, const DataVector &lower_bound, const DataVector &upper_bound, const size_t digits, const size_t max_iterations=50)
Finds the root of the function f with the Newton-Raphson method on each element in a DataVector. More...

double positive_root (double a, double b, double c)
Returns the positive root of a quadratic equation ax^2 + bx + c = 0. More...

std::array< double, 2 > real_roots (double a, double b, double c)
Returns the two real roots of a quadratic equation ax^2 + bx + c = 0. More...

template<typename Function >
double RootFinder::toms748 (const Function &f, const double lower_bound, const double upper_bound, const double absolute_tolerance, const double relative_tolerance, const size_t max_iterations=100)
Finds the root of the function f with the TOMS_748 method. More...

template<typename Function >
DataVector RootFinder::toms748 (const Function &f, const DataVector &lower_bound, const DataVector &upper_bound, const double absolute_tolerance, const double relative_tolerance, const size_t max_iterations=100)
Finds the root of the function f with the TOMS_748 method on each element in a DataVector. More...

template<typename VariableTags , typename MatrixType , size_t Dim>
void apply_matrices (const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension. More...

template<typename VariableTags , typename MatrixType , size_t Dim>
Variables< VariableTags > apply_matrices (const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension. More...

template<typename MatrixType , size_t Dim>
void apply_matrices (const gsl::not_null< DataVector *> result, const std::array< MatrixType, Dim > &matrices, const DataVector &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension. More...

template<typename MatrixType , size_t Dim>
DataVector apply_matrices (const std::array< MatrixType, Dim > &matrices, const DataVector &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension. More...

template<size_t Dim>
void linearize (gsl::not_null< DataVector *> result, const DataVector &u, const Mesh< Dim > &mesh) noexcept
Truncate u to a linear function in each dimension. More...

template<size_t Dim>
DataVector linearize (const DataVector &u, const Mesh< Dim > &mesh) noexcept
Truncate u to a linear function in each dimension. More...

template<size_t Dim>
void linearize (gsl::not_null< DataVector *> result, const DataVector &u, const Mesh< Dim > &mesh, size_t d) noexcept
Truncate u to a linear function in the given dimension. More...

template<size_t Dim>
DataVector linearize (const DataVector &u, const Mesh< Dim > &mesh, size_t d) noexcept
Truncate u to a linear function in the given dimension. More...

template<size_t Dim>
double mean_value_on_boundary (const DataVector &f, const Mesh< Dim > &mesh, size_t d, Side side) noexcept

template<size_t Dim>
double mean_value_on_boundary (gsl::not_null< DataVector *> boundary_buffer, const DataVector &f, const Mesh< Dim > &mesh, size_t d, Side side) noexcept

double mean_value_on_boundary (gsl::not_null< DataVector *>, const DataVector &f, const Mesh< 1 > &mesh, size_t d, Side side) noexcept

template<size_t Dim>
double mean_value_on_boundary (gsl::not_null< DataVector *> boundary_buffer, gsl::span< std::pair< size_t, size_t >> volume_and_slice_indices, const DataVector &f, const Mesh< Dim > &mesh, size_t d, Side) noexcept

double mean_value_on_boundary (gsl::not_null< DataVector *>, gsl::span< std::pair< size_t, size_t >>, const DataVector &f, const Mesh< 1 > &mesh, size_t d, Side side) noexcept

template<typename DerivativeTags , typename VariableTags , size_t Dim>
void logical_partial_derivatives (gsl::not_null< std::array< Variables< DerivativeTags >, Dim > *> logical_partial_derivatives_of_u, const Variables< VariableTags > &u, const Mesh< Dim > &mesh) noexcept
Compute the partial derivatives of each variable with respect to the logical coordinate. More...

template<typename DerivativeTags , typename VariableTags , size_t Dim>
auto logical_partial_derivatives (const Variables< VariableTags > &u, const Mesh< Dim > &mesh) noexcept -> std::array< Variables< DerivativeTags >, Dim >
Compute the partial derivatives of each variable with respect to the logical coordinate. More...

template<typename DerivativeTags , size_t Dim, typename DerivativeFrame >
void partial_derivatives (gsl::not_null< Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >> *> du, const std::array< Variables< DerivativeTags >, Dim > &logical_partial_derivatives_of_u, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept
Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame. More...

template<typename DerivativeTags , typename VariableTags , size_t Dim, typename DerivativeFrame >
void partial_derivatives (gsl::not_null< Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >> *> du, const Variables< VariableTags > &u, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept
Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame. More...

template<typename DerivativeTags , typename VariableTags , size_t Dim, typename DerivativeFrame >
auto partial_derivatives (const Variables< VariableTags > &u, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept -> Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >>
Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame. More...

template<typename U , typename T >
void transpose (const gsl::not_null< T *> result, const U &u, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data. More...

template<typename U , typename T = U>
transpose (const U &u, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data. More...

template<size_t Dim, typename Function , Requires< gsl_multiroot_detail::is_jacobian_callable_v< Function, std::array< double, Dim >>> = nullptr>
std::array< double, Dim > RootFinder::gsl_multiroot (const Function &func, const std::array< double, Dim > &initial_guess, const double absolute_tolerance, const size_t maximum_iterations, const double relative_tolerance=0.0, const Verbosity verbosity=Verbosity::Silent, const double maximum_absolute_tolerance=0.0, const Method method=Method::Newton, const StoppingCondition condition=StoppingCondition::Absolute)
A multidimensional root finder supporting Newton and Hybrid methods, as well as modified methods based on these. More...

## Detailed Description

Generic numerical algorithms.

## ◆ Method

 enum RootFinder::Method
strong

The different options for the rootfinding method of gsl_multiroot.

This enum is for setting the method used the rootfinder. The precise method used by the gsl rootfinder depends on whether or not the function passed to it has a callable jacobian member function. In the case where it doesn't, the jacobian is approximated with a finite difference. For example, if the Method specified is Hybrid, gsl will use the gsl_multiroot_fdfsolver_hybridj method in the case where a jacobian is provided, and gsl_multiroot_fsolver_hybrid in the case where one isn't. See GSL's documentation for multidimensional rootfinding for information on the different methods.

Note
gsl does not provide a finite difference version for the modified Newton method (gsl_multiroot_fdfsolver_gnewton). In the case where a jacobian is not provided the method used will be a non-modified Newton method.
Enumerator
Hybrids

Hybrid of Newton's method along with following the gradient direction.

Note
Sometimes Hybrids works only with the Absolute stopping condition.
Hybrid

"Unscaled version of Hybrids that uses a spherical trust region," see GSL documentation for more details.

Newton

If an analytic jacobian is provided, gsl uses a modification of Newton's method to improve global convergence. Uses vanilla Newton's method if no jacobian is provided.

## ◆ StoppingCondition

 strong

The different options for the convergence criterion of gsl_multiroot.

See GSL's documentation for multidimensional rootfinding for information on the different stopping conditions.

Enumerator
AbsoluteAndRelative

See GSL documentation for gsl_multiroot_test_delta.

Absolute

See GSL documentation for gsl_multiroot_test_residual.

## ◆ Verbosity

 strong

The different options for the verbosity of gsl_multiroot.

Enumerator
Silent

Do not print anything.

Quiet

Print only "success" or "failed" on termination.

Verbose

Print final functions values on termination.

Debug

Print function values on every iteration.

## ◆ apply_matrices() [1/4]

template<typename VariableTags , typename MatrixType , size_t Dim>
 void apply_matrices ( const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > & matrices, const Variables< VariableTags > & u, const Index< Dim > & extents )
noexcept

Multiply by matrices in each dimension.

Multiplies each stripe in the first dimension of u by matrices[0], each stripe in the second dimension of u by matrices[1], and so on. If any of the matrices are empty they will be treated as the identity, but the matrix multiplications will be skipped for increased efficiency.

## ◆ apply_matrices() [2/4]

template<typename VariableTags , typename MatrixType , size_t Dim>
 Variables apply_matrices ( const std::array< MatrixType, Dim > & matrices, const Variables< VariableTags > & u, const Index< Dim > & extents )
noexcept

Multiply by matrices in each dimension.

Multiplies each stripe in the first dimension of u by matrices[0], each stripe in the second dimension of u by matrices[1], and so on. If any of the matrices are empty they will be treated as the identity, but the matrix multiplications will be skipped for increased efficiency.

## ◆ apply_matrices() [3/4]

template<typename MatrixType , size_t Dim>
 void apply_matrices ( const gsl::not_null< DataVector *> result, const std::array< MatrixType, Dim > & matrices, const DataVector & u, const Index< Dim > & extents )
noexcept

Multiply by matrices in each dimension.

Multiplies each stripe in the first dimension of u by matrices[0], each stripe in the second dimension of u by matrices[1], and so on. If any of the matrices are empty they will be treated as the identity, but the matrix multiplications will be skipped for increased efficiency.

## ◆ apply_matrices() [4/4]

template<typename MatrixType , size_t Dim>
 DataVector apply_matrices ( const std::array< MatrixType, Dim > & matrices, const DataVector & u, const Index< Dim > & extents )
noexcept

Multiply by matrices in each dimension.

Multiplies each stripe in the first dimension of u by matrices[0], each stripe in the second dimension of u by matrices[1], and so on. If any of the matrices are empty they will be treated as the identity, but the matrix multiplications will be skipped for increased efficiency.

## ◆ definite_integral()

template<size_t Dim>
 double definite_integral ( const DataVector & integrand, const Mesh< Dim > & mesh )
noexcept

Compute the definite integral of a grid-function over a manifold.

The integral is computed on the reference element by multiplying the DataVector with the Spectral::quadrature_weights() in that dimension.

Requires: number of points in integrand and mesh are equal.

Parameters
 integrand the grid function to integrate. mesh the Mesh of the manifold on which integrand is located.

Returns: the definite integral of integrand on the manifold.

## ◆ find_generalized_eigenvalues()

 void find_generalized_eigenvalues ( gsl::not_null< DataVector *> eigenvalues_real_part, gsl::not_null< DataVector *> eigenvalues_imaginary_part, gsl::not_null< Matrix *> eigenvectors, Matrix matrix_a, Matrix matrix_b )
noexcept

Solve the generalized eigenvalue problem for two matrices.

This function uses the lapack routine dggev (http://www.netlib.org/lapack/explore-3.1.1-html/dggev.f.html) to solve the generalized eigenvalue problem $A v_a =\lambda_a B v_a$ for the generalized eigenvalues $\lambda_a$ and corresponding eigenvectors $v_a$. matrix_a and matrix_b are each a Matrix; they correspond to square matrices $A$ and $B$ that are the same dimension $N$. eigenvalues_real_part is a DataVector of size $N$ that will store the real parts of the eigenvalues, eigenvalues_imaginary_part is a DataVector of size $N$ that will store the imaginary parts of the eigenvalues. Complex eigenvalues always form complex conjugate pairs, and the $j$ and $j+1$ eigenvalues will have the forms $a+ib$ and $a-ib$, respectively. The eigenvectors are returned as the columns of a square Matrix of dimension $N$ called eigenvectors. If eigenvalue $j$ is real, then column $j$ of eigenvectors is the corresponding eigenvector. If eigenvalue $j$ and $j+1$ are complex-conjugate pairs, then the eigenvector for eigenvalue $j$ is (column j) + $i$ (column j+1), and the eigenvector for eigenvalue $j+1$ is (column j) - $i$ (column j+1).

## ◆ gsl_multiroot()

template<size_t Dim, typename Function , Requires< gsl_multiroot_detail::is_jacobian_callable_v< Function, std::array< double, Dim >>> = nullptr>
 std::array< double, Dim > RootFinder::gsl_multiroot ( const Function & func, const std::array< double, Dim > & initial_guess, const double absolute_tolerance, const size_t maximum_iterations, const double relative_tolerance = 0.0, const Verbosity verbosity = Verbosity::Silent, const double maximum_absolute_tolerance = 0.0, const Method method = Method::Newton, const StoppingCondition condition = StoppingCondition::Absolute )

A multidimensional root finder supporting Newton and Hybrid methods, as well as modified methods based on these.

This root finder accepts function objects with and without a callable jacobian member function. This member function both accepts and returns a std::array<double, Dim>, the dimension of the domain and range of the function the root find is being performed on. Whether the jacobian is provided determines the details of the implementation of the root-finding method that is selected by the user using the Method enum. That is, whether the jacobian is computed analytically via the jacobian member function, or whether the jacobian is computed numerically via a finite difference approximation.

Note
GSL does not provide a finite difference version of its modified Newton method, so the unmodified one is used instead when the user uses the Method::Newton method.

The user can select one of two possible criteria for convergence, StoppingCondition::Absolute, where the sum of the absolute values of the components of the residual vector f are compared against the value provided to absolute_tolerance, and StoppingCondition::AbsoluteAndRelative, where the size of the most recent step taken in the root-finding iteration is compared against absolute_tolerance + relative_tolerance * |x_i|, for each component. In either case, a maximum_absolute_tolerance may be specified if the user anticipates that the convergence criterion specified with StoppingCondition will be too strict for a few points out of a population of points found with a sequence of root finds.

See GSL's documentation for multidimensional rootfinding for reference.

Parameters
 func Function whose root is to be found. initial_guess Contains initial guess. absolute_tolerance The absolute tolerance. maximum_iterations The maximum number of iterations. relative_tolerance The relative tolerance. verbosity Whether to print diagnostic messages. maximum_absolute_tolerance Acceptable absolute tolerance when root finder doesn't converge. You may wish to use this if there are only a few "problematic" points where it is difficult to do a precise root find. method The method to use. See the documentation for the Method enum. condition The convergence condition to use. See the documentation for the StoppingCondition enum.

## ◆ linearize() [1/4]

template<size_t Dim>
 void linearize ( gsl::not_null< DataVector *> result, const DataVector & u, const Mesh< Dim > & mesh )
noexcept

Truncate u to a linear function in each dimension.

Ex in 2D: $u^{Lin} = U_0 + U_x x + U_y y + U_{xy} xy$

Warning
the gsl::not_null variant assumes *result is of the correct size.

## ◆ linearize() [2/4]

template<size_t Dim>
 DataVector linearize ( const DataVector & u, const Mesh< Dim > & mesh )
noexcept

Truncate u to a linear function in each dimension.

Ex in 2D: $u^{Lin} = U_0 + U_x x + U_y y + U_{xy} xy$

Warning
the gsl::not_null variant assumes *result is of the correct size.

## ◆ linearize() [3/4]

template<size_t Dim>
 void linearize ( gsl::not_null< DataVector *> result, const DataVector & u, const Mesh< Dim > & mesh, size_t d )
noexcept

Truncate u to a linear function in the given dimension.

Parameters

• u the function to linearize.
• mesh the Mesh of the grid on the manifold on which u is located.
• d the dimension that is to be linearized.
Warning
the gsl::not_null variant assumes *result is of the correct size.

## ◆ linearize() [4/4]

template<size_t Dim>
 DataVector linearize ( const DataVector & u, const Mesh< Dim > & mesh, size_t d )
noexcept

Truncate u to a linear function in the given dimension.

Parameters

• u the function to linearize.
• mesh the Mesh of the grid on the manifold on which u is located.
• d the dimension that is to be linearized.
Warning
the gsl::not_null variant assumes *result is of the correct size.

## ◆ logical_partial_derivatives() [1/2]

template<typename DerivativeTags , typename VariableTags , size_t Dim>
 void logical_partial_derivatives ( gsl::not_null< std::array< Variables< DerivativeTags >, Dim > *> logical_partial_derivatives_of_u, const Variables< VariableTags > & u, const Mesh< Dim > & mesh )
noexcept

Compute the partial derivatives of each variable with respect to the logical coordinate.

Requires: DerivativeTags to be the head of VariableTags

Returns a Variables with a spatial tensor index appended to the front of each tensor within u and each Tag wrapped with a Tags::deriv.

Template Parameters
 DerivativeTags the subset of VariableTags for which derivatives are computed.

## ◆ logical_partial_derivatives() [2/2]

template<typename DerivativeTags , typename VariableTags , size_t Dim>
 auto logical_partial_derivatives ( const Variables< VariableTags > & u, const Mesh< Dim > & mesh ) -> std::array< Variables< DerivativeTags >, Dim >
noexcept

Compute the partial derivatives of each variable with respect to the logical coordinate.

Requires: DerivativeTags to be the head of VariableTags

Returns a Variables with a spatial tensor index appended to the front of each tensor within u and each Tag wrapped with a Tags::deriv.

Template Parameters
 DerivativeTags the subset of VariableTags for which derivatives are computed.

## ◆ mean_value()

template<size_t Dim>
 double mean_value ( const DataVector & f, const Mesh< Dim > & mesh )
noexcept

Compute the mean value of a grid-function over a manifold. $mean value = \int f dV / \int dV$.

Remarks: The mean value is computed on the reference element(s).

Note
The mean w.r.t. a different set of coordinates x can be computed by pre-multiplying the argument f by the Jacobian J = dx/dxi of the mapping from the reference coordinates xi to the coordinates x.

Returns: the mean value of f on the manifold

Parameters
 f the grid function of which to find the mean. mesh the Mesh of the manifold on which f is located.

## ◆ mean_value_on_boundary() [1/5]

template<size_t Dim>
 double mean_value_on_boundary ( const DataVector & f, const Mesh< Dim > & mesh, size_t d, Side side )
noexcept

Compute the mean value of a grid-function on a boundary of a manifold. $mean value = \int f dV / \int dV$

Remarks: The mean value is computed on the reference element(s).

Returns: the mean value of f on the boundary of the manifold

• f the grid function of which to find the mean.
• mesh the Mesh of the manifold on which f is located.
• d the dimension which is sliced away to get the boundary.
• side whether it is the lower or upper boundary in the d-th dimension.
• boundary_buffer is a pointer to a DataVector of size mesh.slice_away(d).number_of_grid_points() used as a temporary buffer when slicing the data to the boundary.
• volume_and_slice_indices a pair of (volume_index_for_point, slice_index_for_point) computed using the SliceIterator. Because SliceIterator is somewhat expensive, if computing the mean value on the same boundary for many different tensor components, prefer computing the slice indices once.

## ◆ mean_value_on_boundary() [2/5]

template<size_t Dim>
 double mean_value_on_boundary ( gsl::not_null< DataVector *> boundary_buffer, const DataVector & f, const Mesh< Dim > & mesh, size_t d, Side side )
noexcept

Compute the mean value of a grid-function on a boundary of a manifold. $mean value = \int f dV / \int dV$

Remarks: The mean value is computed on the reference element(s).

Returns: the mean value of f on the boundary of the manifold

• f the grid function of which to find the mean.
• mesh the Mesh of the manifold on which f is located.
• d the dimension which is sliced away to get the boundary.
• side whether it is the lower or upper boundary in the d-th dimension.
• boundary_buffer is a pointer to a DataVector of size mesh.slice_away(d).number_of_grid_points() used as a temporary buffer when slicing the data to the boundary.
• volume_and_slice_indices a pair of (volume_index_for_point, slice_index_for_point) computed using the SliceIterator. Because SliceIterator is somewhat expensive, if computing the mean value on the same boundary for many different tensor components, prefer computing the slice indices once.

## ◆ mean_value_on_boundary() [3/5]

 double mean_value_on_boundary ( gsl::not_null< DataVector *> , const DataVector & f, const Mesh< 1 > & mesh, size_t d, Side side )
noexcept

Compute the mean value of a grid-function on a boundary of a manifold. $mean value = \int f dV / \int dV$

Remarks: The mean value is computed on the reference element(s).

Returns: the mean value of f on the boundary of the manifold

• f the grid function of which to find the mean.
• mesh the Mesh of the manifold on which f is located.
• d the dimension which is sliced away to get the boundary.
• side whether it is the lower or upper boundary in the d-th dimension.
• boundary_buffer is a pointer to a DataVector of size mesh.slice_away(d).number_of_grid_points() used as a temporary buffer when slicing the data to the boundary.
• volume_and_slice_indices a pair of (volume_index_for_point, slice_index_for_point) computed using the SliceIterator. Because SliceIterator is somewhat expensive, if computing the mean value on the same boundary for many different tensor components, prefer computing the slice indices once.

## ◆ mean_value_on_boundary() [4/5]

template<size_t Dim>
 double mean_value_on_boundary ( gsl::not_null< DataVector *> boundary_buffer, gsl::span< std::pair< size_t, size_t >> volume_and_slice_indices, const DataVector & f, const Mesh< Dim > & mesh, size_t d, Side )
noexcept

Compute the mean value of a grid-function on a boundary of a manifold. $mean value = \int f dV / \int dV$

Remarks: The mean value is computed on the reference element(s).

Returns: the mean value of f on the boundary of the manifold

• f the grid function of which to find the mean.
• mesh the Mesh of the manifold on which f is located.
• d the dimension which is sliced away to get the boundary.
• side whether it is the lower or upper boundary in the d-th dimension.
• boundary_buffer is a pointer to a DataVector of size mesh.slice_away(d).number_of_grid_points() used as a temporary buffer when slicing the data to the boundary.
• volume_and_slice_indices a pair of (volume_index_for_point, slice_index_for_point) computed using the SliceIterator. Because SliceIterator is somewhat expensive, if computing the mean value on the same boundary for many different tensor components, prefer computing the slice indices once.

## ◆ mean_value_on_boundary() [5/5]

 double mean_value_on_boundary ( gsl::not_null< DataVector *> , gsl::span< std::pair< size_t, size_t >> , const DataVector & f, const Mesh< 1 > & mesh, size_t d, Side side )
noexcept

Compute the mean value of a grid-function on a boundary of a manifold. $mean value = \int f dV / \int dV$

Remarks: The mean value is computed on the reference element(s).

Returns: the mean value of f on the boundary of the manifold

• f the grid function of which to find the mean.
• mesh the Mesh of the manifold on which f is located.
• d the dimension which is sliced away to get the boundary.
• side whether it is the lower or upper boundary in the d-th dimension.
• boundary_buffer is a pointer to a DataVector of size mesh.slice_away(d).number_of_grid_points() used as a temporary buffer when slicing the data to the boundary.
• volume_and_slice_indices a pair of (volume_index_for_point, slice_index_for_point) computed using the SliceIterator. Because SliceIterator is somewhat expensive, if computing the mean value on the same boundary for many different tensor components, prefer computing the slice indices once.

## ◆ newton_raphson() [1/2]

template<typename Function >
 double RootFinder::newton_raphson ( const Function & f, const double initial_guess, const double lower_bound, const double upper_bound, const size_t digits, const size_t max_iterations = 50 )

Finds the root of the function f with the Newton-Raphson method.

f is a unary invokable that takes a double which is the current value at which to evaluate f. An example is below.

const size_t digits = 8;
const double correct = sqrt(2.);
const double guess = 1.5;
const double lower = 1.;
const double upper = 2.;
const auto func_and_deriv_lambda = [](double x) noexcept {
return std::make_pair(2. - square(x), -2. * x);
};
const auto root_from_lambda = RootFinder::newton_raphson(
func_and_deriv_lambda, guess, lower, upper, digits);

See the Boost documentation for more details.

Requires: Function f is invokable with a double

Note
The parameter digits specifies the precision of the result in its desired number of base-10 digits.

Throws: convergence_error if the requested precision is not met after max_iterations iterations.

## ◆ newton_raphson() [2/2]

template<typename Function >
 DataVector RootFinder::newton_raphson ( const Function & f, const DataVector & initial_guess, const DataVector & lower_bound, const DataVector & upper_bound, const size_t digits, const size_t max_iterations = 50 )

Finds the root of the function f with the Newton-Raphson method on each element in a DataVector.

f is a binary invokable that takes a double as its first argument and a size_t as its second. The double is the current value at which to evaluate f, and the size_t is the current index into the DataVectors. Below is an example of how to root find different functions by indexing into a lambda-captured DataVector using the size_t passed to f.

const size_t digits = 8;
const DataVector guess{1.6, 1.9, -1.6, -1.9};
const DataVector lower{sqrt(2.), sqrt(2.), -2., -3.};
const DataVector upper{2., 3., -sqrt(2.), -sqrt(2.)};
const DataVector constant{2., 4., 2., 4.};
const auto func_and_deriv_lambda = [&constant](const double x,
const size_t i) noexcept {
return std::make_pair(constant[i] - square(x), -2. * x);
};
const auto root = RootFinder::newton_raphson(func_and_deriv_lambda, guess,
lower, upper, digits);

See the Boost documentation for more details.

Requires: Function f be callable with a double and a size_t

Note
The parameter digits specifies the precision of the result in its desired number of base-10 digits.

Throws: convergence_error if, for any index, the requested precision is not met after max_iterations iterations.

## ◆ partial_derivatives() [1/3]

template<typename DerivativeTags , size_t Dim, typename DerivativeFrame >
 void partial_derivatives ( gsl::not_null< Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >> *> du, const std::array< Variables< DerivativeTags >, Dim > & logical_partial_derivatives_of_u, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > & inverse_jacobian )
noexcept

Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame.

Requires: DerivativeTags to be the head of VariableTags

Returns a Variables with a spatial tensor index appended to the front of each tensor within u and each Tag wrapped with a Tags::deriv.

Template Parameters
 DerivativeTags the subset of VariableTags for which derivatives are computed.

## ◆ partial_derivatives() [2/3]

template<typename DerivativeTags , typename VariableTags , size_t Dim, typename DerivativeFrame >
 void partial_derivatives ( gsl::not_null< Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >> *> du, const Variables< VariableTags > & u, const Mesh< Dim > & mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > & inverse_jacobian )
noexcept

Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame.

Requires: DerivativeTags to be the head of VariableTags

Returns a Variables with a spatial tensor index appended to the front of each tensor within u and each Tag wrapped with a Tags::deriv.

Template Parameters
 DerivativeTags the subset of VariableTags for which derivatives are computed.

## ◆ partial_derivatives() [3/3]

template<typename DerivativeTags , typename VariableTags , size_t Dim, typename DerivativeFrame >
 auto partial_derivatives ( const Variables< VariableTags > & u, const Mesh< Dim > & mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > & inverse_jacobian ) -> Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame >>
noexcept

Compute the partial derivatives of each variable with respect to the coordinates of DerivativeFrame.

Requires: DerivativeTags to be the head of VariableTags

Returns a Variables with a spatial tensor index appended to the front of each tensor within u and each Tag wrapped with a Tags::deriv.

Template Parameters
 DerivativeTags the subset of VariableTags for which derivatives are computed.

## ◆ positive_root()

 double positive_root ( double a, double b, double c )

Returns the positive root of a quadratic equation ax^2 + bx + c = 0.

Returns: The positive root of a quadratic equation.

Requires: That there are two real roots, of which only one is positive.

## ◆ raw_transpose()

template<typename T >
 void raw_transpose ( const gsl::not_null< T *> result, const T *const data, const size_t chunk_size, const size_t number_of_chunks )
noexcept

Function to compute transposed data.

Transpose the data pointed to by data, writing the result to the location pointed to by result. See the transpose function for a safer interface and for the meaning of the other arguments.

## ◆ real_roots()

 std::array real_roots ( double a, double b, double c )

Returns the two real roots of a quadratic equation ax^2 + bx + c = 0.

Returns: An array of the roots of a quadratic equation

Requires: That there are two real roots.

## ◆ toms748() [1/2]

template<typename Function >
 double RootFinder::toms748 ( const Function & f, const double lower_bound, const double upper_bound, const double absolute_tolerance, const double relative_tolerance, const size_t max_iterations = 100 )

Finds the root of the function f with the TOMS_748 method.

f is a unary invokable that takes a double which is the current value at which to evaluate f. An example is below.

const double abs_tol = 1e-15;
const double rel_tol = 1e-15;
const double upper = 2.0;
const double lower = sqrt(2.0) - abs_tol; // bracket surrounds root
const auto f_lambda = [](double x) { return 2.0 - square(x); };
auto root = RootFinder::toms748(f_lambda, lower, upper, abs_tol, rel_tol);

The TOMS_748 algorithm searches for a root in the interval [lower_bound, upper_bound], and will throw if this interval does not bracket a root, i.e. if f(lower_bound) * f(upper_bound) > 0.

See the Boost documentation for more details.

Requires: Function f is invokable with a double

Throws: std::domain_error if the bounds do not bracket a root.

Throws: convergence_error if the requested tolerance is not met after max_iterations iterations.

## ◆ toms748() [2/2]

template<typename Function >
 DataVector RootFinder::toms748 ( const Function & f, const DataVector & lower_bound, const DataVector & upper_bound, const double absolute_tolerance, const double relative_tolerance, const size_t max_iterations = 100 )

Finds the root of the function f with the TOMS_748 method on each element in a DataVector.

f is a binary invokable that takes a double as its first argument and a size_t as its second. The double is the current value at which to evaluate f, and the size_t is the current index into the DataVectors. Below is an example of how to root find different functions by indexing into a lambda-captured DataVector using the size_t passed to f.

const double abs_tol = 1e-15;
const double rel_tol = 1e-15;
const DataVector upper{2.0, 3.0, -sqrt(2.0) + abs_tol, -sqrt(2.0)};
const DataVector lower{sqrt(2.0) - abs_tol, sqrt(2.0), -2.0, -3.0};
const DataVector constant{2.0, 4.0, 2.0, 4.0};
const auto f_lambda = [&constant](const double x, const size_t i) noexcept {
return constant[i] - square(x);
};
const auto root =
RootFinder::toms748(f_lambda, lower, upper, abs_tol, rel_tol);

For each index i into the DataVector, the TOMS_748 algorithm searches for a root in the interval [lower_bound[i], upper_bound[i]], and will throw if this interval does not bracket a root, i.e. if f(lower_bound[i], i) * f(upper_bound[i], i) > 0.

See the Boost documentation for more details.

Requires: Function f be callable with a double and a size_t

Throws: std::domain_error if, for any index, the bounds do not bracket a root.

Throws: convergence_error if, for any index, the requested tolerance is not met after max_iterations iterations.

## ◆ transpose() [1/2]

template<typename U , typename T >
 void transpose ( const gsl::not_null< T *> result, const U & u, const size_t chunk_size, const size_t number_of_chunks )
noexcept

Function to compute transposed data.

The primary use of this function is to rearrange the memory layout so that another function can operate on contiguous chunks of memory.

Requires: result.size() to be the product of number_of_chunks and chunk_size, u.size() to be equal or greater than result.size(), and that both result and u have a data() member function.

### Details

The container u holds a contiguous array of data, treated as a sequence of number_of_chunks contiguous sets of entries of size chunk_size. The output result has its data arranged such that the first number_of_chunks elements in result will be the first element of each chunk of u. The last number_of_chunks elements in result will be the last (i.e. chunk_size-th) element of each chunk of u. If u.size() is greater than result.size() the extra elements of u are ignored.

Note
This is equivalent to treating the first part of u as a matrix and filling result (or the returned object) with the transpose of that matrix.

If u represents a block of data indexed by $(x, y, z, \ldots)$ with the first index varying fastest, transpose serves to rotate the indices. If the extents are $(X, Y, Z, \ldots)$, with product $N$, transpose(u, X, N/X) reorders the data to be indexed $(y, z, \ldots, x)$, transpose(u, X*Y, N/X/Y) reorders the data to be indexed $(z, \ldots, x, y)$, etc.

### Example

const DataVector matrix{ 1., 2., 3.,
4., 5., 6.,
7., 8., 9.,
10., 11., 12.};
CHECK(transpose(matrix, 3, 4) == DataVector{1., 4., 7., 10.,
2., 5., 8., 11.,
3., 6., 9., 12.});
const size_t chunk_size_vars = 8;
const size_t n_grid_pts = 2 * chunk_size_vars;
Variables<two_vars<2>> variables(n_grid_pts, 0.);
for (size_t i = 0; i < variables.size(); ++i) {
// clang-tidy: pointer arithmetic
variables.data()[i] = i * i; // NOLINT
}
const size_t number_of_chunks_vars = variables.size() / chunk_size_vars;
auto transposed_vars = variables;
transpose(make_not_null(&transposed_vars), variables, chunk_size_vars,
number_of_chunks_vars);
for (size_t i = 0; i < chunk_size_vars; ++i) {
for (size_t j = 0; j < number_of_chunks_vars; ++j) {
// clang-tidy: pointer arithmetic
CHECK(variables.data()[i + chunk_size_vars * j] == // NOLINT
transposed_vars.data()[j + number_of_chunks_vars * i]); // NOLINT
}
}
const size_t chunk_size = 8;
const size_t number_of_chunks = 2;
const size_t n_pts = chunk_size * number_of_chunks;
DataVector data(n_pts);
for (size_t i = 0; i < data.size(); ++i) {
data[i] = i * i;
}
DataVector transposed_data(n_pts, 0.);
transposed_data = transpose(data, chunk_size, number_of_chunks);
for (size_t i = 0; i < chunk_size; ++i) {
for (size_t j = 0; j < number_of_chunks; ++j) {
CHECK(data[i + chunk_size * j] ==
transposed_data[j + number_of_chunks * i]);
}
}
Variables<one_var<2>> partial_vars(n_grid_pts, 0.);
get<Var1<2>>(partial_vars) = get<Var1<2>>(variables);
Variables<one_var<2>> partial_transpose(n_grid_pts, 0.);
const size_t partial_number_of_chunks = 2*number_of_chunks_vars / 3;
transpose(make_not_null(&partial_transpose), variables, chunk_size_vars,
partial_number_of_chunks);
for (size_t i = 0; i < chunk_size_vars; ++i) {
for (size_t j = 0; j < partial_number_of_chunks; ++j) {
// clang-tidy: pointer arithmetic
CHECK(partial_transpose
.data()[j + partial_number_of_chunks * i] == // NOLINT
variables.data()[i + chunk_size_vars * j]); // NOLINT
CHECK(partial_transpose
.data()[j + partial_number_of_chunks * i] == // NOLINT
partial_vars.data()[i + chunk_size_vars * j]); // NOLINT
}
}
Template Parameters
 U the type of data to be transposed T the type of the transposed data

## ◆ transpose() [2/2]

template<typename U , typename T = U>
 T transpose ( const U & u, const size_t chunk_size, const size_t number_of_chunks )
noexcept

Function to compute transposed data.

The primary use of this function is to rearrange the memory layout so that another function can operate on contiguous chunks of memory.

Requires: result.size() to be the product of number_of_chunks and chunk_size, u.size() to be equal or greater than result.size(), and that both result and u have a data() member function.

### Details

The container u holds a contiguous array of data, treated as a sequence of number_of_chunks contiguous sets of entries of size chunk_size. The output result has its data arranged such that the first number_of_chunks elements in result will be the first element of each chunk of u. The last number_of_chunks elements in result will be the last (i.e. chunk_size-th) element of each chunk of u. If u.size() is greater than result.size() the extra elements of u are ignored.

Note
This is equivalent to treating the first part of u as a matrix and filling result (or the returned object) with the transpose of that matrix.

If u represents a block of data indexed by $(x, y, z, \ldots)$ with the first index varying fastest, transpose serves to rotate the indices. If the extents are $(X, Y, Z, \ldots)$, with product $N$, transpose(u, X, N/X) reorders the data to be indexed $(y, z, \ldots, x)$, transpose(u, X*Y, N/X/Y) reorders the data to be indexed $(z, \ldots, x, y)$, etc.

### Example

const DataVector matrix{ 1., 2., 3.,
4., 5., 6.,
7., 8., 9.,
10., 11., 12.};
CHECK(transpose(matrix, 3, 4) == DataVector{1., 4., 7., 10.,
2., 5., 8., 11.,
3., 6., 9., 12.});
const size_t chunk_size_vars = 8;
const size_t n_grid_pts = 2 * chunk_size_vars;
Variables<two_vars<2>> variables(n_grid_pts, 0.);
for (size_t i = 0; i < variables.size(); ++i) {
// clang-tidy: pointer arithmetic
variables.data()[i] = i * i; // NOLINT
}
const size_t number_of_chunks_vars = variables.size() / chunk_size_vars;
auto transposed_vars = variables;
transpose(make_not_null(&transposed_vars), variables, chunk_size_vars,
number_of_chunks_vars);
for (size_t i = 0; i < chunk_size_vars; ++i) {
for (size_t j = 0; j < number_of_chunks_vars; ++j) {
// clang-tidy: pointer arithmetic
CHECK(variables.data()[i + chunk_size_vars * j] == // NOLINT
transposed_vars.data()[j + number_of_chunks_vars * i]); // NOLINT
}
}
const size_t chunk_size = 8;
const size_t number_of_chunks = 2;
const size_t n_pts = chunk_size * number_of_chunks;
DataVector data(n_pts);
for (size_t i = 0; i < data.size(); ++i) {
data[i] = i * i;
}
DataVector transposed_data(n_pts, 0.);
transposed_data = transpose(data, chunk_size, number_of_chunks);
for (size_t i = 0; i < chunk_size; ++i) {
for (size_t j = 0; j < number_of_chunks; ++j) {
CHECK(data[i + chunk_size * j] ==
transposed_data[j + number_of_chunks * i]);
}
}
Variables<one_var<2>> partial_vars(n_grid_pts, 0.);
get<Var1<2>>(partial_vars) = get<Var1<2>>(variables);
Variables<one_var<2>> partial_transpose(n_grid_pts, 0.);
const size_t partial_number_of_chunks = 2*number_of_chunks_vars / 3;
transpose(make_not_null(&partial_transpose), variables, chunk_size_vars,
partial_number_of_chunks);
for (size_t i = 0; i < chunk_size_vars; ++i) {
for (size_t j = 0; j < partial_number_of_chunks; ++j) {
// clang-tidy: pointer arithmetic
CHECK(partial_transpose
.data()[j + partial_number_of_chunks * i] == // NOLINT
variables.data()[i + chunk_size_vars * j]); // NOLINT
CHECK(partial_transpose
.data()[j + partial_number_of_chunks * i] == // NOLINT
partial_vars.data()[i + chunk_size_vars * j]); // NOLINT
}
}
Template Parameters
 U the type of data to be transposed T the type of the transposed data