helpers.spectral_helper module

class ChebychevHelper(*args, x0=-1, x1=1, **kwargs)

Bases: SpectralHelper1D

The Chebychev base consists of special kinds of polynomials, with the main advantage that you can easily transform between physical and spectral space by discrete cosine transform. The differentiation in the Chebychev T base is dense, but can be preconditioned to yield a differentiation operator that moves to Chebychev U basis during differentiation, which is sparse. When using this technique, problems need to be formulated in first order formulation.

This implementation is largely based on the Dedalus paper (https://doi.org/10.1103/PhysRevResearch.2.023068).

get_1dgrid()

Generates a 1D grid with Chebychev points. These are clustered at the boundary. You need this kind of grid to use discrete cosine transformation (DCT) to get the Chebychev representation. If you want a different grid, you need to do an affine transformation before any Chebychev business.

Returns:

1D grid

Return type:

numpy.ndarray

get_BC(kind, **kwargs)

Get boundary condition row for boundary bordering. kwargs will be passed on to implementations of the BC of the kind you choose. Specifically, x for ‘dirichlet’ boundary condition, which is the coordinate at which to set the BC.

Parameters:

kind ('integral' or 'dirichlet') – Kind of boundary condition you want

get_Dirichlet_BC_row(x)

Get a row for generating Dirichlet BCs at x with T polynomials. It returns the values of the T polynomials at x.

Parameters:

x (float) – Position of the boundary condition

Returns:

Row to put into a matrix

Return type:

self.xp.ndarray

get_Dirichlet_recombination_matrix()

Get matrix for Dirichlet recombination, which changes the basis to have sparse boundary conditions. This makes for a good right preconditioner.

Returns:

Sparse conversion matrix

Return type:

scipy.sparse

get_Neumann_BC_row(x)

Get a row for generating Neumann BCs at x with T polynomials.

Parameters:

x (float) – Position of the boundary condition

Returns:

Row to put into a matrix

Return type:

self.xp.ndarray

get_basis_change_matrix(conv='T2T', **kwargs)

As the differentiation matrix in Chebychev-T base is dense but is sparse when simultaneously changing base to Chebychev-U, you may need a basis change matrix to transfer the other matrices as well. This function returns a conversion matrix from ChebychevHelper.get_conv. Not that **kwargs are used to absorb arguments for other bases, see documentation of SpectralHelper1D.get_basis_change_matrix.

Parameters:

conv (str) – Conversion code, i.e. T2U

Returns:

Sparse conversion matrix

get_conv(name, N=None)

Get conversion matrix between different kinds of polynomials. The supported kinds are

• T: Chebychev polynomials of first kind

• U: Chebychev polynomials of second kind

• D: Dirichlet recombination.

You get the desired matrix by choosing a name as A2B. I.e. T2U for the conversion matrix from T to U. Once generates matrices are cached. So feel free to call the method as often as you like.

Parameters:

• name (str) – Conversion code, e.g. ‘T2U’

• N (int) – Size of the matrix (optional)

Returns:

Sparse conversion matrix

Return type:

scipy.sparse

get_differentiation_matrix(p=1)

Keep in mind that the T2T differentiation matrix is dense.

Parameters:

p (int) – Derivative you want to compute

Returns:

Differentiation matrix

Return type:

numpy.ndarray

get_integ_BC_row()

Get a row for generating integral BCs with T polynomials. It returns the values of the integrals of T polynomials over the entire interval.

Returns:

Row to put into a matrix

Return type:

self.xp.ndarray

get_integration_matrix(lbnd=0)

Get matrix for integration

Parameters:

lbnd (float) – Lower bound for integration, only 0 is currently implemented

Returns:

Sparse integration matrix

get_norm(N=None)

Get normalization for converting Chebychev coefficients and DCT

Parameters:

N (int, optional) – Resolution

Returns:

Normalization

Return type:

self.xp.array

get_wavenumbers()

Get the domain in spectral space

itransform(u, *args, axes=None, shape=None, **kwargs)

Inverse DCT along axis.

Parameters:

• u – Data you want to transform

• axes (tuple) – Axes you want to transform along

Returns:

Data in physical space

transform(u, *args, axes=None, shape=None, **kwargs)

DCT along axes. kwargs will be passed on to the FFT library.

Parameters:

• u – Data you want to transform

• axes (tuple) – Axes you want to transform along

Returns:

Data in spectral space

class FFTHelper(*args, x0=0, x1=6.283185307179586, **kwargs)

Bases: SpectralHelper1D

distributable = True

get_1dgrid()

We use equally spaced points including the left boundary and not including the right one, which is the left boundary.

get_BC(kind)

Get a sort of boundary condition. You can use kind=integral, to fix the integral, or you can use kind=Nyquist. The latter is not really a boundary condition, but is used to set the Nyquist mode to some value, preferably zero. You should set the Nyquist mode zero when the solution in physical space is real and the resolution is even.

Parameters:

kind ('integral' or 'nyquist') – Kind of BC

Returns:

Boundary condition row

Return type:

self.xp.ndarray

get_Nyquist_mode_index()

Compute the index of the Nyquist mode, i.e. the mode with the lowest wavenumber, which doesn’t have a positive counterpart for even resolution. This means real waves of this wave number cannot be properly resolved and you are best advised to set this mode zero if representing real functions on even-resolution grids is what you’re after.

Returns:

Index of the Nyquist mode

Return type:

int

get_differentiation_matrix(p=1)

This matrix is diagonal, allowing to invert concurrently.

Parameters:

p (int) – Order of the derivative

Returns:

sparse differentiation matrix

get_integ_BC_row()

Only the 0-mode has non-zero integral with FFT basis in periodic BCs

get_integration_matrix(p=1)

Get integration matrix to compute p-th integral over the entire domain.

Parameters:

p (int) – Order of integral you want to compute

Returns:

sparse integration matrix

get_plan(u, forward, *args, **kwargs)

get_wavenumbers()

Be careful that this ordering is very unintuitive.

itransform(u, *args, axes=None, shape=None, **kwargs)

Inverse FFT.

Parameters:

• u – Data you want to transform

• axes (tuple) – Axes over which to transform

Returns:

transformed data

transform(u, *args, axes=None, shape=None, **kwargs)

FFT along axes. kwargs are passed on to the FFT library.

Parameters:

• u – Data you want to transform

• axes (tuple) – Axes you want to transform over

Returns:

transformed data

class SpectralHelper(comm=None, useGPU=False, debug=False)

Bases: object

This class has three functions:

• Easily assemble matrices containing multiple equations

• Direct product of 1D bases to solve problems in more dimensions

• Distribute the FFTs to facilitate concurrency.

comm

MPI communicator

Type:

mpi4py.Intracomm

debug

Perform additional checks at extra computational cost

Type:

bool

useGPU

Whether to use GPUs

Type:

bool

axes

List of 1D bases

Type:

list

components

List of strings of the names of components in the equations

Type:

list

full_BCs

List of Dictionaries containing all information about the boundary conditions

Type:

list

BC_mat

List of lists of sparse matrices to put BCs into and eventually assemble the BC matrix from

Type:

list

BCs

Matrix containing only the BCs

Type:

sparse matrix

fft_cache

Cache FFTs of various shapes here to facilitate padding and so on

Type:

dict

BC_rhs_mask

Mask values that contain boundary conditions in the right hand side

Type:

self.xp.ndarray

BC_zero_index

Indeces of rows in the matrix that are replaced by BCs

Type:

self.xp.ndarray

BC_line_zero_matrix

Matrix that zeros rows where we can then add the BCs in using BCs

Type:

sparse matrix

rhs_BCs_hat

Boundary conditions in spectral space

Type:

self.xp.ndarray

global_shape

Global shape of the solution as in mpi4py-fft

Type:

tuple

fft_obj

When using distributed FFTs, this will be a parallel transform object from mpi4py-fft

init

This is the same init that is used throughout the problem classes

Type:

tuple

init_forward

This is the equivalent of init in spectral space

Type:

tuple

property V

Get domain volume

add_BC(component, equation, axis, kind, v, line=-1, scalar=False, **kwargs)

Add a BC to the matrix. Note that you need to convert the list of lists of BCs that this method generates to a single sparse matrix by calling setup_BCs after adding/removing all BCs. Forward arguments for the boundary conditions using kwargs. Refer to documentation of 1D bases for details.

Parameters:

• component (str) – Name of the component the BC should act on

• equation (str) – Name of the equation for the component you want to put the BC in

• axis (int) – Axis you want to add the BC to

• kind (str) – kind of BC, e.g. Dirichlet

• v – Value of the BC

• line (int) – Line you want the BC to go in

• scalar (bool) – Put the BC in all space positions in the other direction

add_axis(base, *args, **kwargs)

Add an axis to the domain by deciding on suitable 1D base. Arguments to the bases are forwarded using *args and **kwargs. Please refer to the documentation of the 1D bases for possible arguments.

Parameters:

base (str) – 1D spectral method

add_component(name)

Add solution component(s).

Parameters:

name (str or list of strings) – Name(s) of component(s)

add_equation_lhs(A, equation, relations)

Add the left hand part (that you want to solve implicitly) of an equation to a list of lists of sparse matrices that you will convert to an operator later.

Example: Setup linear operator L for 1D heat equation using Chebychev method in first order form and T-to-U preconditioning:

Parameters:

• A (list of lists of sparse matrices) – The operator to be

• equation (str) – The equation of the component you want this in

• relations – (dict): Relations between quantities

check_BCs(u)

Check that the solution satisfies the boundary conditions

Parameters:

u – The solution you want to check

convert_operator_matrix_to_operator(M)

Promote the list of lists of sparse matrices to a single sparse matrix that can be used as linear operator. See documentation of SpectralHelper.add_equation_lhs for an example.

Parameters:

M (list of lists of sparse matrices) – The operator to be

Returns:

sparse linear operator

dtype

alias of mesh

eliminate_zeros(A)

Eliminate zeros from sparse matrix. This can reduce memory footprint of matrices somewhat. Note: At the time of writing, there are memory problems in the cupy implementation of eliminate_zeros. Therefore, this function copies the matrix to host, eliminates the zeros there and then copies back to GPU.

Parameters:

A – sparse matrix to be pruned

Returns:

CSC sparse matrix

expand_matrix_ND(matrix, aligned)

fft_backend = 'scipy'

fft_comm_backend = 'MPI'

fft_lib = <module 'scipy.fft' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/fft/__init__.py'>

get_BC(axis, kind, line=-1, scalar=False, **kwargs)

Use this method for boundary bordering. It gets the respective matrix row and embeds it into a matrix. Pay attention that if you have multiple BCs in a single equation, you need to put them in different lines. Typically, the last line that does not contain a BC is the best choice. Forward arguments for the boundary conditions using kwargs. Refer to documentation of 1D bases for details.

Parameters:

• axis (int) – Axis you want to add the BC to

• kind (str) – kind of BC, e.g. Dirichlet

• line (int) – Line you want the BC to go in

• scalar (bool) – Put the BC in all space positions in the other direction

Returns:

sparse matrix containing the BC

get_Dirichlet_recombination_matrix(axis=-1)

Get Dirichlet recombination matrix along axis. Not that it only makes sense in directions discretized with variations of Chebychev bases.

Parameters:

axis (int) – Axis you discretized with Chebychev

Returns:

sparse matrix

get_Id()

Get identity matrix

Returns:

sparse identity matrix

get_basis_change_matrix(axes=None, **kwargs)

Some spectral bases do a change between bases while differentiating. This method returns matrices that changes the basis to whatever you want. Refer to the methods of the same name of the 1D bases to learn what parameters you need to pass here as kwargs.

Parameters:

axes (tuple) – Axes along which to change basis.

Returns:

sparse basis change matrix

get_differentiation_matrix(axes, **kwargs)

Get differentiation matrix along specified axis. kwargs are forwarded to the 1D base implementation.

Parameters:

axes (tuple) – Axes along which to differentiate.

Returns:

sparse differentiation matrix

get_empty_operator_matrix(diag=False)

Return a matrix of operators to be filled with the connections between the solution components.

Parameters:

diag (bool) – Whether operator is block-diagonal

Returns:

list containing sparse zeros

get_fft(axes=None, direction='object', padding=None, shape=None)

When using MPI, we use PFFT objects generated by mpi4py-fft

Parameters:

• axes (tuple) – Axes you want to transform over

• direction (str) – use “forward” or “backward” to get functions for performing the transforms or “object” to get the PFFT object

• padding (tuple) – Padding for dealiasing

• shape (tuple) – Shape of the transform

Returns:

transform

get_filter_matrix(axis, **kwargs)

Get bandpass filter along axis. See the documentation get_filter_matrix in the 1D bases for what kwargs are admissible.

Returns:

sparse bandpass matrix

get_grid(forward_output=False)

Get grid in physical space

get_indices(forward_output=True)

get_integration_matrix(axes)

Get integration matrix to integrate along specified axis.

Parameters:

axes (tuple) – Axes along which to integrate over.

Returns:

sparse integration matrix

get_local_slice_of_1D_matrix(M, axis)

Get the local version of a 1D matrix. When using distributed FFTs, each rank will carry only a subset of modes, which you can sort out via the SpectralHelper.local_slice() attribute. When constructing a 1D matrix, you can use this method to get the part corresponding to the modes carried by this rank.

Parameters:

• M (sparse matrix) – Global 1D matrix you want to get the local version of

• axis (int) – Direction in which you want the local version. You will get the global matrix in other directions.

Returns:

sparse local matrix

get_pfft(axes=None, padding=None, grid=None)

get_wavenumbers()

Get grid in spectral space

global_slice(forward_output=True)

index(name)

Get the index of component name.

Parameters:

name (str or list of strings) – Name(s) of component(s)

Returns:

Index of the component

Return type:

int

infer_alignment(u, forward_output, padding=None, **kwargs)

itransform(u, *args, axes=None, padding=None, **kwargs)

linalg = <module 'scipy.sparse.linalg' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/sparse/linalg/__init__.py'>

local_slice(forward_output=True)

property ncomponents

property ndim

newDistArray(pfft=None, forward_output=True, val=0, rank=1, view=False)

Get an empty distributed array. This is almost a copy of the function of the same name from mpi4py-fft, but takes care of all the solution components in the tensor.

put_BCs_in_matrix(A)

Put the boundary conditions in a matrix by replacing rows with BCs.

put_BCs_in_rhs(rhs)

Put the BCs in the right hand side for solving. This function will transform along each axis individually and add all BCs in that axis. Consider put_BCs_in_rhs_hat to add BCs with no extra transforms needed.

Parameters:

rhs – Right hand side in physical space

Returns:

rhs in physical space with BCs

put_BCs_in_rhs_hat(rhs_hat)

Put the BCs in the right hand side in spectral space for solving. This function needs no transforms and caches a mask for faster subsequent use.

Parameters:

rhs_hat – Right hand side in spectral space

Returns:

rhs in spectral space with BCs

redistribute(u, axis, forward_output, **kwargs)

remove_BC(component, equation, axis, kind, line=-1, scalar=False, **kwargs)

Remove a BC from the matrix. This is useful e.g. when you add a non-scalar BC and then need to selectively remove single BCs again, as in incompressible Navier-Stokes, for instance. Forwards arguments for the boundary conditions using kwargs. Refer to documentation of 1D bases for details.

Parameters:

• component (str) – Name of the component the BC should act on

• equation (str) – Name of the equation for the component you want to put the BC in

• axis (int) – Axis you want to add the BC to

• kind (str) – kind of BC, e.g. Dirichlet

• v – Value of the BC

• line (int) – Line you want the BC to go in

• scalar (bool) – Put the BC in all space positions in the other direction

setup_BCs()

Convert the list of lists of BCs to the boundary condition operator. Also, boundary bordering requires to zero out all other entries in the matrix in rows containing a boundary condition. This method sets up a suitable sparse matrix to do this.

classmethod setup_CPU(useFFTW=False)

switch to CPU modules

classmethod setup_GPU()

switch to GPU modules

setup_fft(real_spectral_coefficients=False)

This function must be called after all axes have been setup in order to prepare the local shapes of the data. This must also be called before setting up any BCs.

Parameters:

real_spectral_coefficients (bool) – Allow only real coefficients in spectral space

property shape

Get shape of individual solution component

sparse_lib = <module 'scipy.sparse' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/sparse/__init__.py'>

transform(u, *args, axes=None, padding=None, **kwargs)

property u_init

Get empty data container in physical space

property u_init_forward

Get empty data container in spectral space

property u_init_physical

Get empty data container in physical space

xp = <module 'numpy' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/numpy/__init__.py'>

class SpectralHelper1D(N, x0=None, x1=None, useGPU=False, useFFTW=False)

Bases: object

Abstract base class for 1D spectral discretizations. Defines a common interface with parameters and functions that all bases need to have.

When implementing new bases, please take care to use the modules that are supplied as class attributes to enable the code for GPUs.

N

Resolution

Type:

int

x0

Coordinate of left boundary

Type:

float

x1

Coordinate of right boundary

Type:

float

L

Length of the domain

Type:

float

useGPU

Whether to use GPUs

Type:

bool

distributable = False

fft_lib = <module 'scipy.fft' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/fft/__init__.py'>

get_1dgrid()

Get the grid in physical space

Returns:

Grid

Return type:

self.xp.array

get_BC(kind)

To facilitate boundary conditions (BCs) we use either a basis where all functions satisfy the BCs automatically, as is the case in FFT basis for periodic BCs, or boundary bordering. In boundary bordering, specific lines in the matrix are replaced by the boundary conditions as obtained by this method.

Parameters:

• kind (str) – The type of BC you want to implement please refer to the implementations of this method in the

• implemented (individual 1D bases for what is)

Returns:

Boundary condition

Return type:

self.xp.array

get_Id()

Get identity matrix

Returns:

sparse diagonal identity matrix

get_basis_change_matrix(*args, **kwargs)

Some spectral discretization change the basis during differentiation. This method can be used to transfer between the various bases.

This method accepts arbitrary arguments that may not be used in order to provide an easy interface for multi- dimensional bases. For instance, you may combine an FFT discretization with an ultraspherical discretization. The FFT discretization will always be in the same base, but the ultraspherical discretization uses a different base for every derivative. You can then ask all bases for transfer matrices from one ultraspherical derivative base to the next. The FFT discretization will ignore this and return an identity while the ultraspherical discretization will return the desired matrix. After a Kronecker product, you get the 2D version of the matrix you want. This is what the SpectralHelper does when you call the method of the same name on it.

Returns:

sparse bases change matrix

get_differentiation_matrix()

get_empty_operator_matrix(S, O)

Return a matrix of operators to be filled with the connections between the solution components.

Parameters:

• S (int) – Number of components in the solution

• O (sparse matrix) – Zero matrix used for initialization

Returns:

list of lists containing sparse zeros

get_filter_matrix(kmin=0, kmax=None)

Get a bandpass filter.

Parameters:

• kmin (int) – Lower limit of the bandpass filter

• kmax (int) – Upper limit of the bandpass filter

Returns:

sparse matrix

get_integration_matrix()

get_wavenumbers()

Get the grid in spectral space

get_zero()

Get a matrix with all zeros of the correct size.

Returns:

sparse matrix with zeros everywhere

linalg = <module 'scipy.sparse.linalg' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/sparse/linalg/__init__.py'>

classmethod setup_CPU(useFFTW=False)

switch to CPU modules

classmethod setup_GPU()

switch to GPU modules

sparse_lib = <module 'scipy.sparse' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/scipy/sparse/__init__.py'>

xp = <module 'numpy' from '/home/runner/micromamba/envs/pySDC/lib/python3.13/site-packages/numpy/__init__.py'>

class UltrasphericalHelper(*args, x0=-1, x1=1, **kwargs)

Bases: ChebychevHelper

This implementation follows https://doi.org/10.1137/120865458. The ultraspherical method works in Chebychev polynomials as well, but also uses various Gegenbauer polynomials. The idea is that for every derivative of Chebychev T polynomials, there is a basis of Gegenbauer polynomials where the differentiation matrix is a single off-diagonal. There are also conversion operators from one derivative basis to the next that are sparse.

This basis is used like this: For every equation that you have, look for the highest derivative and bump all matrices to the correct basis. If your highest derivative is 2 and you have an identity, it needs to get bumped from 0 to 1 and from 1 to 2. If you have a first derivative as well, it needs to be bumped from 1 to 2. You don’t need the same resulting basis in all equations. You just need to take care that you translate the right hand side to the correct basis as well.

get_S(lmbda)

Get matrix for bumping the derivative base by one from lmbda to lmbda + 1. This is the same language as in https://doi.org/10.1137/120865458.

Parameters:

lmbda (int) – Ingoing derivative base

Returns:

Conversion from derivative base lmbda to lmbda + 1

Return type:

sparse matrix

get_basis_change_matrix(p_in=0, p_out=0, **kwargs)

Get a conversion matrix from derivative base p_in to p_out.

Parameters:

• p_out (int) – Resulting derivative base

• p_in (int) – Ingoing derivative base

get_differentiation_matrix(p=1)

Notice that while sparse, this matrix is not diagonal, which means the inversion cannot be parallelized easily.

Parameters:

p (int) – Order of the derivative

Returns:

sparse differentiation matrix

get_integration_constant(u_hat, axis)

Get integration constant for lower bound of -1. See documentation of UltrasphericalHelper.get_integration_matrix for details.

Parameters:

• u_hat – Solution in spectral space

• axis – Axis you want to integrate over

Returns:

Integration constant, has one less dimension than u_hat

get_integration_matrix()

Get an integration matrix. Please use UltrasphericalHelper.get_integration_constant afterwards to compute the integration constant such that integration starts from x=-1.

Example:

import numpy as np
from pySDC.helpers.spectral_helper import UltrasphericalHelper

N = 4
helper = UltrasphericalHelper(N)
coeffs = np.random.random(N)
coeffs[-1] = 0

poly = np.polynomial.Chebyshev(coeffs)

S = helper.get_integration_matrix()
U_hat = S @ coeffs
U_hat[0] = helper.get_integration_constant(U_hat, axis=-1)

assert np.allclose(poly.integ(lbnd=-1).coef[:-1], U_hat)

Returns:

sparse integration matrix

cache(func)

Decorator for caching return values of functions. This is very similar to functools.cache, but without the memory leaks (see https://docs.astral.sh/ruff/rules/cached-instance-method/).

Example:

num_calls = 0

@cache
def increment(x):
    num_calls += 1
    return x + 1

increment(0)  # returns 1, num_calls = 1
increment(1)  # returns 2, num_calls = 2
increment(0)  # returns 1, num_calls = 2

Parameters:

func (function) – The function you want to cache the return value of

Returns:

return value of func