|
SpECTRE
2021.08.02
|
Generic numerical algorithms. More...
Namespaces | |
| namespace | integration |
| Numerical integration algorithms. | |
| namespace | intrp |
| Contains classes and functions for interpolation. | |
| namespace | intrp::callbacks |
Contains callback functions called by InterpolationTargets. | |
| namespace | OdeIntegration |
| For ODE integration, we suggest using the boost libraries whenever possible. | |
Classes | |
| class | intrp::BarycentricRational |
| A barycentric rational interpolation class. More... | |
| class | intrp::CubicSpline |
| A natural cubic spline interpolation class. More... | |
| class | intrp::Irregular< Dim > |
Interpolates a Variables onto an arbitrary set of points. More... | |
| class | intrp::RegularGrid< Dim > |
| Interpolate data from a Mesh onto a regular grid of points. More... | |
Enumerations | |
| enum class | RootFinder::Method { Method::Hybrids , Method::Hybrid , Method::Newton } |
| The different options for the rootfinding method of gsl_multiroot. More... | |
| enum class | RootFinder::StoppingCondition { StoppingCondition::AbsoluteAndRelative , StoppingCondition::Absolute } |
| The different options for the convergence criterion of gsl_multiroot. More... | |
Functions | |
| 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... | |
| 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 function over a manifold. More... | |
| template<size_t Dim> | |
| double | mean_value (const DataVector &f, const Mesh< Dim > &mesh) noexcept |
| Compute the mean value of a function over a manifold. More... | |
| template<typename... FluxTags, size_t Dim> | |
| void | weak_divergence (const gsl::not_null< Variables< tmpl::list< Tags::div< FluxTags >... > > * > divergence_of_fluxes, const Variables< tmpl::list< FluxTags... > > &fluxes, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > &det_jac_times_inverse_jacobian) noexcept |
| Compute the weak form divergence of fluxes. 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) noexcept |
| Returns the positive root of a quadratic equation \(ax^2 + bx + c = 0\). More... | |
| template<typename T > | |
| T | smallest_root_greater_than_value_within_roundoff (const T &a, const T &b, const T &c, double value) noexcept |
| Returns the smallest root of a quadratic equation \(ax^2 + bx + c = 0\) that is greater than the given value, within roundoff. More... | |
| template<typename T > | |
| T | largest_root_between_values_within_roundoff (const T &a, const T &b, const T &c, double min_value, double max_value) noexcept |
| Returns the largest root of a quadratic equation \(ax^2 + bx + c = 0\) that is between min_value and max_value, within roundoff. More... | |
| std::array< double, 2 > | real_roots (double a, double b, double c) noexcept |
| Returns the two real roots of a quadratic equation \(ax^2 + bx + c = 0\) with the root closer to \(-\infty\) first. More... | |
| template<typename Functor > | |
| void | RootFinder::bracket_possibly_undefined_function_in_interval (const gsl::not_null< double * > lower_bound, const gsl::not_null< double * > upper_bound, const gsl::not_null< double * > f_at_lower_bound, const gsl::not_null< double * > f_at_upper_bound, const Functor &f, const double guess) noexcept |
Brackets the root of the function f, assuming a single root in a given interval \(f[x_\mathrm{lo},x_\mathrm{up}]\) and assuming that f is defined only in an unknown smaller interval \(f[x_a,x_b]\) where \(x_\mathrm{lo} \leq x_a \leq x_b \leq x_\mathrm{hi}\). More... | |
| template<typename Functor > | |
| void | RootFinder::bracket_possibly_undefined_function_in_interval (const gsl::not_null< DataVector * > lower_bound, const gsl::not_null< DataVector * > upper_bound, const gsl::not_null< DataVector * > f_at_lower_bound, const gsl::not_null< DataVector * > f_at_upper_bound, const Functor &f, const DataVector &guess) noexcept |
Brackets the single root of the function f for each element in a DataVector, assuming the root lies in the given interval and that f may be undefined at some points in the interval. More... | |
| template<typename Function > | |
| double | RootFinder::toms748 (const Function &f, const double lower_bound, const double upper_bound, const double f_at_lower_bound, const double f_at_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 > | |
| 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, where function values are not supplied at the lower and upper bounds. | |
| 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 Function > | |
| DataVector | RootFinder::toms748 (const Function &f, const DataVector &lower_bound, const DataVector &upper_bound, const DataVector &f_at_lower_bound, const DataVector &f_at_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, where function values are supplied at the lower and upper bounds. 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 ResultType , typename MatrixType , typename VectorType , size_t Dim> | |
| void | apply_matrices (const gsl::not_null< ResultType * > result, const std::array< MatrixType, Dim > &matrices, const VectorType &u, const Index< Dim > &extents) noexcept |
| Multiply by matrices in each dimension. More... | |
| template<typename MatrixType , typename VectorType , size_t Dim> | |
| VectorType | apply_matrices (const std::array< MatrixType, Dim > &matrices, const VectorType &u, const Index< Dim > &extents) noexcept |
| Multiply by matrices in each dimension. More... | |
| template<typename ResultType , typename MatrixType , typename VectorType , size_t Dim> | |
| ResultType | apply_matrices (const std::array< MatrixType, Dim > &matrices, const VectorType &u, const Index< Dim > &extents) noexcept |
| Multiply by matrices in each dimension. 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> | |
| T | transpose (const U &u, const size_t chunk_size, const size_t number_of_chunks) noexcept |
| Function to compute transposed data. 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<typename... DivTags, typename... FluxTags, size_t Dim, typename DerivativeFrame > | |
| void | divergence (gsl::not_null< Variables< tmpl::list< DivTags... > > * > divergence_of_F, const Variables< tmpl::list< FluxTags... > > &F, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept |
| Compute the (Euclidean) divergence of fluxes. | |
| template<size_t Dim, typename DerivativeFrame > | |
| Scalar< DataVector > | divergence (const tnsr::I< DataVector, Dim, DerivativeFrame > &input, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept |
Compute the divergence of the vector input | |
| template<size_t Dim, typename DerivativeFrame > | |
| void | divergence (gsl::not_null< Scalar< DataVector > * > div_input, const tnsr::I< DataVector, Dim, DerivativeFrame > &input, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept |
Compute the divergence of the vector input | |
| template<size_t Dim, typename VectorType > | |
| void | indefinite_integral (gsl::not_null< VectorType * > integral, const VectorType &integrand, const Mesh< Dim > &mesh, size_t dim_to_integrate) noexcept |
Compute the indefinite integral of a function in the dim_to_integrate, applying a zero boundary condition on each stripe. More... | |
| template<size_t Dim, typename VectorType > | |
| VectorType | indefinite_integral (const VectorType &integrand, const Mesh< Dim > &mesh, size_t dim_to_integrate) noexcept |
Compute the indefinite integral of a function in the dim_to_integrate, applying a zero boundary condition on each stripe. 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 |
| Compute the mean value of a function over a boundary of a manifold. More... | |
| 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 function over a boundary of a manifold. More... | |
| 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 function over a boundary of a manifold. More... | |
| 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 function over a boundary of a manifold. More... | |
| 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 function over a boundary of a manifold. More... | |
| 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 SymmList , typename IndexList , size_t Dim> | |
| void | logical_partial_derivative (gsl::not_null< TensorMetafunctions::prepend_spatial_index< Tensor< DataVector, SymmList, IndexList >, Dim, UpLo::Lo, Frame::Logical > * > logical_derivative_of_u, gsl::not_null< gsl::span< double > * > buffer, const Tensor< DataVector, SymmList, IndexList > &u, const Mesh< Dim > &mesh) noexcept |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b). More... | |
| template<typename SymmList , typename IndexList , size_t Dim> | |
| void | logical_partial_derivative (gsl::not_null< TensorMetafunctions::prepend_spatial_index< Tensor< DataVector, SymmList, IndexList >, Dim, UpLo::Lo, Frame::Logical > * > logical_derivative_of_u, const Tensor< DataVector, SymmList, IndexList > &u, const Mesh< Dim > &mesh) noexcept |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b). More... | |
| template<typename SymmList , typename IndexList , size_t Dim> | |
| auto | logical_partial_derivative (const Tensor< DataVector, SymmList, IndexList > &u, const Mesh< Dim > &mesh) noexcept -> TensorMetafunctions::prepend_spatial_index< Tensor< DataVector, SymmList, IndexList >, Dim, UpLo::Lo, Frame::Logical > |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b). 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<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... | |
Generic numerical algorithms.
|
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.
|
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. |
|
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.
double or std::complex<double>. The matrix, however, must be real. In the case of acting on a vector of complex values, the matrix is treated as having zero imaginary part. This is chosen for efficiency in all use-cases for spectral matrix arithmetic so far encountered.
|
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.
double or std::complex<double>. The matrix, however, must be real. In the case of acting on a vector of complex values, the matrix is treated as having zero imaginary part. This is chosen for efficiency in all use-cases for spectral matrix arithmetic so far encountered.
|
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.
double or std::complex<double>. The matrix, however, must be real. In the case of acting on a vector of complex values, the matrix is treated as having zero imaginary part. This is chosen for efficiency in all use-cases for spectral matrix arithmetic so far encountered.
|
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.
double or std::complex<double>. The matrix, however, must be real. In the case of acting on a vector of complex values, the matrix is treated as having zero imaginary part. This is chosen for efficiency in all use-cases for spectral matrix arithmetic so far encountered.
|
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.
double or std::complex<double>. The matrix, however, must be real. In the case of acting on a vector of complex values, the matrix is treated as having zero imaginary part. This is chosen for efficiency in all use-cases for spectral matrix arithmetic so far encountered.
|
noexcept |
Brackets the single root of the function f for each element in a DataVector, assuming the root lies in the given interval and that f may be undefined at some points in the interval.
f is a binary invokable that takes a double and a size_t as arguments. The double is the current value at which to evaluate f, and the size_t is the index into the DataVectors. f returns a std::optional<double> which evaluates to false if the function is undefined at the supplied point.
Assumes that there is only one root in the interval.
Assumes that if \(f(x_1)\) and \(f(x_2)\) are both defined for some \((x_1,x_2)\), then \(f(x)\) is defined for all \(x\) between \(x_1\) and \(x_2\).
On input, assumes that the root lies in the interval [lower_bound,upper_bound]. Optionally takes a guess for the location of the root.
On return, lower_bound and upper_bound are replaced with values that bracket the root and for which the function is defined, and f_at_lower_bound and f_at_upper_bound are replaced with f evaluated at those bracketing points.
|
noexcept |
Brackets the root of the function f, assuming a single root in a given interval \(f[x_\mathrm{lo},x_\mathrm{up}]\) and assuming that f is defined only in an unknown smaller interval \(f[x_a,x_b]\) where \(x_\mathrm{lo} \leq x_a \leq x_b \leq x_\mathrm{hi}\).
f is a unary invokable that takes a double which is the current value at which to evaluate f. f returns a std::optional<double> which evaluates to false if the function is undefined at the supplied point.
Assumes that there is only one root in the interval.
Assumes that if \(f(x_1)\) and \(f(x_2)\) are both defined for some \((x_1,x_2)\), then \(f(x)\) is defined for all \(x\) between \(x_1\) and \(x_2\).
On input, assumes that the root lies in the interval [lower_bound,upper_bound]. Optionally takes a guess for the location of the root. If guess is supplied, then evaluates the function first at guess and upper_bound before trying lower_bound: this means that it would be optimal if guess underestimates the actual root and if upper_bound was less likely to be undefined than lower_bound.
On return, lower_bound and upper_bound are replaced with values that bracket the root and for which the function is defined, and f_at_lower_bound and f_at_upper_bound are replaced with f evaluated at those bracketing points.
bracket_possibly_undefined_function_in_interval throws an error if all points are valid but of the same sign (because that would indicate multiple roots but we assume only one root), if no root exists, or if the range of a sign change is sufficently small relative to the given interval that the number of iterations to find the root is exceeded.
|
noexcept |
Compute the definite integral of a function over a manifold.
Given a function \(f\), compute its integral \(I\) with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\). E.g., in 1 dimension, \(I = \int_{-1}^1 f d\xi\).
The integral w.r.t. a different set of coordinates \(\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi})\) can be computed by pre-multiplying \(f\) by the Jacobian determinant \(J = \det d\boldsymbol{x}/d\boldsymbol{\xi}\) of the mapping \(\boldsymbol{x}(\boldsymbol{\xi})\). Note that, in the \(\boldsymbol{x}\) coordinates, the domain of integration is the image of the logical cube (square in 2D, interval in 1D) under the mapping.
The integral is computed by quadrature, using the quadrature rule for the basis associated with the collocation points.
| integrand | the function to integrate. |
| mesh | the Mesh defining the grid points on the manifold. |
|
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).
| 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.
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.
| 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. |
|
noexcept |
Compute the indefinite integral of a function in the dim_to_integrate, applying a zero boundary condition on each stripe.
Integrates with respect to one of the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
The integral w.r.t. a different set of coordinates \(\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi})\) can be computed by pre-multiplying integrand by the Jacobian determinant \(J = \det d\boldsymbol{x}/d\boldsymbol{\xi}\) of the mapping \(\boldsymbol{x}(\boldsymbol{\xi})\). The integration is still performed along one logical-coordinate direction, indicated by dim_to_integrate.
Requires: number of points in integrand and mesh are equal.
|
noexcept |
Compute the indefinite integral of a function in the dim_to_integrate, applying a zero boundary condition on each stripe.
Integrates with respect to one of the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
The integral w.r.t. a different set of coordinates \(\boldsymbol{x} = \boldsymbol{x}(\boldsymbol{\xi})\) can be computed by pre-multiplying integrand by the Jacobian determinant \(J = \det d\boldsymbol{x}/d\boldsymbol{\xi}\) of the mapping \(\boldsymbol{x}(\boldsymbol{\xi})\). The integration is still performed along one logical-coordinate direction, indicated by dim_to_integrate.
Requires: number of points in integrand and mesh are equal.
|
noexcept |
Returns the largest root of a quadratic equation \(ax^2 + bx + c = 0\) that is between min_value and max_value, within roundoff.
Returns: A root of a quadratic equation.
Requires: That there are two real roots.
Requires: At least one root is between min_value and max_value, to roundoff.
|
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\)
gsl::not_null variant assumes *result is of the correct size.
|
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.gsl::not_null variant assumes *result is of the correct size.
|
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\)
gsl::not_null variant assumes *result is of the correct size.
|
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.gsl::not_null variant assumes *result is of the correct size.
|
noexcept |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b).
There is an overload that accepts a buffer of size mesh.number_of_grid_points() or larger. When passed this function performs no memory allocations, which helps improve performance.
If you have a Variables with several tensors you need to differentiate you should use the logical_partial_derivatives function that operates on Variables since that'll be more efficient.
|
noexcept |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b).
There is an overload that accepts a buffer of size mesh.number_of_grid_points() or larger. When passed this function performs no memory allocations, which helps improve performance.
If you have a Variables with several tensors you need to differentiate you should use the logical_partial_derivatives function that operates on Variables since that'll be more efficient.
|
noexcept |
Computes the logical partial derivative of a tensor, prepending the spatial derivative index, e.g. for \(\partial_i T_{a}{}^{b}\) the C++ call is get(i, a, b).
There is an overload that accepts a buffer of size mesh.number_of_grid_points() or larger. When passed this function performs no memory allocations, which helps improve performance.
If you have a Variables with several tensors you need to differentiate you should use the logical_partial_derivatives function that operates on Variables since that'll be more efficient.
|
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.
| DerivativeTags | the subset of VariableTags for which derivatives are computed. |
|
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.
| DerivativeTags | the subset of VariableTags for which derivatives are computed. |
|
noexcept |
Compute the mean value of a function over a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\) with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\). E.g., in 1 dimension, \(\bar{f} = \int_{-1}^1 f d\xi \Big/ \int_{-1}^1 d\xi\).
mean_value, \(f\) must be pre-multiplied by the Jacobian determinant \(J = \det d\boldsymbol{x}/d\boldsymbol{\xi}\) of the mapping \(\boldsymbol{x}(\boldsymbol{\xi})\). Additionally, the output of mean_value must be multiplied by a factor \(2^{\text{d}} / \int J d^{\text{d}}\xi\) (in \(d\) dimensions), to account for the different volume of the manifold in the \(\boldsymbol{x}\) coordinates.| f | the function to average. |
| mesh | the Mesh defining the grid points on the manifold. |
|
noexcept |
Compute the mean value of a function over a boundary of a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\), over a boundary, with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
mean_value for notes about means w.r.t. other coordinates.f the function to average.mesh the Mesh defining the grid points on the manifold.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.
|
noexcept |
Compute the mean value of a function over a boundary of a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\), over a boundary, with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
mean_value for notes about means w.r.t. other coordinates.f the function to average.mesh the Mesh defining the grid points on the manifold.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.
|
noexcept |
Compute the mean value of a function over a boundary of a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\), over a boundary, with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
mean_value for notes about means w.r.t. other coordinates.f the function to average.mesh the Mesh defining the grid points on the manifold.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.
|
noexcept |
Compute the mean value of a function over a boundary of a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\), over a boundary, with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
mean_value for notes about means w.r.t. other coordinates.f the function to average.mesh the Mesh defining the grid points on the manifold.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.
|
noexcept |
Compute the mean value of a function over a boundary of a manifold.
Given a function \(f\), compute its mean value \(\bar{f}\), over a boundary, with respect to the logical coordinates \(\boldsymbol{\xi} = (\xi, \eta, \zeta)\).
mean_value for notes about means w.r.t. other coordinates.f the function to average.mesh the Mesh defining the grid points on the manifold.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. | 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. f must return a std::pair<double, double> where the first element is the function value and the second element is the derivative of the function. 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.
See the Boost documentation for more details.
Requires: Function f be callable with a double and a size_t
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.
| 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. f must return a std::pair<double, double> where the first element is the function value and the second element is the derivative of the function. An example is below.
See the Boost documentation for more details.
Requires: Function f is invokable with a double
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.
|
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.
| DerivativeTags | the subset of VariableTags for which derivatives are computed. |
|
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.
| DerivativeTags | the subset of VariableTags for which derivatives are computed. |
|
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.
| DerivativeTags | the subset of VariableTags for which derivatives are computed. |
|
noexcept |
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.
|
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.
|
noexcept |
Returns the two real roots of a quadratic equation \(ax^2 + bx + c = 0\) with the root closer to \(-\infty\) first.
Returns: An array of the roots of a quadratic equation
Requires: That there are two real roots.
|
noexcept |
Returns the smallest root of a quadratic equation \(ax^2 + bx + c = 0\) that is greater than the given value, within roundoff.
Returns: A root of a quadratic equation.
Requires: That there are two real roots.
Requires: At least one root is greater than the given value, to roundoff.
| DataVector RootFinder::toms748 | ( | const Function & | f, |
| const DataVector & | lower_bound, | ||
| const DataVector & | upper_bound, | ||
| const DataVector & | f_at_lower_bound, | ||
| const DataVector & | f_at_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, where function values are supplied at the lower and upper bounds.
Supplying function values is an optimization that saves two function calls per point. The function values are often available because one often checks if the root is bracketed before calling toms748.
| 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.
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.
| double RootFinder::toms748 | ( | const Function & | f, |
| const double | lower_bound, | ||
| const double | upper_bound, | ||
| const double | f_at_lower_bound, | ||
| const double | f_at_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.
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.
The arguments f_at_lower_bound and f_at_upper_bound are optional, and are the function values at lower_bound and upper_bound. These function values are often known because the user typically checks if a root is bracketed before calling toms748; passing the function values here saves two function evaluations.
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.
|
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.
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.
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.
| U | the type of data to be transposed |
| T | the type of the transposed data |
|
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.
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.
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.
| U | the type of data to be transposed |
| T | the type of the transposed data |
|
noexcept |
Compute the weak form divergence of fluxes.
In a discontinuous Galerkin scheme we integrate the equations against the basis functions over the element. For the flux divergence term this gives:
\begin{align*} \int_{\Omega}d^n x \phi_{\breve{\imath}}\partial_i F^i, \end{align*}
where the basis functions are denoted by \(\phi_{\breve{\imath}}\).
Integrating by parts we get
\begin{align*} \int_{\Omega}d^n x\, \phi_{\breve{\imath}}\partial_i F^i = -\int_{\Omega}d^n x\, F^i \partial_i \phi_{\breve{\imath}} + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i \phi_{\breve{\imath}} \end{align*}
Next we expand the flux \(F^i\) in terms of the basis functions, yielding
\begin{align*} - \int_{\Omega}d^n x\,F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i \phi_{\breve{\imath}} + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \phi_{\breve{\imath}} \end{align*}
This function computes the volume term:
\begin{align*} \int_{\Omega}d^n x\,F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i \phi_{\breve{\imath}} \end{align*}