# shenfun package¶

## shenfun.la module¶

This module contains linear algebra solvers for SparseMatrixes

class shenfun.la.NeumannSolve(A, test)[source]

Bases: object

Solver class for matrix created by Neumann bases

Assuming Neumann test- and trialfunction, where index k=0 is used only to fix the mean value.

Parameters
• A (SparseMatrix)

• test (BasisFunction)

__call__(b, u=None, axis=0)[source]

Solve matrix problem A u = b

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

If u is not provided, then b is overwritten with the solution and returned

class shenfun.la.PDMA(mat, solver='cython')[source]

Bases: object

Parameters
• mat (SparseMatrix) – Symmetric pentadiagonal matrix with diagonals in offsets -4, -2, 0, 2, 4

• solver (str, optional) – Choose implementation

• cython - Use efficient cython implementation

• python - Use python/scipy

static PDMA_SymLU(*args, **kwargs)[source]

Symmetric LU decomposition

static PDMA_SymSolve(*args, **kwargs)[source]

Symmetric solve (for testing only)

__call__(b, u=None, axis=0)[source]

Solve matrix problem self u = b

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

Note

If u is not provided, then b is overwritten with the solution and returned

init()[source]

Initialize and allocate solver

class shenfun.la.Solve(A, test)[source]

Bases: object

Solver class for matrix created by Dirichlet bases

Possibly with inhomogeneous boundary values

Parameters
• A (SparseMatrix)

• test (BasisFunction)

__call__(b, u=None, axis=0)[source]

Solve matrix problem Au = b

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

Note

If u is not provided, then b is overwritten with the solution and returned

class shenfun.la.Solver2D(tpmats)[source]

Bases: object

Generic solver for tensorproductspaces in 2D

Parameters

mats (sequence) – sequence of instances of TPMatrix

__call__(b, u=None)[source]

Call self as a function.

matvec(u, c)[source]
class shenfun.la.SolverGeneric1ND(mats)[source]

Bases: object

Generic solver for tensorproduct matrices consisting of non-diagonal matrices along only one axis.

Parameters

mats (sequence) – sequence of instances of TPMatrix

Note

In addition to the one non-diagonal direction, the solver can also handle up to two diagonal (Fourier) directions. Also, this Python version of the solver is not very efficient. Consider implementing in Cython.

__call__(b, u=None)[source]

Call self as a function.

class shenfun.la.SolverGeneric2ND(tpmats)[source]

Bases: object

Generic solver for problems consisting of tensorproduct matrices containing two non-diagonal submatrices.

Parameters

mats (sequence) – sequence of instances of TPMatrix

Note

In addition to two non-diagonal matrices, the solver can also handle one additional diagonal matrix (one Fourier matrix).

__call__(b, u=None, format='csr')[source]

Call self as a function.

diags(i, format='csr')[source]

Return matrix for given index i in diagonal direction

get_diagonal_axis()[source]
matvec(u, c)[source]
class shenfun.la.TDMA(mat)[source]

Bases: object

Tridiagonal matrix solver

Parameters

mat (SparseMatrix) – Symmetric tridiagonal matrix with diagonals in offsets -2, 0, 2

static TDMA_SymLU(*args, **kwargs)[source]
static TDMA_SymSolve(*args, **kwargs)[source]
__call__(b, u=None, axis=0)[source]

Solve matrix problem self u = b

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

Note

If u is not provided, then b is overwritten with the solution and returned

init()[source]

Initialize and allocate solver

class shenfun.la.TDMA_O(mat)[source]

Bases: object

Tridiagonal matrix solver

Parameters

mat (SparseMatrix) – Symmetric tridiagonal matrix with diagonals in offsets -1, 0, 1

static TDMA_O_SymLU(*args, **kwargs)[source]
static TDMA_O_SymSolve(*args, **kwargs)[source]
__call__(b, u=None, axis=0)[source]

Solve matrix problem self u = b

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

Note

If u is not provided, then b is overwritten with the solution and returned

init()[source]

Initialize and allocate solver

## shenfun.matrixbase module¶

This module contains classes for working with sparse matrices.

The sparse matrices are computed as inner products of forms containing test and trial functions. These basis functions are chosen from the following, where we denote the $$k$$’th basis function of basis V as $$\phi_k$$:

Chebyshev basis:

Chebyshev basis of first kind

$\begin{split}\phi_k &= T_k \\ V &= span\{\phi_k\}_{k=0}^{N}\end{split}$

For homogeneous Dirichlet boundary conditions:

$\begin{split}\phi_k &= T_k - T_{k+2} \\ V &= span\{\phi_k\}_{k=0}^{N-2}\end{split}$

For homogeneous Neumann boundary conditions:

$\begin{split}\phi_k &= T_k - \left(\frac{k}{k+2}\right)^2T_{k+2} \\ V &= span\{\phi_k\}_{k=1}^{N-2}\end{split}$

For Biharmonic basis with both homogeneous Dirichlet and Neumann:

$\begin{split}\phi_k &= T_k - 2 \frac{k+2}{k+3} T_{k+2} + \frac{k+1}{k+3} T_{k+4} \\ V &= span\{\phi_k\}_{k=0}^{N-4}\end{split}$

The scalar product is computed as a weighted inner product with

$$w=1/\sqrt{1-x^2}$$ the weights.

Legendre basis:

Regular Legendre

$\begin{split}\phi_k &= L_k \\ V &= span\{\phi_k\}_{k=0}^{N}\end{split}$

Dirichlet boundary conditions

$\begin{split}\phi_k &= L_k-L_{k+2} \\ V &= span\{\phi_k\}_{k=0}^{N-2}\end{split}$

Homogeneous Neumann boundary conditions:

$\begin{split}\phi_k &= L_k - \frac{k(k+1)}{(k+2)(k+3)}L_{k+2} \\ V &= span\{\phi_k\}_{k=1}^{N-2}\end{split}$

Both homogeneous Dirichlet and Neumann:

$\begin{split}\psi_k &= L_k - 2 \frac{2k+5}{2k+7} L_{k+2} + \frac{2k+3}{2k+7} L_{k+4} \\ V &= span\{\phi_k\}_{k=0}^{N-4}\end{split}$
Laguerre basis:

Regular Laguerre function

$\begin{split}\phi_k &= L_k \cdot \exp(-x) \\ V &= span\{\phi_k\}_{k=0}^{N}\end{split}$

Homogeneous Dirichlet boundary conditions

$\begin{split}\phi_k &= (L_k-L_{k+2}) \cdot \exp(-x) \\ V &= span\{\phi_k\}_{k=0}^{N-1}\end{split}$

where $$L_k$$ is the Laguerre polynomial of order k.

Hermite basis:

Regular Hermite function

$\begin{split}\phi_k &= H_k \cdot \exp(-x^2/2)/(\pi^{0.25}\sqrt{2^k k!}) \\ V &= span\{\phi_k\}_{k=0}^{N}\end{split}$

where $$K_k$$ is the Hermite polynomial of order k. Homogeneous Dirichlet boundary conditions on (-inf, inf)

Fourier basis:

$\begin{split}\phi_k &= exp(ikx) \\ V &= span\{\phi_k\}_{k=-N/2}^{N/2-1}\end{split}$

Jacobi basis:

Regular Jacobi polynomials

$\begin{split}\phi_k(x) &= J_k(x, \alpha, \beta) \\ V &= span\{\phi_k\}_{k=0}^{N-1}\end{split}$

where $$\alpha > -1$$ and $$\beta > -1$$. $$J_k$$ is the regular Jacobi polynomial

Homogeneous Dirichlet boundary conditions

$\phi_k &= j_k(x, -1, -1) V &= span\{\phi_k\}_{k=0}^{N-3}$

where $$j_k$$ is the generalized Jacobi polynomial

Homogeneous Dirichlet and Neumann boundary conditions

$\phi_k &= j_k(x, -2, -2) V &= span\{\phi_k\}_{k=0}^{N-5}$

Homogeneous Dirichlet and first and second order derivatives (for 6th order equation)

$\phi_k &= j_k(x, -3, -3) V &= span\{\phi_k\}_{k=0}^{N-7}$
class shenfun.matrixbase.BlockMatrix(tpmats)[source]

Bases: object

A class for block matrices

Parameters

tpmats (sequence of TPMatrix or SparseMatrix) – The individual blocks for the matrix

Example

Stokes equations, periodic in x and y-directions

$\begin{split}-\nabla^2 u - \nabla p &= 0 \\ \nabla \cdot u &= 0 \\ u(x, y, z=\pm 1) &= 0\end{split}$

We use for the z-direction a Dirichlet basis (SD) and a regular basis with no boundary conditions (ST). This is combined with Fourier in the x- and y-directions (K0, K1), such that we get two TensorProductSpaces (TD, TT) that are the Cartesian product of these bases

$\begin{split}TD &= K0 \times K1 \times SD \\ TT &= K0 \times K1 \times ST\end{split}$

We choose trialfunctions $$u \in [TD]^3$$ and $$p \in TT$$, and then solve the weak problem

$\begin{split}\left( \nabla v, \nabla u\right) + \left(\nabla \cdot v, p \right) = 0\\ \left( q, \nabla \cdot u\right) = 0\end{split}$

for all $$v \in [TD]^3$$ and $$q \in TT$$.

To solve the problem we need to assemble a block matrix

$\begin{split}\begin{bmatrix} \left( \nabla v, \nabla u\right) & \left(\nabla \cdot v, p \right) \\ \left( q, \nabla \cdot u\right) & 0 \end{bmatrix}\end{split}$

This matrix is assemble below

>>> from shenfun import *
>>> from mpi4py import MPI
>>> comm = MPI.COMM_WORLD
>>> N = (24, 24, 24)
>>> K0 = FunctionSpace(N[0], 'Fourier', dtype='d')
>>> K1 = FunctionSpace(N[1], 'Fourier', dtype='D')
>>> SD = FunctionSpace(N[2], 'Legendre', bc=(0, 0))
>>> ST = FunctionSpace(N[2], 'Legendre')
>>> TD = TensorProductSpace(comm, (K0, K1, SD), axes=(2, 1, 0))
>>> TT = TensorProductSpace(comm, (K0, K1, ST), axes=(2, 1, 0))
>>> VT = VectorTensorProductSpace(TD)
>>> Q = MixedTensorProductSpace([VT, TD])
>>> up = TrialFunction(Q)
>>> vq = TestFunction(Q)
>>> u, p = up
>>> v, q = vq
>>> A01 = inner(div(v), p)
>>> A10 = inner(q, div(u))
>>> M = BlockMatrix(A00+A01+A10)

static apply_constraint(A, b, offset, i, constraint)[source]
diags(it=0, format='csr')[source]

Return global block matrix in scipy sparse format

For multidimensional forms the returned matrix is constructed for given indices in the periodic directions.

Parameters
• it (n-tuple of ints) – where n is dimensions. These are the indices into the scale arrays of the TPMatrices in various blocks. Should be zero along the non- periodic direction.

• format (str) – The format of the returned matrix. See Scipy sparse matrices

get_mats(return_first=False)[source]

Return flattened list of matrices in self

get_offset(i, axis=0)[source]
matvec(v, c)[source]

Compute matrix vector product

c = self * v

Parameters
Returns

c

Return type

Function

solve(b, u=None, constraints=(), return_system=False, Alu=None)[source]

Solve matrix system Au = b

where A is the current BlockMatrix (self)

Parameters
• b (array) – Array of right hand side

• u (array, optional) – Output array

• constraints (sequence of 3-tuples of (int, int, number)) – Any 3-tuple describe a dof to be constrained. The first int represents the block number of the function to be constrained. The second int gives which degree of freedom to constrain and the number gives the value it should obtain. For example, for the global restriction that

$\frac{1}{V}\int p dx = number$

where we have

$p = \sum_{k=0}^{N-1} \hat{p}_k \phi_k$

it is sufficient to fix the first dof of p, hat{p}_0, since the bases are created such that all basis functions except the first integrates to zero. So in this case the 3-tuple can be (2, 0, 0) if p is found in block 2 of the mixed basis.

The constraint can only be applied to bases with no given explicit boundary condition, like the pure Chebyshev or Legendre bases.

Other Parameters
• return_system (bool, optional) – If True then return the assembled block matrix as well as the solution in a 2-tuple (solution, matrix). This is helpful for repeated solves, because the returned matrix may then be factorized once and reused. Only for non-periodic problems

• Alu (pre-factorized matrix, optional) – Computed with Alu = splu(self), where self is the assembled block matrix. Only for non-periodic problems.

class shenfun.matrixbase.Identity(shape, scale=1)[source]

The identity matrix in SparseMatrix form

Parameters
• shape (2-tuple of ints) – The shape of the matrix

• scale (number, optional) – Scalar multiple of the matrix, defaults to unity

solve(b, u=None, axis=0)[source]

Solve matrix system Au = b

where A is the current matrix (self)

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multi- dimensional

Note

Vectors may be one- or multidimensional.

class shenfun.matrixbase.SparseMatrix(d, shape, scale=1.0)[source]

Bases: dict

Base class for sparse matrices.

The data is stored as a dictionary, where keys and values are, respectively, the offsets and values of the diagonals. In addition, each matrix is stored with a coefficient that is used as a scalar multiple of the matrix.

Parameters
• d (dict) – Dictionary, where keys are the diagonal offsets and values the diagonals

• shape (two-tuple of ints)

• scale (number, optional) – Scale matrix with this number

Examples

A tridiagonal matrix of shape N x N could be created as

>>> from shenfun import SparseMatrix
>>> import numpy as np
>>> N = 4
>>> d = {-1: 1, 0: -2, 1: 1}
>>> SparseMatrix(d, (N, N))
{-1: 1, 0: -2, 1: 1}


In case of variable values, store the entire diagonal. For an N x N matrix use

>>> d = {-1: np.ones(N-1),
...       0: -2*np.ones(N),
...       1: np.ones(N-1)}
>>> SparseMatrix(d, (N, N))
{-1: array([1., 1., 1.]), 0: array([-2., -2., -2., -2.]), 1: array([1., 1., 1.])}

clean_diagonals(reltol=1e-08)[source]

Eliminate essentially zerovalued diagonals

Parameters

reltol (number) – Relative tolerance

diags(format='dia')[source]

Return a regular sparse matrix of specified format

Parameters

format (str, optional) – Choice of matrix type (see scipy.sparse.diags)

• dia - Sparse matrix with DIAgonal storage

• csr - Compressed sparse row

Note

This method returns the matrix scaled by self.scale.

get_key()[source]
isdiagonal()[source]
isidentity()[source]
matvec(v, c, format='dia', axis=0)[source]

Matrix vector product

Returns c = dot(self, v)

Parameters
• v (array) – Numpy input array of ndim>=1

• c (array) – Numpy output array of same shape as v

• format (str, optional) – Choice for computation

• csr - Compressed sparse row format

• dia - Sparse matrix with DIAgonal storage

• python - Use numpy and vectorization

• self - To be implemented in subclass

• cython - Cython implementation that may be implemented in subclass

• axis (int, optional) – The axis over which to take the matrix vector product

same_keys(a)[source]
scale_array(c)[source]
solve(b, u=None, axis=0)[source]

Solve matrix system Au = b

where A is the current matrix (self)

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multi- dimensional

Note

Vectors may be one- or multidimensional.

class shenfun.matrixbase.SpectralMatrix(d, test, trial, scale=1.0, measure=1)[source]

Base class for inner product matrices.

Parameters
• d (dict) – Dictionary, where keys are the diagonal offsets and values the diagonals

• trial (2-tuple of (basis, int)) – The basis is an instance of a class for one of the bases in

The int represents the number of times the trial function should be differentiated. Representing matrix column.

• test (2-tuple of (basis, int)) – As trial, but representing matrix row.

• scale (number, optional) – Scale matrix with this number

Examples

Mass matrix for Chebyshev Dirichlet basis:

$(\phi_k, \phi_j)_w = \int_{-1}^{1} \phi_k(x) \phi_j(x) w(x) dx$

Stiffness matrix for Chebyshev Dirichlet basis:

$(\phi_k'', \phi_j)_w = \int_{-1}^{1} \phi_k''(x) \phi_j(x) w(x) dx$

The matrices can be automatically created using, e.g., for the mass matrix of the Dirichlet space:

SD = ShenDirichlet
N = 16
M = SpectralMatrix({}, (SD(N), 0), (SD(N), 0))


where the first (SD(N), 0) represents the test function and the second the trial function. The stiffness matrix can be obtained as:

A = SpectralMatrix({}, (SD(N), 0), (SD(N), 2))


where (SD(N), 2) signals that we use the second derivative of this trial function. The number N is the number of quadrature points used for the basis.

The automatically created matrices may be overloaded with more exactly computed diagonals.

Note that matrices with the Neumann basis are stored using index space $$k = 0, 1, ..., N-2$$, i.e., including the zero index for a nonzero average value.

property axis

Return the axis of the TensorProductSpace this matrix is created for

get_key()[source]
matvec(v, c, format='csr', axis=0)[source]

Matrix vector product

Returns c = dot(self, v)

Parameters
• v (array) – Numpy input array of ndim>=1

• c (array) – Numpy output array of same shape as v

• format (str, optional) – Choice for computation

• csr - Compressed sparse row format

• dia - Sparse matrix with DIAgonal storage

• python - Use numpy and vectorization

• self - To be implemented in subclass

• cython - Cython implementation that may be implemented in subclass

• axis (int, optional) – The axis over which to take the matrix vector product

simplify_diagonal_matrices()[source]
solve(b, u=None, axis=0)[source]

Solve matrix system Au = b

where A is the current matrix (self)

Parameters
• b (array) – Array of right hand side on entry and solution on exit unless u is provided.

• u (array, optional) – Output array

• axis (int, optional) – The axis over which to solve for if b and u are multidimensional

Note

Vectors may be one- or multidimensional.

property tensorproductspace

Return the TensorProductSpace this matrix has been computed for

class shenfun.matrixbase.TPMatrix(mats, testspace, trialspace, scale=1.0, global_index=None, mixedbase=None)[source]

Bases: object

Tensorproduct matrix

A TensorProductSpace is the outer product of D bases. A matrix assembled from test and trialfunctions on TensorProductSpaces will, as such, be represented as the outer product of D smaller matrices, one for each basis. This class represents the complete matrix.

Parameters
property dimensions

Return dimension of TPMatrix

get_key()[source]

Return key of the one nondiagonal matrix in the TPMatrix

Note

Raises an error of there are more than one single nondiagonal matrix in TPMatrix.

is_bc_matrix()[source]
isdiagonal()[source]
isidentity()[source]
matvec(v, c)[source]
simplify_diagonal_matrices()[source]
solve(b, u=None)[source]
shenfun.matrixbase.check_sanity(A, test, trial, measure=1)[source]

Sanity check for matrix.

Test that automatically created matrix agrees with overloaded one

Parameters
• A (matrix)

• test (2-tuple of (basis, int)) – The basis is an instance of a class for one of the bases in

The int represents the number of times the test function should be differentiated. Representing matrix row.

• trial (2-tuple of (basis, int)) – As test, but representing matrix column.

• measure (sympy function of coordinate, optional)

shenfun.matrixbase.extract_diagonal_matrix(M, abstol=1e-08, reltol=1e-08)[source]

Return SparseMatrix version of dense matrix M

Parameters
• M (Numpy array of ndim=2)

• abstol (float) – Tolerance. Only diagonals with max($$|d|$$) < tol are kept in the returned SparseMatrix, where $$d$$ is the diagonal

• reltol (float) – Relative tolerance. Only diagonals with max($$|d|$$)/max($$|M|$$) > reltol are kept in the returned SparseMatrix

shenfun.matrixbase.get_dense_matrix(test, trial, measure=1)[source]

Return dense matrix automatically computed from basis

Parameters
• test (2-tuple of (basis, int)) – The basis is an instance of a class for one of the bases in

The int represents the number of times the test function should be differentiated. Representing matrix row.

• trial (2-tuple of (basis, int)) – As test, but representing matrix column.

• measure (Sympy expression of coordinate, or number, optional) – Additional weight to integral. For example, in cylindrical coordinates an additional measure is the readius r.

shenfun.matrixbase.get_dense_matrix_sympy(test, trial, measure=1)[source]

Return dense matrix automatically computed from basis

Parameters
• test (2-tuple of (basis, int)) – The basis is an instance of a class for one of the bases in

The int represents the number of times the test function should be differentiated. Representing matrix row.

• trial (2-tuple of (basis, int)) – As test, but representing matrix column.

• measure (Sympy expression of coordinate, or number, optional) – Additional weight to integral. For example, in cylindrical coordinates an additional measure is the radius r.

## shenfun.spectralbase module¶

This module contains classes for working with the spectral-Galerkin method

There are currently classes for 11 bases and corresponding function spaces

All bases have expansions

$$u(x) = \sum_{k\in\mathcal{I}}\hat{u}_{k} \phi_k(x)$$

where $$\mathcal{I}$$ the index set of the basis, and the index set differs from base to base, see function space definitions below. $$\phi_k$$ is the k’t basis function in the basis. It is also called a test function, whereas $$u(x)$$ often is called a trial function.

Chebyshev:
Orthogonal:

basis functions: $$\phi_k = T_k$$

basis: $$span(T_k, k=0,1,..., N-1)$$

ShenDirichlet:
basis functions:

$$\phi_k = T_k-T_{k+2}$$

$$\phi_{N-2} = 0.5(T_0+T_1)$$ for Poisson’s equation

$$\phi_{N-1} = 0.5(T_0-T_1)$$ for Poisson’s equation

space: $$span(\phi_k, k=0,1,...,N-1)$$

where $$u(1)=a, u(-1)=b$$, such that $$\hat{u}_{N-2}=a, \hat{u}_{N-1}=b$$.

Note that there are only N-2 unknown coefficients, $$\hat{u}_k$$, since $$\hat{u}_{N-2}$$ and $$\hat{u}_{N-1}$$ are determined by boundary conditions. Inhomogeneous boundary conditions are possible for the Poisson equation, because $$\phi_{N-1}$$ and $$\phi_{N-2}$$ are in the kernel of the Poisson operator. For homogeneous boundary conditions $$\phi_{N-2}$$ and $$\phi_{N-1}$$ are simply ignored.

ShenNeumann:
basis function:

$$\phi_k = T_k-\left(\frac{k}{k+2}\right)^2T_{k+2}$$

space:

$$span(\phi_k, k=1,2,...,N-3)$$

Homogeneous Neumann boundary conditions, $$u'(\pm 1) = 0$$, and zero weighted mean: $$\int_{-1}^{1}u(x)w(x)dx = 0$$.

ShenBiharmonic:
basis function:

$$\phi_k = T_k - \frac{2(k+2)}{k+3}T_{k+2} + \frac{k+1}{k+3}T_{k+4}$$

space:

$$span(\phi_k, k=0,1,...,N-5)$$

Homogeneous Dirichlet and Neumann, $$u(\pm 1)=0$$ and $$u'(\pm 1)=0$$.

Legendre:
Orthogonal:
basis function:

$$\phi_k = L_k$$

space:

$$span(L_k, k=0,1,...N-1)$$

ShenDirichlet:
basis function:

$$\phi_k = L_k-L_{k+2}$$

$$\phi_{N-2} = 0.5(L_0+L_1)$$, for Poisson’s equation

$$\phi_{N-1} = 0.5(L_0-L_1)$$, for Poisson’s equation

space:

$$span(\phi_k, k=0,1,...,N-1)$$

where $$u(1)=a, u(-1)=b$$, such that $$\hat{u}_{N-2}=a, \hat{u}_{N-1}=b$$

Note that there are only N-2 unknown coefficients, $$\hat{u}_k$$, since $$\hat{u}_{N-2}$$ and $$\hat{u}_{N-1}$$ are determined by boundary conditions. Inhomogeneous boundary conditions are possible for the Poisson equation, because $$\phi_{N-1}$$ and $$\phi_{N-2}$$ are in the kernel of the Poisson operator. For homogeneous boundary conditions $$\phi_{N-2}$$ and $$\phi_{N-1}$$ are simply ignored.

ShenNeumann:
basis function:

$$\phi_k = L_k-\frac{k(k+1)}{(k+2)(k+3)}L_{k+2}$$

space:

$$span(\phi_k, k=1,2,...,N-3)$$

Homogeneous Neumann boundary conditions, $$u'(\pm 1) = 0$$, and zero mean: $$\int_{-1}^{1}u(x)dx = 0$$.

ShenBiharmonic:
basis function:

$$\phi_k = L_k - \frac{2(2k+5)}{2k+7}L_{k+2} + \frac{2k+3}{2k+7}L_{k+4}$$

space:

$$span(\phi_k, k=0,1,...,N-5)$$

Homogeneous Dirichlet and Neumann, $$u(\pm 1)=0$$ and $$u'(\pm 1)=0$$.

Laguerre:
Orthogonal:
basis function:

$$\phi_k(x) = L_k(x) \cdot \exp(-x)$$

space:

$$span(L_k, k=0,1,...N-1)$$

where $$L_k$$ is the Laguerre polynomial of order k.

ShenDirichlet:
basis function:

$$\phi_k = (L_k-L_{k+1})\cdot \exp(-x)$$

space:

$$span(\phi_k, k=0,1,...,N-2)$$

Homogeneous Dirichlet for domain [0, inf).

Hermite:
Orthogonal:
basis function:

$$\phi_k(x) = H_k(x) \cdot \exp(-x^2/2)/(\pi^{0.25}\sqrt{2^k k!})$$

space:

$$span(\phi_k, k=0,1,...N-1)$$

where $$H_k$$ is the Hermite polynomial of order k.

Jacobi:
Orthogonal:
basis function:

$$\phi_k = J_k(\alpha, \beta)$$

space:

$$span(L_k, k=0,1,...N-1)$$

where $$J_k(\alpha, \beta)$$ is the regular Jacobi polynomial and $$\alpha > -1$$ and $$\beta > -1$$.

ShenDirichlet:
basis function:

$$\phi_k = j_k(-1, -1)$$

space:

$$span(\phi_k, k=0,1,...,N-3)$$

where $$j_k$$ is the generalized Jacobi polynomial

ShenBiharmonic:
basis function:

$$\phi_k = j_k(-2, -2)$$

space:

$$span(\phi_k, k=0,1,...,N-5)$$

Homogeneous Dirichlet and Neumann, $$u(\pm 1)=0$$ and $$u'(\pm 1)=0$$.

ShenOrder6:
basis function:

$$\phi_k = j_k(-3, -3)$$

space:

$$span(\phi_k, k=0,1,...,N-7)$$

Homogeneous $$u(\pm 1)=u'(\pm 1)=u''(\pm 1)=0$$.

Fourier:
R2C and C2C:
basis function:

$$\phi_k = exp(ikx)$$

space:

$$span(\phi_k, k=-N/2, -N/2+1, ..., N/2-1)$$

Note that if N is even, then the Nyquist frequency (-N/2) requires special attention for the R2C space. We should then have been using a Fourier interpolator that symmetrizes by adding and extra $$\phi_{N/2}$$ and by multiplying $$\phi_{N/2}$$ and $$\phi_{-N/2}$$ with 0.5. This effectively sets the Nyquist frequency (-N/2) to zero for odd derivatives. This is not done automatically by shenfun. Instead we recommend using the SpectralBase.mask_nyquist() function that effectively sets the Nyquist frequency to zero (if N is even).

R2C and C2C are the same, but R2C is used on real physical data and it takes advantage of Hermitian symmetry, $$\hat{u}_{-k} = conj(\hat{u}_k)$$, for $$k = 1, ..., N/2$$

Each class has methods for moving (fast) between spectral and physical space, and for computing the (weighted) scalar product.

class shenfun.spectralbase.FuncWrap(func, input_array, output_array)[source]

Bases: object

property func
property input_array
property output_array
class shenfun.spectralbase.MixedFunctionSpace(bases)[source]

Bases: object

Class for composite bases in 1D

Parameters

spaces (list) – List of bases

dim()[source]

Return dimension of self (degrees of freedom)

property dimensions

Return dimension of scalar space

dims()[source]

Return dimensions (degrees of freedom) for MixedFunctionSpace

flatten()[source]
property is_composite_space
num_components()[source]

Return number of bases in mixed basis

property rank
shape(forward_output=False)[source]

Return shape of arrays for MixedFunctionSpace

Parameters

forward_output (bool, optional) – If True then return shape of an array that is the result of a forward transform. If False then return shape of physical space, i.e., the input to a forward transform.

class shenfun.spectralbase.SpectralBase(N, quad='', padding_factor=1, domain=- 1.0, 1.0, dtype=None, dealias_direct=False, coordinates=None)[source]

Bases: object

Abstract base class for all spectral function spaces

Parameters
• N (int) – Number of quadrature points

• GL - Chebyshev-Gauss-Lobatto or Legendre-Gauss-Lobatto

• GC - Chebyshev-Gauss

• LG - Legendre-Gauss or Laguerre-Gauss

• HG - Hermite-Gauss

• domain (2-tuple of floats, optional) – The computational domain

• dtype (data-type, optional) – Type of input data in real physical space. Will be overloaded when basis is part of a TensorProductSpace.

• dealias_direct (bool, optional) – If True, then set all upper 2/3 wavenumbers to zero before backward transform. Cannot be used if padding_factor is different than 1.

• coordinates (2- or 3-tuple (coordinate, position vector (, sympy assumptions)), optional) – Map for curvilinear coordinatesystem. The new coordinate variable in the new coordinate system is the first item. Second item is a tuple for the Cartesian position vector as function of the new variable in the first tuple. Example:

theta = sp.Symbols('x', real=True, positive=True)
rv = (sp.cos(theta), sp.sin(theta))


where theta and rv are the first and second items in the tuple. If a third item is provided with the tuple, then this third item is used as an additional assumption. For example, it is necessary to provide the assumption sympy.Q.positive(sympy.sin(theta)), such that sympy can evaluate sqrt(sympy.sin(theta)**2) to sympy.sin(theta) and not Abs(sympy.sin(theta)). Different coordinates may require different assumptions to help sympy when computing basis functions etc.

apply_inverse_mass(array)[source]

Apply inverse mass matrix

Parameters

array (array (input/output)) – Expansion coefficients. Overwritten by applying the inverse mass matrix, and returned.

Returns

Return type

array

backward(input_array=None, output_array=None, fast_transform=True)[source]

Compute backward (inverse) transform

Parameters
• input_array (array, optional) – Expansion coefficients

• output_array (array, optional) – Function values on quadrature mesh

• fast_transform (bool, optional) – If True use fast transforms (if implemented), if False use Vandermonde type

Note

If input_array/output_array are not given, then use predefined arrays as planned with self.plan

backward_uniform(input_array=None, output_array=None)[source]

Evaluate function on uniform mesh

Parameters
• input_array (array, optional) – Expansion coefficients

• output_array (array, optional) – Function values on quadrature mesh

Note

If input_array/output_array are not given, then use predefined arrays as planned with self.plan

static boundary_condition()[source]
broadcast_to_ndims(x)[source]

Return 1D array x as an array of shape according to the dimensions() of the TensorProductSpace class that this base (self) belongs to.

Parameters

x (1D array)

Note

The returned array has shape one in all ndims-1 dimensions apart from self.axis.

Example

>>> import numpy as np
>>> from shenfun import Basis, TensorProductSpace
>>> from mpi4py import MPI
>>> K0 = FunctionSpace(8, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> T = TensorProductSpace(MPI.COMM_WORLD, (K0, K1))
>>> x = np.arange(4)
>>> print(y.shape)
(4, 1)

compatible_base(space)[source]
curvilinear_mesh(uniform=False)[source]

Return curvilinear mesh of basis

Parameters

uniform (bool, optional) – Use uniform mesh

dim()[source]

Return the dimension of self (the number of degrees of freedom)

property dimensions

Return the dimensions (the number of bases) of the TensorProductSpace class this basis is planned for.

dims()[source]

Return tuple (length one since a basis only has one dim) containing self.dim()

domain_factor()[source]

Return scaling factor for domain

property dtype

Return datatype basis is planned for

eval(x, u, output_array=None)[source]

Evaluate Function u at position x

Parameters
• x (float or array of floats)

• u (array) – Expansion coefficients or instance of Function

• output_array (array, optional) – Function values at points

Returns

output_array

Return type

array

evaluate_basis(x, i=0, output_array=None)[source]

Evaluate basis i at points x

Parameters
• x (float or array of floats)

• i (int, optional) – Basis number

• output_array (array, optional) – Return result in output_array if provided

Returns

output_array

Return type

array

evaluate_basis_all(x=None, argument=0)[source]

Evaluate basis at x or all quadrature points

Parameters
• x (float or array of floats, optional) – If not provided use quadrature points of self

• argument (int) – Zero for test and 1 for trialfunction

Returns

Vandermonde matrix

Return type

array

evaluate_basis_derivative(x=None, i=0, k=0, output_array=None)[source]

Evaluate k’th derivative of basis i at x or all quadrature points

Parameters
• x (float or array of floats, optional) – If not provided use quadrature points of self

• i (int, optional) – Basis number

• k (int, optional) – k’th derivative

• output_array (array, optional) – return array

Returns

output_array

Return type

array

evaluate_basis_derivative_all(x=None, k=0, argument=0)[source]

Return k’th derivative of basis evaluated at x or all quadrature points as a Vandermonde matrix.

Parameters
• x (float or array of floats, optional) – If not provided use quadrature points of self

• k (int, optional) – k’th derivative

• argument (int) – Zero for test and 1 for trialfunction

Returns

Vandermonde matrix

Return type

array

evaluate_expansion_all(input_array, output_array, fast_transform=False)[source]

Evaluate expansion on entire mesh

$u(x_j) = \sum_{k\in\mathcal{I}} \hat{u}_k T_k(x_j) \quad \text{ for all} \quad j = 0, 1, ..., N$
Parameters
• input_array ($$\hat{u}_k$$) – Expansion coefficients, or instance of Function

• output_array ($$u(x_j)$$) – Function values on quadrature mesh, instance of Array

• fast_transform (bool, optional) – Whether to use fast transforms (if implemented)

evaluate_scalar_product(input_array, output_array, fast_transform=False)[source]

Evaluate scalar product

Parameters
• input_array (array, optional) – Function values on quadrature mesh

• output_array (array, optional) – Expansion coefficients

• fast_transform (bool, optional) – If True use fast transforms (if implemented), if False use Vandermonde type

static family()[source]
forward(input_array=None, output_array=None, fast_transform=True)[source]

Compute forward transform

Parameters
• input_array (array, optional) – Function values on quadrature mesh

• output_array (array, optional) – Expansion coefficients

• fast_transform (bool, optional) – If True use fast transforms, if False use Vandermonde type

Note

If input_array/output_array are not given, then use predefined arrays as planned with self.plan

get_bcmass_matrix(dx=1)[source]
get_dealiased(padding_factor=1.5, dealias_direct=False)[source]

Return space (otherwise as self) to be used for dealiasing

Parameters

• dealias_direct (bool, optional) – If True, set upper 2/3 of wavenumbers to zero before backward transform. Cannot be used together with padding_factor different than 1.

Returns

The space to be used for dealiasing

Return type

SpectralBase

get_mass_matrix()[source]
get_measured_array(array)[source]

Return weights times measure

Parameters
• N (integer, optional) – The number of quadrature points

• measure (None or sympy.Expr)

Note

If basis is part of a TensorProductSpace, then the array will be measured there. So in that case, just return the array unchanged.

get_measured_weights(N=None, measure=1)[source]

Return weights times measure

Parameters
• N (integer, optional) – The number of quadrature points

• measure (None or sympy.Expr)

get_normalization()[source]
get_refined(N)[source]

Return space (otherwise as self) with N quadrature points

Parameters

N (int) – The number of quadrature points for returned space

Returns

A new space with new number of quadrature points, otherwise as self.

Return type

SpectralBase

property has_nonhomogeneous_bcs
property is_composite_space
property is_orthogonal
map_reference_domain(x)[source]

Return true point x mapped to reference domain

map_true_domain(x)[source]

Return reference point x mapped to true domain

mesh(bcast=True, map_true_domain=True, uniform=False)[source]

Return the computational mesh

Parameters
mpmath_points_and_weights(N=None, map_true_domain=False, weighted=True, **kw)[source]

Return points and weights of quadrature using extended precision mpmath if available

$\int_{\Omega} f(x) w(x) dx \approx \sum_{i} f(x_i) w_i$
Parameters
• N (int, optional) – Number of quadrature points

• map_true_domain (bool, optional) – Whether or not to map points to true domain

• weighted (bool, optional) – Whether to use quadrature weights for a weighted inner product (default), or a regular, non-weighted inner product.

Note

If not implemented, or if mpmath/quadpy are not available, then simply returns the regular numpy points_and_weights().

property ndim

Return ndim of basis

num_components()[source]

Return number of components for basis

plan(shape, axis, dtype, options)[source]

Plan transform

Allocate work arrays for transforms and set up methods forward, backward and scalar_product with or without padding

Parameters
• shape (array) – Local shape of global array

• axis (int) – This base’s axis in global TensorProductSpace

• dtype (numpy.dtype) – Type of array

• options (dict) – Options for planning transforms

points_and_weights(N=None, map_true_domain=False, weighted=True, **kw)[source]

Return points and weights of quadrature for weighted integral

$\int_{\Omega} f(x) w(x) dx \approx \sum_{i} f(x_i) w_i$
Parameters
• N (int, optional) – Number of quadrature points

• map_true_domain (bool, optional) – Whether or not to map points to true domain

• weighted (bool, optional) – Whether to use quadrature weights for a weighted inner product (default), or a regular, non-weighted inner product.

property rank

Return tensor rank of basis

reference_domain()[source]

Return reference domain of basis

scalar_product(input_array=None, output_array=None, fast_transform=True)[source]

Compute weighted scalar product

Parameters
• input_array (array, optional) – Function values on quadrature mesh

• output_array (array, optional) – Expansion coefficients

• fast_transform (bool, optional) – If True use fast transforms, if False use Vandermonde type

Note

If input_array/output_array are not given, then use predefined arrays as planned with self.plan

shape(forward_output=True)[source]

Return the allocated shape of arrays used for self

Parameters

forward_output (bool, optional) – If True then return allocated shape of spectral space (the result of a forward transform). If False then return allocated shape of physical space (the input to a forward transform).

slice()[source]

Return index set of current basis

sympy_basis(i=0, x=None)[source]

Return basis function i as sympy function

Parameters
• i (int, optional) – The degree of freedom of the basis function

• x (sympy Symbol, optional)

sympy_basis_all(x=None)[source]

Return all basis functions as sympy functions

sympy_reference_domain()[source]
sympy_weight(x=None)[source]
property tensorproductspace

Return the last TensorProductSpace this basis has been planned for (if planned)

Note

A basis may be part of several :class:.TensorProductSpaces, but they all need to be of the same global shape.

to_ortho(input_array, output_array=None)[source]

Project to orthogonal basis

Parameters
• input_array (array) – Expansion coefficients of input basis

• output_array (array, optional) – Expansion coefficients in orthogonal basis

Returns

output_array

Return type

array

vandermonde(x)[source]

Return Vandermonde matrix based on the primary basis of the family.

Evaluates basis $$\psi_k(x)$$ for all wavenumbers, and all x. Returned Vandermonde matrix is an N x M matrix with N the length of x and M the number of bases.

$\begin{split}\begin{bmatrix} \psi_0(x_0) & \psi_1(x_0) & \ldots & \psi_{M-1}(x_0)\\ \psi_0(x_1) & \psi_1(x_1) & \ldots & \psi_{M-1}(x_1)\\ \vdots & \ldots \\ \psi_{0}(x_{N-1}) & \psi_1(x_{N-1}) & \ldots & \psi_{M-1}(x_{N-1}) \end{bmatrix}\end{split}$
Parameters

x (array of floats) – points for evaluation

Note

This function returns a matrix of evaluated primary basis functions for either family. That is, it is using either pure Chebyshev, Legendre or Fourier exponentials. The true Vandermonde matrix of a basis is obtained through SpectralBase.evaluate_basis_all().

vandermonde_evaluate_expansion(points, input_array, output_array)[source]

Evaluate expansion at certain points, possibly different from the quadrature points

Parameters
• points (array) – Points for evaluation

• input_array (array) – Expansion coefficients

• output_array (array) – Function values on points

vandermonde_evaluate_expansion_all(input_array, output_array, x=None)[source]

Naive implementation of evaluate_expansion_all

Parameters
• input_array (array) – Expansion coefficients

• output_array (array) – Function values on quadrature mesh

• x (mesh or None, optional)

vandermonde_scalar_product(input_array, output_array)[source]

Naive implementation of scalar product

Parameters
• input_array (array) – Function values on quadrature mesh

• output_array (array) – Expansion coefficients

wavenumbers(bcast=True, **kw)[source]

Return the wavenumbermesh

Parameters

bcast (bool) – Whether or not to broadcast to SpectralBase.dimensions() if basis belongs to a TensorProductSpace

class shenfun.spectralbase.Transform(func, xfftn, input_array, tmp_array, output_array)[source]
property tmp_array
property xfftn
class shenfun.spectralbase.VectorBasisTransform(transforms)[source]

Bases: object

shenfun.spectralbase.inner_product(test, trial, measure=1)[source]

Return 1D weighted inner product of bilinear form

Parameters
• test (2-tuple of (Basis, integer)) – Basis is any of the classes from

The integer determines the number of times the basis is differentiated. The test represents the matrix row

• trial (2-tuple of (Basis, integer)) – Like test

• measure (function of coordinate, optional)

Note

This function only performs 1D inner products and is unaware of any TensorProductSpace

Example

Compute mass matrix of Shen’s Chebyshev Dirichlet basis:

>>> from shenfun.spectralbase import inner_product
>>> from shenfun.chebyshev.bases import ShenDirichletBasis
>>> SD = ShenDirichletBasis(6)
>>> B = inner_product((SD, 0), (SD, 0))
>>> d = {-2: np.array([-np.pi/2]),
...       0: np.array([1.5*np.pi, np.pi, np.pi, np.pi]),
...       2: np.array([-np.pi/2])}
>>> [np.all(B[k] == v) for k, v in d.items()]
[True, True, True]

class shenfun.spectralbase.islicedict(axis=0, dimensions=1)[source]

Bases: dict

Return a tuple of slices, broadcasted to dimensions number of dimensions, and with integer a along axis.

Parameters

Example

>>> from shenfun.spectralbase import islicedict
>>> s = islicedict(axis=1, dimensions=3)
>>> print(s[0])
(slice(None, None, None), 0, slice(None, None, None))

class shenfun.spectralbase.slicedict(axis=0, dimensions=1)[source]

Bases: dict

Return a tuple of slices, broadcasted to dimensions number of dimensions, and with slice a along axis.

Parameters

Example

>>> from shenfun.spectralbase import slicedict
>>> s = slicedict(axis=1, dimensions=3)
>>> print(s[slice(0, 5)])
(slice(None, None, None), slice(0, 5, None), slice(None, None, None))


## shenfun.tensorproductspace module¶

Module for implementation of the TensorProductSpace class and related methods.

class shenfun.tensorproductspace.Convolve(padding_space)[source]

Bases: object

Class for convolving without truncation.

The convolution of $$\hat{a}$$ and $$\hat{b}$$ is computed by first transforming backwards with padding:

a = Tp.backward(a_hat)
b = Tp.backward(b_hat)


and then transforming the product a*b forward without truncation:

ab_hat = T.forward(a*b)


where Tp is a TensorProductSpace for regular padding, and T is a TensorProductSpace with no padding, but using the shape of the padded a and b arrays.

For convolve with truncation forward, use just the convolve method of the Tp space instead.

Parameters

padding_space (TensorProductSpace) – Space with regular padding backward and truncation forward.

__call__(a_hat, b_hat, ab_hat=None)[source]

Compute convolution of a_hat and b_hat without truncation

Parameters
class shenfun.tensorproductspace.MixedTensorProductSpace(spaces)[source]

Bases: object

Class for composite tensorproductspaces.

The mixed spaces are Cartesian products of TensorproductSpaces or other MixedTensorProductSpaces.

Parameters

spaces (list) – List of spaces

compatible_base(space)[source]
convolve(a_hat, b_hat, ab_hat)[source]

Convolution of a_hat and b_hat

Parameters
• a_hat (array) – Input array of shape and type as output array from self.forward, or instance of Function

• b_hat (array) – Input array of shape and type as output array from self.forward, or instance of Function

• ab_hat (array) – Return array of same type and shape as a_hat and b_hat

Note

Note that self should have bases with padding for this method to give a convolution without aliasing. The padding is specified when creating instances of bases for the TensorProductSpace.

FIXME Efficiency due to allocation

dim()[source]

Return dimension of self (degrees of freedom)

property dimensions

Return dimension of scalar space

dims()[source]

Return dimensions (degrees of freedom) of all bases in self

eval(points, coefficients, output_array=None, method=0)[source]

Evaluate Function at points, given expansion coefficients

Parameters
• points (float or array of floats)

• coefficients (array) – Expansion coefficients

• output_array (array, optional) – Return array, function values at points

• method (int, optional) – Chooses implementation. The default 0 is a low-memory cython version. Using method = 1 leads to a faster cython implementation that, on the downside, uses more memory. The final, method = 2, is a python implementation used only for verification.

flatten()[source]
get_dealiased(padding_factor=1.5, dealias_direct=False)[source]
get_orthogonal()[source]
get_refined(N)[source]
global_shape(forward_output=False)[source]

Return global shape for MixedTensorProductSpace

Parameters

forward_output (bool, optional) – If True then return shape of an array that is the result of a forward transform. If False then return shape of physical space, i.e., the input to a forward transform.

property is_composite_space
local_slice(forward_output=False)[source]

The local view into the global data

Parameters

forward_output (bool, optional) – Return local slices of output array (spectral space) if True, else return local slices of input array (physical space)

num_components()[source]

Return number of spaces in mixed space

property rank

Return rank of space

shape(forward_output=False)[source]

Return shape of arrays for MixedTensorProductSpace

Parameters

forward_output (bool, optional) – If True then return shape of an array that is the result of a forward transform. If False then return shape of physical space, i.e., the input to a forward transform.

Note

A MixedTensorProductSpace may contain tensor product spaces of different shape in spectral space. Hence this function returns a list of shapes and not one single tuple.

size(forward_output=False)[source]

Return number of elements in MixedTensorProductSpace

slice()[source]

The slices of dofs for each dimension

class shenfun.tensorproductspace.TensorProductSpace(comm, bases, axes=None, dtype=None, slab=False, collapse_fourier=False, backward_from_pencil=False, coordinates=None, **kw)[source]

Class for multidimensional tensorproductspaces.

The tensorproductspaces are created as Cartesian products from a set of 1D bases. The 1D bases are subclassed instances of the SpectralBase class.

Parameters
• comm (MPI communicator)

• bases (list) – List of 1D bases

• axes (tuple of ints, optional) – A tuple containing the order of which to perform transforms. Last item is transformed first. Defaults to range(len(bases))

• dtype (data-type, optional) – Type of input data in real physical space. If not provided it will be inferred from the bases.

• slab (bool, optional) – Use 1D slab decomposition.

• collapse_fourier (bool, optional) – Collapse axes for Fourier bases if possible

• backward_from_pencil (False or Pencil) – In case of Pencil configure the transform by starting from the Pencil distribution in spectral space. This is primarily intended for a padded space, where the spectral distribution must be equal to the non-padded space.

• coordinates (2- or 3-tuple, optional) – Map for curvilinear coordinatesystem. First tuple are the coordinate variables in the new coordinate system Second tuple are the Cartesian coordinates as functions of the variables in the first tuple. Example:

psi = (theta, r) = sp.symbols('x,y', real=True, positive=True)
rv = (r*sp.cos(theta), r*sp.sin(theta))


where psi and rv are the first and second tuples, respectively. If a third item is provided with the tuple, then this third item is used as an additional assumption. For example, it is necessary to provide the assumption sympy.Q.positive(sympy.sin(theta)), such that sympy can evaluate sqrt(sympy.sin(theta)**2) to sympy.sin(theta) and not Abs(sympy.sin(theta)). Different coordinates may require different assumptions to help sympy when computing basis functions etc.

• kw (dict, optional) – Dictionary that can be used to plan transforms. Input to method plan for the bases.

compatible_base(space)[source]

Return whether space is compatible with self.

Parameters

space (TensorProductSpace) – The space compared to

Note

Two spaces are deemed compatible if the underlying bases, along each direction, belong to the same family, has the same number of quadrature points, and the same quadrature scheme.

configure_backwards(pencil, dtype, kw)[source]

Configure transforms starting from spectral space

Parameters
• pencil (Pencil) – The distribution in spectral space

• dtype (Numpy dtype) – The type of data in spectral space

• kw (dict) – Any parameters for planning transforms

Note

To ensure the same distribution in spectral space, the padded space must be configured by moving from the spectral space towards the physical space. The distribution does not have to agree in physical space, because the padding is done only in the spectral.

convolve(a_hat, b_hat, ab_hat)[source]

Convolution of a_hat and b_hat

Parameters
• a_hat (array) – Input array of shape and type as output array from self.forward, or instance of Function

• b_hat (array) – Input array of shape and type as output array from self.forward, or instance of Function

• ab_hat (array) – Return array of same type and shape as a_hat and b_hat

Note

The return array ab_hat is truncated to the shape of a_hat and b_hat.

Also note that self should have bases with padding for this method to give a convolution without aliasing. The padding is specified when creating instances of bases for the TensorProductSpace.

FIXME Efficiency due to allocation

curvilinear_mesh(uniform=False)[source]

Return curvilinear mesh

Parameters

uniform (bool, optional) – Use uniform mesh for non-periodic bases

dim()[source]

Return dimension of self (degrees of freedom)

property dimensions

Return dimension of TensorProductSpace

dims()[source]

Return dimensions (degrees of freedom) of all bases in self

eval(points, coefficients, output_array=None, method=2)[source]

Evaluate Function at points, given expansion coefficients

Parameters
• points (float or array of floats) – Array must be of shape (D, N), for N points in D dimensions

• coefficients (array) – Expansion coefficients, or instance of Function

• output_array (array, optional) – Return array, function values at points

• method (int, optional) – Chooses implementation. The default 0 is a low-memory cython version. Using method = 1 leads to a faster cython implementation that, on the downside, uses more memory. The final, method = 2, is a python implementation.

get_dealiased(padding_factor=1.5, dealias_direct=False)[source]
get_mask_nyquist()[source]

Return an array with zeros for Nyquist coefficients and one otherwise

get_measured_array(u)[source]

Weigh Array u with integral measure

Parameters

u (Array)

get_nondiagonal_axes()[source]

Return list of axes that may contain non-diagonal matrices

get_nonperiodic_axes()[source]

Return list of axes that are not periodic

get_orthogonal()[source]
get_refined(N)[source]
global_shape(forward_output=False)[source]

Return global shape of arrays for TensorProductSpace

Parameters

forward_output (bool, optional) – If True then return shape of an array that is the result of a forward transform. If False then return shape of physical space, i.e., the input to a forward transform.

property is_composite_space
property is_orthogonal
local_curvilinear_mesh(uniform=False)[source]

Return curvilinear mesh

Parameters

uniform (bool, optional) – Use uniform mesh for non-periodic bases

local_mesh(broadcast=False, uniform=False)[source]

Return list of local 1D physical meshes for each dimension of TensorProductSpace

Parameters
• broadcast (bool, optional) – Broadcast each 1D mesh to real shape of TensorProductSpace

• uniform (bool, optional) – Use uniform mesh for non-periodic bases

local_wavenumbers(broadcast=False, scaled=False, eliminate_highest_freq=False)[source]

Return list of local wavenumbers of TensorProductSpace

Parameters
• broadcast (bool, optional) – Broadcast returned wavenumber arrays to actual dimensions of TensorProductSpace

• scaled (bool, optional) – Scale wavenumbers with size of box

• eliminate_highest_freq (bool, optional) – Set Nyquist frequency to zero for evenly shaped axes of Fourier bases

mask_nyquist(u_hat, mask=None)[source]

Return Function u_hat with zero Nyquist coefficients

Parameters
• u_hat (array) – Function to be masked

• mask (array or None, optional) – mask array, if not provided then get the mask by calling get_mask_nyquist()

mesh(uniform=False)[source]

Return list of 1D physical meshes for each dimension of TensorProductSpace

Parameters

uniform (bool, optional) – Use uniform mesh for non-periodic bases

num_components()[source]

Return number of scalar spaces in TensorProductSpace

property rank

Return tensor rank of TensorProductSpace

size(forward_output=False)[source]

Return number of elements in TensorProductSpace

slice()[source]

The slices of dofs for each dimension

wavenumbers(scaled=False, eliminate_highest_freq=False)[source]

Return list of wavenumbers of TensorProductSpace

Parameters
• scaled (bool, optional) – Scale wavenumbers with size of box

• eliminate_highest_freq (bool, optional) – Set Nyquist frequency to zero for evenly shaped axes of Fourier bases.

class shenfun.tensorproductspace.VectorTensorProductSpace(space)[source]

A special MixedTensorProductSpace where the number of spaces must equal the geometrical dimension of the problem.

For example, a TensorProductSpace created by a tensorproduct of 2 1D bases, will have vectors of length 2. A TensorProductSpace created from 3 1D bases will have vectors of length 3.

Parameters

space (TensorProductSpace or list of ndim :class:.TensorProductSpaces) – Spaces to create vector from

get_dealiased(padding_factor=1.5, dealias_direct=False)[source]
get_orthogonal()[source]
get_refined(N)[source]
num_components()[source]

Return number of spaces in mixed space

property rank

Return tensor rank of space

shape(forward_output=False)[source]

Return shape of arrays for VectorTensorProductSpace

Parameters

forward_output (bool, optional) – If True then return shape of an array that is the result of a forward transform. If False then return shape of physical space, i.e., the input to a forward transform.

## Module contents¶

This is the shenfun package

### What is shenfun?¶

Shenfun is a high performance computing platform for solving partial differential equations (PDEs) by the spectral Galerkin method. The user interface to shenfun is very similar to FEniCS, but applications are limited to multidimensional tensor product grids. The code is parallelized with MPI through the mpi4py-fft package.

Shenfun enables fast development of efficient and accurate PDE solvers (spectral order and accuracy), in the comfortable high-level Python language. The spectral accuracy is ensured from using high-order global orthogonal basis functions (Fourier, Legendre, Chebyshev, Laguerre, Hermite and Jacobi), as opposed to finite element codes like FEniCS that are using low-order local basis functions. Efficiency is ensured through vectorization (Numpy), parallelization (mpi4py) and by moving critical routines to Cython.

Shenfun has support for solving scalar and vector equations in curvilinear coordinates, like polar or spherical coordinates.