Spectral Toolbox¶
Objects for implementing spectral methods.
This module provides several methods contained within the SpectralToolbox class for performing standard vector calculus operations in spectral space. It is intended for use on a two-dimensional domain of equally-spaced gridpoints.
Classes¶
- SpectralToolbox – Superclass to implement standard techniques of spectral methods
Miscellaneous Functions for Unit Testing¶
- test_func – Test function for unit testing
- check_derivatives
- check_jacobian
- check_multiply
- check_laplacian
All of the above compare the computed solution to an analytical solution.
Notes¶
- Invoking this module from the command line without any arguments will perform unit tests on the methods contained within.
- The Fourier representation used is:
The spectrum is stored in the logical-to-mathematicians order, i.e. the wavenumbers are stored as:
| (-N/2)(-N/2) | (-N/2)(-N/2 + 1) | ... | (-N/2)(N/2 - 1) |
| (-N/2 + 1)(-N/2) | (-N/2 + 1)(-N/2 + 1) | ... | (-N/2 + 1)(N/2 -1) |
| ... | ... | ... | ... |
| (N/2 - 1)(-N/2) | (N/2 - 1)(-N/2 + 1) | ... | (N/2 - 1)(N/2 -1) |
With the Nyquist frequencies (corresponding to the most negative frequencies) located where one or both index = -N/2, and the constant mode corresponding to \(k_{1} = 0\), \(k_{2} = 0\) in the (N, N) position in the array.
2) The sign matrix, sign_mat, is used to relate the truncated Fourier series of a function, with wavenumbers that run from [-N/2, N/2), with the FFT, with wavenumbers in [0, N). The point of doing this is to make the operations in spectral space more intuitive from a mathematical perspective. We proceed by discretising the desired function on an equispaced grid as above, noting that the sum runs from k = -N/2 to N/2 -1. In order to make the sum correspond with the FFT, we perform the change of variables \(p = k + N/2\):
Note that a factor of \((-1)^{p}\) is obtained. These alternating signs are implemented in the sign matrix to permit a single multiplication operation.
See also¶
numpy.fft, numpy
-
class
spectral_toolbox.SpectralToolbox(N, L)[source]¶ Implements the required vector calculus techniques for spectral methods.
This class stores the required constants as attributes to be used by its methods for spectral methods in 2 spatial dimensions.
Attributes
- N : The number of grid points in the domain, such that there are N X N equally-spaced points
- L : The side length of the domain
- factor : Factor arising in spectral differentiation on non-2pi domain (see below)
- sign_mat : Matrix to reorder spectral space (see below)
- deriv_mat_x1, deriv_mat_x2` : matrices for spectral differentiation in x1 and x2 directions (see calc_derivative method)
- laplace_op : Array to implement Laplacian (see laplacian method)
Methods
- dealias_pad – Pads the input array with zeros to prevent aliasing
- dealias_unpad – Removes the padding from dealias_pad
- forward_fft – Wrapper to perform Fast Fourier Transform with desired normalisation
- inverse_fft – Wrapper to invert FFT with desired normalisation
- calc_derivative – Computes 1st or 2nd order derivatives along x and/or y directions
- multiply_nonlinear – Multiplies two functions in a pseudo-spectral fashion
- solve_inverse_laplacian – Solves the inverse Laplacian problem with a chosen constant of integration
- jacobian1 – Computes Jacobian of a single function
- jacobian – Computes Jacobian of two functions
- laplacian – Computes the Laplacian of a given function
Example
>> st = SpectralToolbox(128, 2*np.pi)>> du_dx = st.calc_derivative(A, 'x') # Compute x-derivative-
calc_derivative(array_in, direction1, direction2=False)[source]¶ Computes 1st or 2nd order derivatives in x and/or y directions.
This method implements spectral (spatial) differentiation through the use of the derivative matrices, stored as attributes of the SpectralToolbox class. The Nyquist frequency is treated specially for odd orders of multiplication.
Parameters
- array_in : the array of spectral coefficients to be differentiated
- direction1 : the first direction to be differentiated (‘x’ or ‘y’)
- direction2 : the optional second direction to be differentiated (‘x’ or ‘y’)
Returns
- out : the differentiated spectrum
Notes
Recall that the function, u, is represented in Fourier space as:
\[u(x, y, t) = \sum_{k_{1} = -N/2}^{N/2 -1} \sum_{k_{2} = -N/2}^{N/2 - 1} \hat{u}(t) \exp((2i\pi/L)(k_{1}x + k_{2}y))\]Then the first derivative in the x-direction, as \(\hat{u}(t)\) is not a function of space, is simply:
\[\frac{\partial u(x, y, t)}{\partial x} = \sum_{k_{1} = -N/2}^{N/2 -1} \sum_{k_{2} = -N/2}^{N/2 - 1} \hat{u}(t) (2i\pi k_{1}/L) \exp((2i\pi/L)(k_{1}x + k_{2}y))\]Thus, differentiation in spectral space is a matter of multiplying each Fourier mode by the length factor (\(B=2\pi/L\)) times i times the associated wavenumber. To reduce loops and make the implementation simpler, these modes are contained in the derivatrive matrices, which take the form:
\[\begin{split}\text{deriv\_mat\_x1} = \begin{array}{cccc} iBk_{1} & iBk_{1} & \cdots & iBk_{1} \\ iBk_{2} & iBk_{2} & \cdots & iBk_{2} \\ \vdots & \vdots & \ddots & \vdots \\ iBk_{n} & iBk_{n} & \cdots & iBk_{n} \end{array}\end{split}\]and:
\[\begin{split}\text{deriv\_mat\_x2} = \begin{array}{cccc} iBk_{1} & iBk_{2} & \cdots & iBk_{n} \\ iBk_{1} & iBk_{2} & \cdots & iBk_{n}\\ \vdots & \vdots & \ddots & \vdots \\ iBk_{1} & iBk_{2} & \cdots & iBk_{n} \end{array}\end{split}\]Example
>> A = array([[ 1.+0.j, 1.+0.j],[ 1.+0.j, 1.+0.j]])>> st.calc_derivative(A, 'x')>> array([[ 0.+0.j, 0.+0.j],[ 0.+0.j, 0.+0.j]])
-
dealias_pad(A)[source]¶ Pads the spectrum with zeros in high frequencies for aliasing control.
This method doubles the size of the spectrum and sets all high frequencies to zero for aliasing control. This causes the resolution in real space (after inverse FFT) to be doubled. The higher frequencies are then discarded upon returning to spectral space. The array is doubled in a symmetric fashion.
Parameters
- A : the spectrum to be padded
Returns
- B : the padded spectrum
See also
dealias_unpad
Example
>> A = np.ones((2,2), dtype = complex)>> st.dealias_pad(A)>> array([[ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],[ 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],[ 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],[ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])
-
dealias_unpad(A)[source]¶ Unpads the spectrum, removing high frequencies, for aliasing control.
This method reverses the padding set up in the dealias_pad method, returning the size of the spectrum to the initial value upon returning to spectral space. The highest half of the frequencies are discarded.
Parameters
- A : the spectrum to be unpadded
Returns
- B : the unpadded spectrum
See also
dealias_pad
Example
>> A = array([[ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],[ 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],[ 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],[ 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])>> st.dealias_unpad(A)>> array([[ 1.+0.j, 1.+0.j],[ 1.+0.j, 1.+0.j]])
-
forward_fft(array_in)[source]¶ Wrapper to implement the forward Fast Fourier Transform.
This method implements the forward FFT (i.e. from real-space to spectral space). The values are normalised such that the spectral coefficient corresponds to a wave of exactly that amplitude, rather than the normalised value returned by the standard FFT. The sign matrix is used to yield a spectrum running from -N/2 to N/2 -1.
Parameters
- array_in : the array of real values
Returns
- out : the spectral coefficients, in the chosen framework
See Also
Note 2 in the class header (sign matrix), inverse_fft
-
inverse_fft(array_in)[source]¶ Wrapper to implement the inverse Fast Fourier Transform.
This method implements the inverse FFT (i.e. from spectral space to real space). The values are normalised such that the spectral coefficient corresponds to a wave of exactly that amplitude, rather than the normalised value returned by the standard FFT. The sign matrix is used to yield a spectrum running from -N/2 to N/2 -1.
Parameters
- array_in : the array of spectral coefficients
Returns
- out : the real values at gridpoints, in the chosen framework
See Also
Note 2 in the class header (sign matrix), forward_fft
-
jacobian(A, B)[source]¶ Computes the Jacobian for a two scalar-valued functions.
This method implements the Jacobian for two scalar-valued functions. The Jacobian may be written as:
\[J(A, B) = A_{x}B_{y} - A_{y}B_{x}\]where subscripts denote partial derivatives.
Parameters
- A, B : The Fourier representations of the functions A and B (above). Must have the same dimensions.
Returns
- J : The Fourier modes of the Jacobian
-
jacobian1(A)[source]¶ Computes the Jacobian for a single scalar-valued function.
This method implements the Jacobian for a single function, such as a streamfunction or velocity potential. The Jacobian in this case takes the form:
\[J(u) = u_{xx}u_{yy} - 2u_{xy}\]where subscripts denote partial derivatives.
Parameters
- A : The function of which the Jacobian is desired
Returns
- J : The Fourier modes of the Jacobian
-
laplacian(array_in)[source]¶ Computes the Laplacian of a given function in spectral space.
The Laplacian may be written as:
\[f = \nabla^{2}u = \nabla \cdot \nabla u\]Parameters
- array_in : the value of interest, in Fourier space
Returns
- array_out : the Laplacian of the input array
See Also
calc_derivative
-
multiply_nonlinear(array1, array2)[source]¶ Simple method for multiplying two sets of data when nonlinearities are involved and transformation into real space is necessary.
This method computes the product of two quantities defined in spectral space using a pseudo-spectral method, i.e. the multiplication is performed in realspace. The spectra are padded for aliasing control and must be of the same dimension.
Parameters
- array_1, array_2 : the input spectra to be multiplied
Returns
- f3_hat : the spectrum obtainined from the multiplication
See Also
dealias_pad, dealias_unpad
-
solve_inverse_laplacian(array_in, zero_mode_val)[source]¶ Solves the inverse Laplacian problem in spectral space.
This method implements the inverse Laplacian, i.e. it finds u in:
\[\nabla^{2}u = f\]with a constant of integration which must be specified by the user.
Parameters
- array_in : the known spectrum, f
- zero_mode_val : the constant of integration, corresponding to the constant Fourier mode
Returns
- array_out : the desired quantity, u
See Also
calc_derivative
