PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces > Class Template Reference

Sum Factorization derived class. More...

#include <operators.h>

Inheritance diagram for PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >:

Collaboration diagram for PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >:

Public Member Functions
	SumFactorizedOperators (const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
	Precompute 1D operator in constructor.
template<typename real >
void	matrix_vector_mult (const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const bool adding=false, const double factor=1.0)
	Computes a matrix-vector product using sum-factorization. Pass the one-dimensional basis, where x runs the fastest, then y, and z runs the slowest. Also, assume each one-dimensional basis is the same size. More...
template<typename real >
void	divergence_matrix_vector_mult (const dealii::Tensor< 1, dim, std::vector< real >> &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const dealii::FullMatrix< double > &gradient_basis_x, const dealii::FullMatrix< double > &gradient_basis_y, const dealii::FullMatrix< double > &gradient_basis_z)
	Computes the divergence using the sum factorization matrix-vector multiplication. More...
template<typename real >
void	divergence_matrix_vector_mult_1D (const dealii::Tensor< 1, dim, std::vector< real >> &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis, const dealii::FullMatrix< double > &gradient_basis)
	Computes the divergence using sum-factorization where the basis are the same in each direction.
template<typename real >
void	gradient_matrix_vector_mult (const std::vector< real > &input_vect, dealii::Tensor< 1, dim, std::vector< real >> &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const dealii::FullMatrix< double > &gradient_basis_x, const dealii::FullMatrix< double > &gradient_basis_y, const dealii::FullMatrix< double > &gradient_basis_z)
	Computes the gradient of a scalar using sum-factorization.
template<typename real >
void	gradient_matrix_vector_mult_1D (const std::vector< real > &input_vect, dealii::Tensor< 1, dim, std::vector< real >> &output_vect, const dealii::FullMatrix< double > &basis, const dealii::FullMatrix< double > &gradient_basis)
	Computes the gradient of a scalar using sum-factorization where the basis are the same in each direction.
template<typename real >
void	inner_product (const std::vector< real > &input_vect, const std::vector< double > &weight_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const bool adding=false, const double factor=1.0)
	Computes the inner product between a matrix and a vector multiplied by some weight function. More...
void	divergence_two_pt_flux_Hadamard_product (const dealii::Tensor< 1, dim, dealii::FullMatrix< double >> &input_mat, std::vector< double > &output_vect, const std::vector< double > &weights, const dealii::FullMatrix< double > &basis, const double scaling=2.0)
	Computes the divergence of the 2pt flux Hadamard products, then sums the rows. More...
void	surface_two_pt_flux_Hadamard_product (const dealii::FullMatrix< double > &input_mat, std::vector< double > &output_vect_vol, std::vector< double > &output_vect_surf, const std::vector< double > &weights, const std::array< dealii::FullMatrix< double >, 2 > &surf_basis, const unsigned int iface, const unsigned int dim_not_zero, const double scaling=2.0)
	Computes the surface cross Hadamard products for skew-symmetric form from Eq. (15) in Chan, Jesse. "Skew-symmetric entropy stable modal discontinuous Galerkin formulations." Journal of Scientific Computing 81.1 (2019): 459-485.
void	two_pt_flux_Hadamard_product (const dealii::FullMatrix< double > &input_mat, dealii::FullMatrix< double > &output_mat, const dealii::FullMatrix< double > &basis, const std::vector< double > &weights, const int direction)
	Computes the Hadamard product ONLY for 2pt flux calculations. More...
void	sum_factorized_Hadamard_sparsity_pattern (const unsigned int rows_size, const unsigned int columns_size, std::vector< std::array< unsigned int, dim >> &rows, std::vector< std::array< unsigned int, dim >> &columns)
	Computes the rows and columns vectors with non-zero indices for sum-factorized Hadamard products.
void	sum_factorized_Hadamard_basis_assembly (const unsigned int rows_size_1D, const unsigned int columns_size_1D, const std::vector< std::array< unsigned int, dim >> &rows, const std::vector< std::array< unsigned int, dim >> &columns, const dealii::FullMatrix< double > &basis, const std::vector< double > &weights, std::array< dealii::FullMatrix< double >, dim > &basis_sparse)
	Constructs the \( n^d \times n\) basis operator storing all non-zero entries for a "sum-factorized" Hadamard product.
void	sum_factorized_Hadamard_surface_sparsity_pattern (const unsigned int rows_size, const unsigned int columns_size, std::vector< unsigned int > &rows, std::vector< unsigned int > &columns, const int dim_not_zero)
	Computes the rows and columns vectors with non-zero indices for surface sum-factorized Hadamard products.
void	sum_factorized_Hadamard_surface_basis_assembly (const unsigned int rows_size, const unsigned int columns_size_1D, const std::vector< unsigned int > &rows, const std::vector< unsigned int > &columns, const dealii::FullMatrix< double > &basis, const std::vector< double > &weights, dealii::FullMatrix< double > &basis_sparse, const int dim_not_zero)
	Constructs the \( n^{d-1} \times n\) basis operator storing all non-zero entries for a "sum-factorized" surface Hadamard product.
template<typename real >
void	matrix_vector_mult_1D (const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const bool adding=false, const double factor=1.0)
	Apply the matrix vector operation using the 1D operator in each direction. More...
template<typename real >
void	inner_product_1D (const std::vector< real > &input_vect, const std::vector< double > &weight_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const bool adding=false, const double factor=1.0)
	Apply the inner product operation using the 1D operator in each direction.
template<typename real >
void	matrix_vector_mult_surface_1D (const std::vector< bool > face_orientation, const unsigned int face_number, const std::vector< real > &input_vect, std::vector< real > &output_vect, const std::array< dealii::FullMatrix< double >, 2 > &basis_surf, const dealii::FullMatrix< double > &basis_vol, const bool adding=false, const double factor=1.0)
	Apply sum-factorization matrix vector multiplication on a surface. More...
template<typename real >
void	inner_product_surface_1D (const std::vector< bool > face_orientation, const unsigned int face_number, const std::vector< real > &input_vect, const std::vector< double > &weight_vect, std::vector< real > &output_vect, const std::array< dealii::FullMatrix< double >, 2 > &basis_surf, const dealii::FullMatrix< double > &basis_vol, const bool adding=false, const double factor=1.0)
	Apply sum-factorization inner product on a surface.
template<typename real >
void	face_orientation_tensor_product (const std::vector< bool > face_orientation, const unsigned int face_number, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis)
	These function correct the face orientation of a face, in case deal.ii changes it. More...
template<typename real >
void	face_orientation_inner_product (const std::vector< bool > face_orientation, const unsigned int face_number, const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis)
void	Hadamard_product (const dealii::FullMatrix< double > &input_mat1, const dealii::FullMatrix< double > &input_mat2, dealii::FullMatrix< double > &output_mat)
	Computes a single Hadamard product. More...
template<typename real >
void	Hadamard_product_AD_vector (const dealii::FullMatrix< double > &input_mat1, const std::vector< real > &input_mat2, std::vector< real > &output_mat)
	Computes a single Hadamard product for AD type. More...
Public Member Functions inherited from PHiLiP::OPERATOR::OperatorsBase< dim, n_faces >
virtual	~OperatorsBase ()=default
	Destructor.
	OperatorsBase (const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
	Constructor.
dealii::FullMatrix< double >	tensor_product (const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z)
	Returns the tensor product of matrices passed.
dealii::FullMatrix< double >	tensor_product_state (const int nstate, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z)
	Returns the tensor product of matrices passed, but makes it sparse diagonal by state. More...
double	compute_factorial (double n)
	Standard function to compute factorial of a number.

Public Attributes
dealii::FullMatrix< double >	oneD_vol_operator
	Stores the one dimensional volume operator.
std::array< dealii::FullMatrix< double >, 2 >	oneD_surf_operator
	Stores the one dimensional surface operator. More...
dealii::FullMatrix< double >	oneD_grad_operator
	Stores the one dimensional gradient operator.
std::array< dealii::FullMatrix< double >, 2 >	oneD_surf_grad_operator
	Stores the one dimensional surface gradient operator.
Public Attributes inherited from PHiLiP::OPERATOR::OperatorsBase< dim, n_faces >
const unsigned int	max_degree
	Max polynomial degree.
const unsigned int	max_grid_degree
	Max grid degree.
const int	nstate
	Number of states.

Additional Inherited Members
Protected Attributes inherited from PHiLiP::OPERATOR::OperatorsBase< dim, n_faces >
unsigned int	max_grid_degree_check
	Check to see if the metrics used are a higher order then the initialized grid.
const MPI_Comm	mpi_communicator
	MPI communicator.
dealii::ConditionalOStream	pcout
	Parallel std::cout that only outputs on mpi_rank==0.

Detailed Description

template<int dim, int n_faces>
class PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >

Sum Factorization derived class.

Definition at line 112 of file operators.h.

Member Function Documentation

divergence_matrix_vector_mult()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::divergence_matrix_vector_mult< RadFadType >	(	const dealii::Tensor< 1, dim, std::vector< real >> &	input_vect,
		std::vector< real > &	output_vect,
		const dealii::FullMatrix< double > &	basis_x,
		const dealii::FullMatrix< double > &	basis_y,
		const dealii::FullMatrix< double > &	basis_z,
		const dealii::FullMatrix< double > &	gradient_basis_x,
		const dealii::FullMatrix< double > &	gradient_basis_y,
		const dealii::FullMatrix< double > &	gradient_basis_z
	)

Computes the divergence using the sum factorization matrix-vector multiplication.

Often, we compute a dot product in dim, where each matrix multiplictaion uses sum factorization. Example, consider taking the reference divergence of the reference flux:

\[ \nabla^r \cdot f^r = \frac{\partial f^r}{\partial \xi} + \frac{\partial f^r}{\partial \eta} + \frac{\partial f^r}{\partial \zeta} = \left( \frac{d \mathbf{\chi}(\mathbf{\xi})}{d\xi} \otimes \mathbf{\chi}(\mathbf{\eta}) \otimes \mathbf{\chi}(\mathbf{\zeta}) \right) \left(\hat{\mathbf{f}^r}\right)^T + \left( \mathbf{\chi}(\mathbf{\xi}) \otimes \frac{d\mathbf{\chi}(\mathbf{\eta})}{d\eta} \otimes \mathbf{\chi}(\mathbf{\zeta}) \right) \left(\hat{\mathbf{f}^r}\right)^T + \left( \mathbf{\chi}(\mathbf{\xi}) \otimes \mathbf{\chi}(\mathbf{\eta} \otimes \frac{d\mathbf{\chi}(\mathbf{\zeta})}{d\zeta})\right) \left(\hat{\mathbf{f}^r}\right)^T, \]

where we use sum factorization to evaluate each matrix-vector multiplication in each dim direction.

Definition at line 494 of file operators.cpp.

divergence_two_pt_flux_Hadamard_product()

template<int dim, int n_faces>

void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::divergence_two_pt_flux_Hadamard_product	(	const dealii::Tensor< 1, dim, dealii::FullMatrix< double >> &	input_mat,
		std::vector< double > &	output_vect,
		const std::vector< double > &	weights,
		const dealii::FullMatrix< double > &	basis,
		const double	scaling = 2.0
	)

Computes the divergence of the 2pt flux Hadamard products, then sums the rows.

Note that we have the factor of 2.0 definied in this function. We also make use of the structure of the flux basis to get the matrix vector product after the Hadamard product to be \( \mathcal{O}(n^{d+1})\).

Definition at line 652 of file operators.cpp.

face_orientation_tensor_product()

template<int dim, int n_faces>

template<typename real >

void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::face_orientation_tensor_product	(	const std::vector< bool >	face_orientation,
		const unsigned int	face_number,
		std::vector< real > &	output_vect,
		const dealii::FullMatrix< double > &	basis
	)

These function correct the face orientation of a face, in case deal.ii changes it.

When we have to interpolate to the face, the reference quadrature points that we interpolate to may have changed orientation. The functions below account for this and change the ordering of the quadrature points. After we compute the fluxes (i.e., when the points need to match between the int and ext cells), we need to account for this again when we perform inner products.

Definition at line 209 of file operators.cpp.

Hadamard_product()

template<int dim, int n_faces>

void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::Hadamard_product	(	const dealii::FullMatrix< double > &	input_mat1,
		const dealii::FullMatrix< double > &	input_mat2,
		dealii::FullMatrix< double > &	output_mat
	)

Computes a single Hadamard product.

For input mat1 \( A \) and input mat2 \( B \), this computes \( A \circ B = C \implies \left( C \right)_{ij} = \left( A \right)_{ij}\left( B \right)_{ij}\).

Definition at line 911 of file operators.cpp.

Hadamard_product_AD_vector()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::Hadamard_product_AD_vector< RadFadType >	(	const dealii::FullMatrix< double > &	input_mat1,
		const std::vector< real > &	input_mat2,
		std::vector< real > &	output_mat
	)

Computes a single Hadamard product for AD type.

For input mat1 \( A \) and input mat2 \( B \), this computes \( A \circ B = C \implies \left( C \right)_{ij} = \left( A \right)_{ij}\left( B \right)_{ij}\). The input_mat1 is a matrix of double, input_mat2 is a matrix stored as a vector in contiguous memory, and outputs a matrix stored as a vector in contiguous memory.

Definition at line 931 of file operators.cpp.

inner_product()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::inner_product< RadFadType >	(	const std::vector< real > &	input_vect,
		const std::vector< double > &	weight_vect,
		std::vector< real > &	output_vect,
		const dealii::FullMatrix< double > &	basis_x,
		const dealii::FullMatrix< double > &	basis_y,
		const dealii::FullMatrix< double > &	basis_z,
		const bool	adding = false,
		const double	factor = 1.0
	)

Computes the inner product between a matrix and a vector multiplied by some weight function.

That is, we compute \( \int Awu d\mathbf{\Omega}_r = \mathbf{A}^T \text{diag}(w) \mathbf{u}^T \). When using this function, pass \( \mathbf{A} \) and NOT it's transpose–the function transposes it in the first few lines.

Definition at line 575 of file operators.cpp.

matrix_vector_mult()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::matrix_vector_mult< RadFadType >	(	const std::vector< real > &	input_vect,
		std::vector< real > &	output_vect,
		const dealii::FullMatrix< double > &	basis_x,
		const dealii::FullMatrix< double > &	basis_y,
		const dealii::FullMatrix< double > &	basis_z,
		const bool	adding = false,
		const double	factor = 1.0
	)

Computes a matrix-vector product using sum-factorization. Pass the one-dimensional basis, where x runs the fastest, then y, and z runs the slowest. Also, assume each one-dimensional basis is the same size.

Uses sum-factorization with BLAS techniques to solve the the matrix-vector multiplication, where the matrix is the tensor product of three one-dimensional matrices. We use the standard notation that x runs the fastest, then y, and z runs the slowest. For an operator \(\mathbf{A}\) of size \(n^d\) with \(n\) the one dimensional dofs and \( d\) the dim, for ouput row vector \( v \) and input row vector \( u \), we compute \(v^T=\mathbf{A}u^T\). Lastly, the adding allows the result to add onto the previous output_vect scaled by "factor".

Definition at line 294 of file operators.cpp.

matrix_vector_mult_1D()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::matrix_vector_mult_1D< RadFadType >	(	const std::vector< real > &	input_vect,
		std::vector< real > &	output_vect,
		const dealii::FullMatrix< double > &	basis_x,
		const bool	adding = false,
		const double	factor = 1.0
	)

Apply the matrix vector operation using the 1D operator in each direction.

This is for the case where the operator of size dim is the dyadic product of the same 1D operator in each direction

Definition at line 402 of file operators.cpp.

matrix_vector_mult_surface_1D()

template<int dim, int n_faces>

template<typename real >

template void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::matrix_vector_mult_surface_1D< RadFadType >	(	const std::vector< bool >	face_orientation,
		const unsigned int	face_number,
		const std::vector< real > &	input_vect,
		std::vector< real > &	output_vect,
		const std::array< dealii::FullMatrix< double >, 2 > &	basis_surf,
		const dealii::FullMatrix< double > &	basis_vol,
		const bool	adding = false,
		const double	factor = 1.0
	)

Apply sum-factorization matrix vector multiplication on a surface.

Often times we have to interpolate to a surface, where in multiple dimensions, that's the tensor product of a surface operator with volume operators. This simplifies the function call. Explicitly, this passes basis_surf in the direction by face_number, and basis_vol in all other directions.

Definition at line 414 of file operators.cpp.

two_pt_flux_Hadamard_product()

template<int dim, int n_faces>

void PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::two_pt_flux_Hadamard_product	(	const dealii::FullMatrix< double > &	input_mat,
		dealii::FullMatrix< double > &	output_mat,
		const dealii::FullMatrix< double > &	basis,
		const std::vector< double > &	weights,
		const int	direction
	)

Computes the Hadamard product ONLY for 2pt flux calculations.

The Hadamard product comes up naturally in calculcating 2-point fluxes, so we needed an efficient way to compute it. Using the commutative property of Hadamard products: \( (A \otimes B) \circ ( C \otimes D) = (A\circ C) \otimes (B\circ D) \), we can find a "sum-factorization" type expression for \( A \circ U \), where here \( A = A_x \otimes A_y \otimes A_z \) and \( U \) is an \( n \times n \) matrix.
We make use of the flux basis being collocated on flux nodes, so the directions that aren't the derivative are identity. (Note that weights can be a diagonal matrix not necessarily identity). This results in the Hadamard product only needing \( \mathcal{O}(n^{d+1})\) flops to compute non-zero entries.
This is NOT for GENERAL Hadamard products since those are \( \mathcal{O}(n^{2d})\) .

Definition at line 798 of file operators.cpp.

Member Data Documentation

oneD_surf_operator

template<int dim, int n_faces>

std::array<dealii::FullMatrix<double>,2> PHiLiP::OPERATOR::SumFactorizedOperators< dim, n_faces >::oneD_surf_operator

Stores the one dimensional surface operator.

Note that in 1D there are only 2 faces

Definition at line 385 of file operators.h.

---

The documentation for this class was generated from the following files:

• src/operators/operators.h
• src/operators/operators.cpp