SpECTRE  v2024.06.18
fd Namespace Reference

Functions and classes for finite difference methods. More...

Namespaces

namespace  reconstruction
 Variable and flux vector splitting reconstruction schemes for finite difference methods.
 

Enumerations

enum class  DerivativeOrder : int {
  OneHigherThanRecons = -1 , OneHigherThanReconsButFiveToFour = -2 , Two = 2 , Four = 4 ,
  Six = 6 , Eight = 8 , Ten = 10
}
 Controls which FD derivative order is used. More...
 

Functions

std::ostreamoperator<< (std::ostream &os, DerivativeOrder der_order)
 
template<size_t Dim>
void low_pass_filter (gsl::not_null< gsl::span< double > * > filtered_data, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Mesh< Dim > &volume_mesh, size_t number_of_variables, size_t fd_order, double epsilon)
 Apply a low-pass filter to the data. More...
 
template<size_t Dim>
void kreiss_oliger_filter (gsl::not_null< gsl::span< double > * > filtered_data, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Mesh< Dim > &volume_mesh, size_t number_of_variables, size_t fd_order, double epsilon)
 Apply Kreiss-Oliger dissipation [111] of order fd_order to the variables. More...
 
template<size_t Dim, typename FluxesTags >
void set_cartesian_neighbor_cell_centered_fluxes (const gsl::not_null< DirectionMap< Dim, Variables< FluxesTags > > * > flux_neighbor_data, const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &all_ghost_data, const Mesh< Dim > &subcell_mesh, const size_t ghost_zone_size, const size_t number_of_rdmp_values_in_ghost_data)
 Fill the flux_neighbor_data with pointers into the all_ghost_data. More...
 
template<size_t Dim, typename... EvolvedVarsTags, typename FluxesTags = tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>
void cartesian_high_order_flux_corrections (const gsl::not_null< std::optional< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > > * > high_order_corrections, const std::optional< Variables< FluxesTags > > &cell_centered_fluxes, const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &second_order_boundary_corrections, const fd::DerivativeOrder &fd_derivative_order, const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &all_ghost_data, const Mesh< Dim > &subcell_mesh, const size_t ghost_zone_size, const std::array< gsl::span< std::uint8_t >, Dim > &reconstruction_order={}, const size_t number_of_rdmp_values_in_ghost_data=0)
 Computes the high-order Cartesian flux corrections if necessary. More...
 
template<size_t Dim, typename ReconstructionTags >
void neighbor_data_as_variables (const gsl::not_null< DirectionalIdMap< Dim, Variables< ReconstructionTags > > * > vars_neighbor_data, const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &all_ghost_data, const size_t ghost_zone_size, const Mesh< Dim > &subcell_mesh)
 Given the type-erased neighbor data for reconstruction stored in a DataVector, have Variables point into them. More...
 
template<size_t StencilSize>
std::array< std::array< double, StencilSize >, StencilSize > non_uniform_1d_weights (const std::deque< double > &times)
 Returns the weights for a 1D non-uniform finite difference stencil. More...
 
template<size_t Dim>
void logical_partial_derivatives (const gsl::not_null< std::array< gsl::span< double >, Dim > * > logical_derivatives, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Mesh< Dim > &volume_mesh, size_t number_of_variables, size_t fd_order)
 Compute the logical partial derivatives using cell-centered finite difference derivatives. 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 > > * > partial_derivatives, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double > > &ghost_cell_vars, const Mesh< Dim > &volume_mesh, size_t number_of_variables, size_t fd_order, const InverseJacobian< DataVector, Dim, Frame::ElementLogical, DerivativeFrame > &inverse_jacobian)
 Compute the partial derivative on the DerivativeFrame using the inverse_jacobian. More...
 
template<DerivativeOrder DerivOrder, size_t Dim, typename... EvolvedVarsTags>
void cartesian_high_order_fluxes_using_nodes (const gsl::not_null< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > * > high_order_boundary_corrections_in_logical_direction, const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &second_order_boundary_corrections_in_logical_direction, const Variables< tmpl::list< ::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > &cell_centered_inertial_flux, const DirectionMap< Dim, Variables< tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > > &ghost_cell_inertial_flux, const Mesh< Dim > &subcell_mesh, const size_t number_of_ghost_cells, const std::array< gsl::span< std::uint8_t >, Dim > &reconstruction_order={})
 Computes a high-order boundary correction \(G\) at the FD interface. More...
 
template<size_t Dim, typename... EvolvedVarsTags>
void cartesian_high_order_fluxes_using_nodes (const gsl::not_null< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > * > high_order_boundary_corrections_in_logical_direction, const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &second_order_boundary_corrections_in_logical_direction, const Variables< tmpl::list< ::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > &cell_centered_inertial_flux, const DirectionMap< Dim, Variables< tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > > &ghost_cell_inertial_flux, const Mesh< Dim > &subcell_mesh, const size_t number_of_ghost_cells, const DerivativeOrder derivative_order, const std::array< gsl::span< std::uint8_t >, Dim > &reconstruction_order={})
 Computes a high-order boundary correction \(G\) at the FD interface. More...
 

Detailed Description

Functions and classes for finite difference methods.

Enumeration Type Documentation

◆ DerivativeOrder

enum class fd::DerivativeOrder : int
strong

Controls which FD derivative order is used.

Enumerator
OneHigherThanRecons 

Use one order high derivative.

For example, if fifth order reconstruction is used, then a sixth-order derivative is used.

OneHigherThanReconsButFiveToFour 

Same as OneHigherThanRecons except uses a fourth-order derivative if fifth-order reconstruction was used.

Two 

Use 2nd order derivatives.

Four 

Use 4th order derivatives.

Six 

Use 6th order derivatives.

Eight 

Use 8th order derivatives.

Ten 

Use 10th order derivatives.

Function Documentation

◆ cartesian_high_order_flux_corrections()

template<size_t Dim, typename... EvolvedVarsTags, typename FluxesTags = tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>
void fd::cartesian_high_order_flux_corrections ( const gsl::not_null< std::optional< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > > * >  high_order_corrections,
const std::optional< Variables< FluxesTags > > &  cell_centered_fluxes,
const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &  second_order_boundary_corrections,
const fd::DerivativeOrder fd_derivative_order,
const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &  all_ghost_data,
const Mesh< Dim > &  subcell_mesh,
const size_t  ghost_zone_size,
const std::array< gsl::span< std::uint8_t >, Dim > &  reconstruction_order = {},
const size_t  number_of_rdmp_values_in_ghost_data = 0 
)

Computes the high-order Cartesian flux corrections if necessary.

The cell_centered_fluxes is stored in the tag evolution::dg::subcell::Tags::CellCenteredFlux, fd_derivative_order is from evolution::dg::subcell::Tags::SubcellOptions (.finite_difference_derivative_order()), the all_ghost_data is stored in the tag evolution::dg::subcell::Tags::GhostDataForReconstruction, the ghost_zone_size should come from the FD reconstructor.

By default we assume no RDMP data is in the ghost_data buffer. In the future we will want to update how we store the data in order to eliminate more memory allocations and copies, in which case that value will be non-zero.

Note
high_order_corrections must either not have a value or have all elements be of the same size as second_order_boundary_corrections[0].number_of_grid_points(), where we've assumed second_order_boundary_corrections is the same in all directions.

◆ cartesian_high_order_fluxes_using_nodes() [1/2]

template<size_t Dim, typename... EvolvedVarsTags>
void fd::cartesian_high_order_fluxes_using_nodes ( const gsl::not_null< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > * >  high_order_boundary_corrections_in_logical_direction,
const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &  second_order_boundary_corrections_in_logical_direction,
const Variables< tmpl::list< ::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > &  cell_centered_inertial_flux,
const DirectionMap< Dim, Variables< tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > > &  ghost_cell_inertial_flux,
const Mesh< Dim > &  subcell_mesh,
const size_t  number_of_ghost_cells,
const DerivativeOrder  derivative_order,
const std::array< gsl::span< std::uint8_t >, Dim > &  reconstruction_order = {} 
)

Computes a high-order boundary correction \(G\) at the FD interface.

The correction to the second-order boundary correction is given by

\begin{align*} G=G^{(2)}-G^{(4)}+G^{(6)}-G^{(8)}+G^{(10)}, \end{align*}

where

\begin{align*} G^{(4)}_{j+1/2}&=\frac{1}{6}\left(G_j -2 G^{(2)} + G_{j+1}\right), \\ G^{(6)}_{j+1/2}&=\frac{1}{180}\left(G_{j-1} - 9 G_j + 16 G^{(2)} -9 G_{j+1} + G_{j+2}\right), \\ G^{(8)}_{j+1/2}&=\frac{1}{2100}\left(G_{j-2} - \frac{25}{3} G_{j-1} + 50 G_j - \frac{256}{3} G^{(2)} + 50 G_{j+1} - \frac{25}{3} G_{j+2} +G_{j+3}\right), \\ G^{(10)}_{j+1/2}&=\frac{1}{17640} \left(G_{j-3} - \frac{49}{5} G_{j-2} + 49 G_{j-1} - 245 G_j + \frac{2048}{5} G^{(2)}\right. \nonumber \\ &\left.- 245 G_{j+1}+ 49 G_{j+2} - \frac{49}{5} G_{j+3} + G_{j+4}\right), \end{align*}

where

\begin{align*} G_{j} &= F^i_j n_i^{j+1/2}, \\ G_{j\pm1} &= F^i_{j\pm1} n_i^{j+1/2}, \\ G_{j\pm2} &= F^i_{j\pm2} n_i^{j+1/2}, \\ G_{j\pm3} &= F^i_{j\pm3} n_i^{j+1/2}, \\ G_{j\pm4} &= F^i_{j\pm4} n_i^{j+1/2}. \end{align*}

This is a generalization of the correction presented in [37].

This high-order flux can be fed into a flux limiter, e.g. to guarantee positivity.

Note
This implementation should be profiled and optimized.
Warning
This documentation is for the general case. In the restricted Cartesian case we use the cell-centered flux as opposed to G^{(4)}, which differs by a minus sign. This amounts to a minus sign change in front of the \(G^{(k)}\) terms in computing \(G\) for \(k>2\), and also a sign change in front of \(G^{(2)}\) in all \(G^{(k)}\) for \(k>2\).

◆ cartesian_high_order_fluxes_using_nodes() [2/2]

template<DerivativeOrder DerivOrder, size_t Dim, typename... EvolvedVarsTags>
void fd::cartesian_high_order_fluxes_using_nodes ( const gsl::not_null< std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > * >  high_order_boundary_corrections_in_logical_direction,
const std::array< Variables< tmpl::list< EvolvedVarsTags... > >, Dim > &  second_order_boundary_corrections_in_logical_direction,
const Variables< tmpl::list< ::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > &  cell_centered_inertial_flux,
const DirectionMap< Dim, Variables< tmpl::list<::Tags::Flux< EvolvedVarsTags, tmpl::size_t< Dim >, Frame::Inertial >... > > > &  ghost_cell_inertial_flux,
const Mesh< Dim > &  subcell_mesh,
const size_t  number_of_ghost_cells,
const std::array< gsl::span< std::uint8_t >, Dim > &  reconstruction_order = {} 
)

Computes a high-order boundary correction \(G\) at the FD interface.

The correction to the second-order boundary correction is given by

\begin{align*} G=G^{(2)}-G^{(4)}+G^{(6)}-G^{(8)}+G^{(10)}, \end{align*}

where

\begin{align*} G^{(4)}_{j+1/2}&=\frac{1}{6}\left(G_j -2 G^{(2)} + G_{j+1}\right), \\ G^{(6)}_{j+1/2}&=\frac{1}{180}\left(G_{j-1} - 9 G_j + 16 G^{(2)} -9 G_{j+1} + G_{j+2}\right), \\ G^{(8)}_{j+1/2}&=\frac{1}{2100}\left(G_{j-2} - \frac{25}{3} G_{j-1} + 50 G_j - \frac{256}{3} G^{(2)} + 50 G_{j+1} - \frac{25}{3} G_{j+2} +G_{j+3}\right), \\ G^{(10)}_{j+1/2}&=\frac{1}{17640} \left(G_{j-3} - \frac{49}{5} G_{j-2} + 49 G_{j-1} - 245 G_j + \frac{2048}{5} G^{(2)}\right. \nonumber \\ &\left.- 245 G_{j+1}+ 49 G_{j+2} - \frac{49}{5} G_{j+3} + G_{j+4}\right), \end{align*}

where

\begin{align*} G_{j} &= F^i_j n_i^{j+1/2}, \\ G_{j\pm1} &= F^i_{j\pm1} n_i^{j+1/2}, \\ G_{j\pm2} &= F^i_{j\pm2} n_i^{j+1/2}, \\ G_{j\pm3} &= F^i_{j\pm3} n_i^{j+1/2}, \\ G_{j\pm4} &= F^i_{j\pm4} n_i^{j+1/2}. \end{align*}

This is a generalization of the correction presented in [37].

This high-order flux can be fed into a flux limiter, e.g. to guarantee positivity.

Note
This implementation should be profiled and optimized.
Warning
This documentation is for the general case. In the restricted Cartesian case we use the cell-centered flux as opposed to G^{(4)}, which differs by a minus sign. This amounts to a minus sign change in front of the \(G^{(k)}\) terms in computing \(G\) for \(k>2\), and also a sign change in front of \(G^{(2)}\) in all \(G^{(k)}\) for \(k>2\).

◆ kreiss_oliger_filter()

template<size_t Dim>
void fd::kreiss_oliger_filter ( gsl::not_null< gsl::span< double > * >  filtered_data,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Mesh< Dim > &  volume_mesh,
size_t  number_of_variables,
size_t  fd_order,
double  epsilon 
)

Apply Kreiss-Oliger dissipation [111] of order fd_order to the variables.

Define the operators \(D_+\) and \(D_-\) be defined as:

\begin{align*} D_+f_i&=\frac{(f_{i+1}-f_i)}{\Delta x} \\ D_-f_i&=\frac{(f_i-f_{i-1})}{\Delta x} \end{align*}

where the subscript \(i\) refers to the grid index. The dissipation operators are generally applied dimension-by-dimension, and so we have:

\begin{align*} \mathcal{D}^{(2m)}=-\frac{(-1)^m}{2^{2m}}\Delta x^{2m-1} \epsilon(D_{+})^m(D_{-})^m \end{align*}

where \(\epsilon\) controls the amount of dissipation and is restricted to \(0\leq\epsilon\leq1\), and \(m\) is the order of the finite difference derivative so as not to spoil the accuracy of the scheme. That is, for second order FD, \(m=2\) and one should use \(\mathcal{D}^{(4)}\). However, this choice requires a larger stencil and whether or not this is necessary also depends on when and how in the algorithm the operators are applied.

We arrive at the following operators:

\begin{align*} \mathcal{D}^{(2)}f_i&=-\frac{\epsilon}{\Delta x} (f_{i+1} - 2f_i+f_{i-1}) \\ \mathcal{D}^{(4)}f_i&=-\frac{\epsilon}{16\Delta x} (f_{i+2}-4f_{i+1}+6f_i-4f_{i-1}+f_{i-2}) \\ \mathcal{D}^{(6)}f_i&=\frac{\epsilon}{64\Delta x}(f_{i+3}-6f_{i+2}+15f_{i+1} -20f_i+15f_{i-1}-6f_{i-2}+f_{i-3}) \\ \mathcal{D}^{(8)}f_i&=-\frac{\epsilon}{256\Delta x} (f_{i+4}-8f_{i+3}+28f_{i+2}-56f_{i+1}+70f_{i}-56f_{i-1}+28f_{i-2} -8f_{i-3}+f_{i-4}) \\ \mathcal{D}^{(10)}f_i&=\frac{\epsilon}{1024\Delta x} (f_{i+5}-10f_{i+4}+45f_{i+3} -120f_{i+2}+210f_{i+1}-252f_{i}+210f_{i-1}-120f_{i-2}+45f_{i-3} -10f_{i-4}+f_{i-5}) \end{align*}

Note
This function applies \(\Delta x \mathcal{D}^{(2m)}\).

◆ logical_partial_derivatives()

template<size_t Dim>
void fd::logical_partial_derivatives ( const gsl::not_null< std::array< gsl::span< double >, Dim > * >  logical_derivatives,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Mesh< Dim > &  volume_mesh,
size_t  number_of_variables,
size_t  fd_order 
)

Compute the logical partial derivatives using cell-centered finite difference derivatives.

Up to 8th order stencils are supported.

Note
Currently the stride is always one because we transpose the data before reconstruction. However, it may be faster to have a non-unit stride without the transpose. We have the stride parameter in the derivative stencils to make testing performance easier in the future.
This code does not do any explicit SIMD vectorization. We will want to profile and decide if optimization are possible. The Vc SIMD library has an example of vectorizing single-precision FD derivatives. There is also a paper "Optimization of Finite-Differencing Kernels for Numerical Relativity Applications" by Alfieri, Bernuzzi, Perego, and Radice that uses compiler auto-vectorization.

◆ low_pass_filter()

template<size_t Dim>
void fd::low_pass_filter ( gsl::not_null< gsl::span< double > * >  filtered_data,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Mesh< Dim > &  volume_mesh,
size_t  number_of_variables,
size_t  fd_order,
double  epsilon 
)

Apply a low-pass filter to the data.

The filter fits a Legendre polynomial of degree equal to fd_order to each grid point and its neighboring points, then subtracts out the highest mode contribution at the grid point. This is inspired by a Heaviside filter for Legendre polynomial-based spectral methods.

The filter at different orders is given by:

\begin{align*} F^{(9)} u_i&= -\frac{35 \times 531441}{128} \left( \frac{1}{6406400} (u_{i-4} + u_{i+4}) - \frac{1}{800800} (u_{i-3} + u_{i+3}) + \frac{1}{228800} (u_{i-2} + u_{i+2}) - \frac{1}{114400} (u_{i-1} + u_{i+1}) + \frac{1}{91520} u_i\right) \\ F^{(7)} u_i&= \frac{5 \times 16807}{16} \left( \frac{1}{95040} (u_{i-3} + u_{i+3}) - \frac{1}{15840} (u_{i-2} + u_{i+2}) + \frac{1}{6336} (u_{i-1} + u_{i+1}) - \frac{1}{4572} u_i\right) \\ F^{(5)} u_i&= -\frac{3 \times 125}{8} \left( \frac{1}{336} (u_{i-2} + u_{i+2}) - \frac{1}{84} (u_{i-1} + u_{i+1}) + \frac{1}{56} u_i\right) \\ F^{(3)} u_i&= \frac{1 \times 3}{2} \left( \frac{1}{4} (u_{i-1} + u_{i+1}) - \frac{1}{2} u_i\right) \end{align*}

Note
The \(F^{(11)}\) filter isn't implemented yet.
The argument \(\epsilon\) controls how much of the mode is filter out. \(\epsilon=1\) means the highest mode is completely filtered while \(\epsilon=0\) means it's not at all filtered.

◆ neighbor_data_as_variables()

template<size_t Dim, typename ReconstructionTags >
void fd::neighbor_data_as_variables ( const gsl::not_null< DirectionalIdMap< Dim, Variables< ReconstructionTags > > * >  vars_neighbor_data,
const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &  all_ghost_data,
const size_t  ghost_zone_size,
const Mesh< Dim > &  subcell_mesh 
)

Given the type-erased neighbor data for reconstruction stored in a DataVector, have Variables point into them.

This function is helpful for reconstruction, especially when wanting to apply different reconstruction methods to different tags. This can happen, for example, when doing positivity-preserving reconstruction. The density should remain positive, but negative velocities are fine.

◆ partial_derivatives()

template<typename DerivativeTags , size_t Dim, typename DerivativeFrame >
void fd::partial_derivatives ( gsl::not_null< Variables< db::wrap_tags_in< Tags::deriv, DerivativeTags, tmpl::size_t< Dim >, DerivativeFrame > > * >  partial_derivatives,
const gsl::span< const double > &  volume_vars,
const DirectionMap< Dim, gsl::span< const double > > &  ghost_cell_vars,
const Mesh< Dim > &  volume_mesh,
size_t  number_of_variables,
size_t  fd_order,
const InverseJacobian< DataVector, Dim, Frame::ElementLogical, DerivativeFrame > &  inverse_jacobian 
)

Compute the partial derivative on the DerivativeFrame using the inverse_jacobian.

Logical partial derivatives are first computed using the fd::logical_partial_derivatives() function.

◆ set_cartesian_neighbor_cell_centered_fluxes()

template<size_t Dim, typename FluxesTags >
void fd::set_cartesian_neighbor_cell_centered_fluxes ( const gsl::not_null< DirectionMap< Dim, Variables< FluxesTags > > * >  flux_neighbor_data,
const DirectionalIdMap< Dim, evolution::dg::subcell::GhostData > &  all_ghost_data,
const Mesh< Dim > &  subcell_mesh,
const size_t  ghost_zone_size,
const size_t  number_of_rdmp_values_in_ghost_data 
)

Fill the flux_neighbor_data with pointers into the all_ghost_data.

The all_ghost_data is stored in the tag evolution::dg::subcell::Tags::GhostDataForReconstruction, and the ghost_zone_size should come from the FD reconstructor.