shenfun package

Subpackages

Submodules

shenfun.la module

This module contains linear algebra solvers for SparseMatrices, TPMatrices and BlockMatrices.

class shenfun.la.BandedMatrixSolver(mats)[source]

Bases: SparseMatrixSolver

static LU(data)[source]

LU-decomposition using either Cython or Numba

Parameters:

data (2D-array) – Storage for dia-matrix on entry and L and U matrices on exit.

static Solve(u, data, axis)[source]

Fast solve using either Cython or Numba

Parameters:

  • u (array) – rhs on entry, solution on exit

  • data (2D-array) – Storage for dia-matrix containing L and U matrices

  • axis (int) – The axis we are solving over

solve(b, u, axis, lu)[source]

Solve Au=b

Solve along axis if b and u are multidimensional arrays.

Parameters:

  • b, u (arrays of rhs and output) – Both can be multidimensional

  • axis (int) – The axis we are solving over

  • lu (LU-decomposition) – Can be either the output from splu, or a dia-matrix containing the L and U matrices. The latter is used in subclasses.

class shenfun.la.BlockMatrixSolver(mats)[source]

Bases: object

__call__(b, u=None, constraints=())[source]

Call self as a function.

static apply_constraint(A, b, offset, i, constraint)[source]
class shenfun.la.DiagMA(mats)[source]

Bases: BandedMatrixSolver

Diagonal matrix solver

Parameters:

mats (Diagonal SparseMatrix or list of diagonal SparseMatrix instances)

Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, lu)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.FDMA(mats)[source]

Bases: BandedMatrixSolver

4-diagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – 4-diagonal matrix with diagonals in offsets -2, 0, 2, 4

LU = <shenfun.optimization.runtimeoptimizer object>
Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.HeptaDMA(mats)[source]

Bases: BandedMatrixSolver

Heptadiagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – Heptadiagonal matrix with diagonals in offsets -4, -2, 0, 2, 4, 6, 8

LU = <shenfun.optimization.runtimeoptimizer object>
Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.PDMA(mats)[source]

Bases: BandedMatrixSolver

Pentadiagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – Pentadiagonal matrix with diagonals in offsets -4, -2, 0, 2, 4

LU = <shenfun.optimization.runtimeoptimizer object>
Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.Solve(mats, format=None)[source]

Bases: SparseMatrixSolver

Generic solver class for :SparseMatrix

Possibly with inhomogeneous boundary values

Parameters:

  • mats (SparseMatrix or list of SparseMatrix instances)

  • format (str, optional) – The format of the scipy.sparse.spmatrix to convert into before solving. Default is Compressed Sparse Column csc.

Note

This solver converts the matrix to a Scipy sparse matrix of choice and uses scipy.sparse methods splu and spsolve.

shenfun.la.Solver(mats)[source]

Return appropriate solver for mats

Parameters:

mats (SparseMatrix or list of SparseMatrix instances)

Return type:

Matrix solver (SparseMatrixSolver)

Note

The list of matrices may include boundary matrices. The returned solver will incorporate these boundary matrices automatically on the right hand side of the equation system.

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

Bases: object

Generic solver for tensorproductspaces in 2D

Parameters:

mats (sequence) – sequence of instances of TPMatrix

Note

If there are boundary matrices in the list of mats, then these matrices are used to modify the right hand side before solving. If this is not the desired behaviour, then use extract_bc_matrices() on tpmats before using this class.

__call__(b, u=None, constraints=())[source]

Solve generic problem for sum of TPMatrix instances

Parameters:

  • b (array, right hand side)

  • u (array, solution)

  • constraints (tuple of 2-tuples) – Each 2-tuple (row, value) is a constraint set for the non-periodic direction, for Fourier index 0 in 2D and (0, 0) in 3D

static apply_constraints(A, b, constraints)[source]
matvec(u, c)[source]
class shenfun.la.Solver3D(tpmats)[source]

Bases: Solver2D

Generic solver for tensorproductspaces in 3D

Parameters:

mats (sequence) – sequence of instances of TPMatrix

Note

If there are boundary matrices in the list of mats, then these matrices are used to modify the right hand side before solving. If this is not the desired behaviour, then use extract_bc_matrices() on mats before using this class.

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

Bases: object

Solver for purely diagonal matrices, like Fourier in Cartesian coordinates.

Parameters:

tpmats (sequence) – sequence of instances of TPMatrix

__call__(b, u=None, constraints=())[source]

Solve problem with TPMatrix consisting only of diagonal matrices

Parameters:

  • b (array, right hand side)

  • u (array, solution)

  • constraints (tuple of 2-tuples) – Each 2-tuple (row, value) is a constraint set for the non-periodic direction, for Fourier index 0 in 2D and (0, 0) in 3D

class shenfun.la.SolverGeneric1ND(mats)[source]

Bases: object

Generic solver for tensorproduct matrices consisting of non-diagonal matrices along only one axis and Fourier along the others.

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 note that if there are boundary matrices in the list of mats, then these matrices are used to modify the right hand side before solving. If this is not the desired behaviour, then use extract_bc_matrices() on mats before using this class.

__call__(b, u=None, constraints=(), fast=True)[source]

Solve problem with one non-diagonal direction

Parameters:

  • b (array, right hand side)

  • u (array, solution)

  • constraints (tuple of 2-tuples) – Each 2-tuple (row, value) is a constraint set for the non-periodic direction, for Fourier index 0 in 2D and (0, 0) in 3D

  • fast (bool) – Use fast routine if possible. A fast routine is possible for any system of matrices with a tailored solver, like the TDMA, PDMA, FDMA and TwoDMA.

apply_constraints(b, constraints=())[source]
assemble()[source]
fast_solve(u, b, solvers1D, naxes)[source]
get_data(is_rank_zero)[source]
matvec(u, c)[source]
perform_lu()[source]
solve(u, b, solvers1D, naxes)[source]
solve_data = <shenfun.optimization.runtimeoptimizer object>
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, constraints=())[source]

Solve generic problem

Parameters:

  • b (array, right hand side)

  • u (array, solution)

  • constraints (tuple of 2-tuples) – Each 2-tuple (row, value) is a constraint set for the non-periodic direction, for Fourier index 0 in 2D and (0, 0) in 3D

apply_constraints(b, constraints)[source]
assemble()[source]
diags(i)[source]

Return matrix for given index i in diagonal direction

Parameters:

i (int) – Fourier wavenumber

get_diagonal_axis()[source]
matvec(u, c)[source]
perform_lu()[source]
class shenfun.la.SolverND(tpmats)[source]

Bases: Solver2D

Generic solver for tensorproductspaces in N dimensions

Parameters:

mats (sequence) – sequence of instances of TPMatrix

Note

If there are boundary matrices in the list of mats, then these matrices are used to modify the right hand side before solving. If this is not the desired behaviour, then use extract_bc_matrices() on mats before using this class.

class shenfun.la.SparseMatrixSolver(mats)[source]

Bases: object

Solver for SparseMatrix matrices.

Parameters:

mats (SparseMatrix or list of SparseMatrix instances)

Note

The list of matrices may include boundary matrices. The returned solver will incorporate these boundary matrices automatically on the right hand side of the equation system.

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

Solve matrix problem Au = b along axis

This routine also applies boundary conditions and constraints, and performes LU-decomposition on the fully assembled matrix.

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

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

Note

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

apply_bcs(b, u, axis=0)[source]
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, lu)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
solve(b, u, axis, lu)[source]

Solve Au=b

Solve along axis if b and u are multidimensional arrays.

Parameters:

  • b, u (arrays of rhs and output) – Both can be multidimensional

  • axis (int) – The axis we are solving over

  • lu (LU-decomposition) – Can be either the output from splu, or a dia-matrix containing the L and U matrices. The latter is used in subclasses.

class shenfun.la.TDMA(mats)[source]

Bases: BandedMatrixSolver

Tridiagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – Tridiagonal matrix with diagonals in offsets -2, 0, 2

LU = <shenfun.optimization.runtimeoptimizer object>
Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.TDMA_O(mats)[source]

Bases: BandedMatrixSolver

Tridiagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – Symmetric tridiagonal matrix with diagonals in offsets -1, 0, 1

LU = <shenfun.optimization.runtimeoptimizer object>
Solve = <shenfun.optimization.runtimeoptimizer object>
static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.ThreeDMA(mats)[source]

Bases: BandedMatrixSolver

3-diagonal matrix solver - all diagonals upper

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – 3-diagonal matrix with diagonals in offsets 0, 2, 4

Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]
class shenfun.la.TwoDMA(mats)[source]

Bases: BandedMatrixSolver

2-diagonal matrix solver

Parameters:

mats (SparseMatrix or list of SparseMatrix instances) – 2-diagonal matrix with diagonals in offsets 0, 2

Solve = <shenfun.optimization.runtimeoptimizer object>
apply_constraints(b, constraints, axis=0)[source]

Apply constraints to matrix self.mat and rhs vector b

Parameters:

  • b (array)

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

  • axis (int) – The axis we are solving over

static inner_solve(u, data)[source]

Solve Au=b for one-dimensional u

On entry u is the rhs b, on exit it contains the solution.

Parameters:

  • u (array 1D) – rhs on entry and solution on exit

  • lu (LU-decomposition) – Can be either a 2-tuple with (output from splu, dtype), or a scipy dia-matrix containing the L and U matrices. The latter is used in subclasses.

perform_lu()[source]

shenfun.matrixbase module

This module contains classes for working with sparse matrices.

shenfun.matrixbase.BlockMatrices(tpmats)[source]

Return two instances of the BlockMatrix class.

Parameters:

tpmats (sequence of TPMatrix’es or single BlockMatrix) – There can be both boundary matrices from inhomogeneous Dirichlet or Neumann conditions, as well as regular matrices.

Note

Use BlockMatrix directly if you do not have any inhomogeneous boundary conditions.

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

Note

The tensor product matrices may be either boundary matrices, regular matrices, or a mixture of both.

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 tensor products of these bases

\[\begin{split}TD &= K0 \otimes K1 \otimes SD \\ TT &= K0 \otimes K1 \otimes 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 assembled 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 = VectorSpace(TD)
>>> Q = CompositeSpace([VT, TD])
>>> up = TrialFunction(Q)
>>> vq = TestFunction(Q)
>>> u, p = up
>>> v, q = vq
>>> A00 = inner(grad(v), grad(u))
>>> A01 = inner(div(v), p)
>>> A10 = inner(q, div(u))
>>> M = BlockMatrix(A00+A01+A10)
assemble(format=None)[source]

Assemble matrices in scipy sparse format

Parameters:

format (str or None, optional) – The format of the sparse scipy matrix. Using config['matrix']['block']['assemble'] if None.

diags(it=None, format=None)[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 or None, optional) – where n is dimensions-1. These are the indices into the diagonal axes, or the axes with Fourier bases.

  • format (str or None, optional) – The format of the returned matrix. See Scipy sparse matrices If None, then use default for TPMatrix.

get_mats(return_first=False)[source]

Return flattened list of matrices in self

Parameters:

return_first (bool, optional) – Return just the first matrix in the loop if True

matvec(v, c, format=None, use_scipy=None)[source]

Compute matrix vector product

\[c = A v\]

where \(A\) is the self block matrix and \(v,c\) are flattened instances of the class Function.

Parameters:

  • v (Function)

  • c (Function)

  • format (str, optional) – The format of the matrices used for the matvec. See Scipy sparse matrices

  • use_scipy (boolean, optional) – Whether to assemble and use scipy’s bmat for the matvec, or to use the matvec methods of this BlockMatrix’s TPMatrices. Using config['matrix']['block']['use_scipy'] if use_scipy is None

Returns:

c

Return type:

Function

solve(b, u=None, constraints=())[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.

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

Bases: SparseMatrix

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, constraints=())[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

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

Note

Vectors may be one- or multidimensional.

class shenfun.matrixbase.ScipyMatrix(mats)[source]

Bases: csr_matrix

matvec(v, c, 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

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

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

Bases: MutableMapping

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

Note

The matrix format and storage is similar to Scipy’s dia_matrix. The format is chosen because spectral matrices often are computed by hand and presented in the literature as banded matrices. Note that a SparseMatrix can easily be transformed to any of Scipy’s formats using the diags method. However, Scipy’s matrices are not implemented to act along different axes of multidimensional arrays, which is required for tensor product matrices, see TPMatrix. Hence the need for this SparseMatrix class.

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}
>>> S = SparseMatrix(d, (N, N))
>>> dict(S)
{-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)}
>>> S = SparseMatrix(d, (N, N))
>>> dict(S)
{-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

copy()[source]

Return SparseMatrix deep copy of self

diags(format=None, scaled=True)[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

    • csc - Compressed sparse column

    Using config['matrix']['sparse']['diags'] setting if format is None

  • scaled (bool, optional) – Return matrix scaled by the constant self.scale if True

Note

This method returns the matrix scaled by self.scale if keyword scaled is True.

get_solver()[source]

Return appropriate solver for self

Note

Fall back on generic Solve, which is using Scipy sparse matrices with splu/spsolve. This is still pretty fast.

incorporate_scale()[source]

Modifies matrix such that self.scale = 1

matvec(v, c, format=None, 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

    • numba - Numba implementation that may be implemented in subclass

    Using config['matrix']['sparse']['matvec'] setting if format is None

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

solve(b, u=None, axis=0, constraints=())[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

  • constraints (tuple of 2-tuples) – The 2-tuples represent (row, val) The constraint indents the matrix row and sets b[row] = val

Note

Vectors may be one- or multidimensional.

class shenfun.matrixbase.SpectralMatDict[source]

Bases: dict

Dictionary for inner product matrices

Matrices are looked up with keys that are one of:

((test, k), (trial, l))
((test, k), (trial, l), measure)

where test and trial are classes subclassed from SpectralBase and k and l are integers >= 0 that determines how many times the test or trial functions should be differentiated. The measure is optional.

class shenfun.matrixbase.SpectralMatrix(test, trial, scale=1.0, measure=1, assemble=None, kind=None, fixed_resolution=None)[source]

Bases: SparseMatrix

Base class for inner product matrices.

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 column.

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

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

  • measure (number or Sympy expression, optional) – A function of the reference coordinate.

  • assemble (None or str, optional) – Determines how to perform the integration,

    • ‘quadrature’ (default)

    • ‘exact’

    • ‘adaptive’

    Exact and adaptive should result in the same matrix. Exact computes the integral using Sympy integrate, whereas adaptive makes use of adaptive quadrature through scipy.

  • kind (None or str, optional) – Alternative kinds of methods.

    • ‘implemented’ - Hardcoded implementations

    • ‘stencil’ - Use orthogonal bases and stencil-matrices

    • ‘vandermonde’ - Use Vandermonde matrix

    The default is to first try to look for implemented kind, and if that fails try first ‘stencil’ and then finally fall back on vandermonde. Vandermonde creates a dense matrix of size NxN, so it should be avoided (e.g., by implementing the matrix) for large N.

  • fixed_resolution (None or str, optional) – A fixed number of quadrature points used to compute the matrix. If ‘fixed_resolution’ is set, then assemble is set to ‘quadrature’ and kind is set to ‘vandermonde’.

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:

>>> from shenfun import FunctionSpace, SpectralMatrix
>>> SD = FunctionSpace(16, 'C', bc=(0, 0))
>>> M = SpectralMatrix((SD, 0), (SD, 0))

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

>>> A = SpectralMatrix((SD, 0), (SD, 2))

where (SD, 2) signals that we use the second derivative of this trial function. A more natural way to do the same thing is to

>>> from shenfun import TrialFunction, TestFunction, inner, Dx
>>> u = TrialFunction(SD)
>>> v = TestFunction(SD)
>>> A = inner(v, Dx(u, 0, 2))

where Dx() is a partial (or ordinary) derivative.

assemble(method)[source]

Return diagonals of SpectralMatrix

Parameters:

method (str) – Type of integration

  • ‘exact’

  • ‘quadrature’

Note

Subclass SpectralMatrix and overload this method in order to provide a fast and accurate implementation of the matrix representing an inner product. See the matrix modules in either one of

  • legendre.matrix

  • chebyshev.matrix

  • chebyshevu.matrix

  • ultraspherical.matrix

  • fourier.matrix

  • laguerre.matrix

  • hermite.matrix

  • jacobi.matrix

Example

The mass matrix for Chebyshev polynomials is

\[(T_j, T_i)_{\omega} = \frac{c_i \pi}{2}\delta_{ij},\]

where \(c_0=2\) and \(c_i=1\) for integer \(i>0\). We can implement this as

>>> from shenfun import SpectralMatrix
>>> class Bmat(SpectralMatrix):
...     def assemble(self, method):
...         test, trial = self.testfunction, self.trialfunction
...         ci = np.ones(test[0].N)
...         ci[0] = 2
...         if test[0].quad == 'GL' and method != 'exact':
...             # Gauss-Lobatto quadrature inexact at highest polynomial order
...             ci[-1] = 2
...         return {0: ci*np.pi/2}

Here {0: ci*np.pi/2} is the 0’th diagonal of the matrix. Note that test and trial are two-tuples of (instance of :class:.SpectralBase`, number)`, where the number represents the number of derivatives. For the mass matrix the number will be 0. Also note that the length of the diagonal must be correct.

property axis

Return the axis of the TensorProductSpace this matrix is created for

matvec(v, c, format=None, 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

    • numba - Numba implementation that may be implemented in subclass

    Using config['matrix']['sparse']['matvec'] setting if format is None

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

property tensorproductspace

Return the TensorProductSpace this matrix has been computed for

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

Bases: object

Tensor product matrix

A TensorProductSpace is the tensor product of D univariate function spaces. A normal matrix (a second order tensor) is assembled from bilinear forms (i.e., forms containing both test and trial functions) on one univariate function space. A bilinear form on a tensor product space will assemble to D outer products of such univariate matrices. That is, for a two-dimensional tensor product you get fourth order tensors (outer product of two matrices), and three-dimensional tensor product spaces leads to a sixth order tensor (outer product of three matrices). This class contains D second order matrices. The complete matrix is as such the outer product of these D matrices.

Note that the outer product of two matrices often is called the Kronecker product.

Parameters:

  • mats (sequence, or sequence of sequence of matrices) – Instances of SpectralMatrix or SparseMatrix The length of mats is the number of dimensions of the TensorProductSpace

  • testspace (Function space) – The test TensorProductSpace

  • trialspace (Function space) – The trial TensorProductSpace

  • scale (array, optional) – Scalar multiple of matrices. Must have ndim equal to the number of dimensions in the TensorProductSpace, and the shape must be 1 along any directions with a nondiagonal matrix.

  • global_index (2-tuple, optional) – Indices (test, trial) into mixed space CompositeSpace.

  • testbase (CompositeSpace, optional) – Instance of the base test space

  • trialbase (CompositeSpace, optional) – Instance of the base trial space

copy()[source]

Return TPMatrix deep copy of self

property dimensions

Return dimension of TPMatrix

shenfun.matrixbase.assemble_sympy(test, trial, measure=1, implicit=True, assemble='exact')[source]

Return sympy representation of mass matrix

Parameters:

  • test (TestFunction or 2-tuple of SpectralBase, int) – If 2-tuple, then the integer represents the number of derivatives, which should be zero for this function

  • trial (Like test but representing trial function)

  • measure (Number or Sympy function) – Function of coordinate

  • implicit (bool, optional) – Whether to use unevaluated Sympy functions instead of the actual values of the diagonals.

  • assemble (str, optional) –

    • ‘exact’

    • ‘quadrature’

Example

>>> from shenfun import assemble_sympy, TrialFunction, TestFunction
>>> N = 8
>>> D = FunctionSpace(N, 'C', bc=(0, 0))
>>> v = TestFunction(D)
>>> u = TrialFunction(D)
>>> assemble_sympy(v, u)
(KroneckerDelta(i, j) - KroneckerDelta(i, j + 2))*h(i) - (KroneckerDelta(j, i + 2) - KroneckerDelta(i + 2, j + 2))*h(i + 2)

Note that when implicit is True, then h(i) represents the l2-norm, or the L2-norm (if exact) of the orthogonal basis.

>>> D = FunctionSpace(N, 'C', bc={'left': {'N': 0}, 'right': {'N': 0}})
>>> u = TrialFunction(D)
>>> v = TestFunction(D)
>>> assemble_sympy(v, u, implicit=True)
(KroneckerDelta(i, j) + KroneckerDelta(i, j + 2)*s2(j))*h(i) + (KroneckerDelta(j, i + 2) + KroneckerDelta(i + 2, j + 2)*s2(j))*h(i + 2)*k2(i)
>>> assemble_sympy(v, u, implicit=False)
-i**2*(-j**2*KroneckerDelta(i + 2, j + 2)/(j + 2)**2 + KroneckerDelta(j, i + 2))*h(i + 2)/(i + 2)**2 + (-j**2*KroneckerDelta(i, j + 2)/(j + 2)**2 + KroneckerDelta(i, j))*h(i)

Here the implicit version uses ‘k’ for the diagonals of the test function, and ‘s’ for the trial function. The number represents the location of the diagonal, so ‘s2’ is the second upper diagonal of the stencil matrix of the trial function.

You can get the diagonals like this: >>> import sympy as sp >>> i, j = sp.symbols(‘i,j’, integer=True) >>> M = assemble_sympy(v, u, implicit=False) >>> M.subs(j, i) # main diagonal pi*i**4/(2*(i + 2)**4) + pi/2 >>> M.subs(j, i+2) # second upper diagonal -pi*i**2/(2*(i + 2)**2) >>> M.subs(j, i-2) # second lower diagonal -pi*(i - 2)**2/(2*i**2)

i is the row number, so the last one starts for i=2.

shenfun.matrixbase.check_sanity(A, test, trial, measure=1, assemble='quadrature', kind='vandermonde', fixed_resolution=None)[source]

Sanity check for matrix.

Test that created matrix agrees with quadrature computed using a memory-consuming Vandermonde implementation.

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) – Function in the physical coordinate. Gets mapped to reference domain.

  • assemble (str, optional) – Determines how to perform the integration we compare with

    • ‘exact’

    • ‘adaptive’

    • ‘quadrature’

  • kind (str, optional) – The type of quadrature to compare with.

    • ‘vandermonde’

    • ‘stencil’

  • fixed_resolution (None or str, optional) – A fixed number of quadrature points used to compute the matrix. If ‘fixed_resolution’ is set, then assemble is set to ‘quadrature’ and kind is set to ‘vandermonde’.

shenfun.matrixbase.extract_bc_matrices(mats)[source]

Extract boundary matrices from list of mats

Parameters:

mats (list of list of instances of TPMatrix or) – SparseMatrix

Returns:

list of boundary matrices.

Return type:

list

Note

The mats list is modified in place since boundary matrices are extracted.

shenfun.matrixbase.extract_diagonal_matrix(M, lowerband=None, upperband=None, abstol=1e-10, reltol=1e-10)[source]

Return SparseMatrix version of dense matrix M

Parameters:

  • M (Numpy array of ndim=2 or sparse scipy matrix)

  • lowerband (int or None) – Assumed lower bandwidth of M

  • upperband (int or None) – Assumed upper bandwidth of M

  • 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_simplified_tpmatrices(tpmats)[source]

Return copy of tpmats list, where diagonal matrices have been simplified and placed in scale arrays.

Parameters:

tpmats – Instances of TPMatrix

Returns:

List of TPMatrix’es, that have been simplified

Return type:

List[TPMatrix]

shenfun.spectralbase module

This module contains classes for working with global spectral Galerkin methods in one dimension. All spectral functions have the following expansions on the real line

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

where \(\mathcal{I}\) is a index set of the basis, which differs from space to space. \(\phi_k\) is the k’t basis function and the function space is the span of these basis functions.

See the [documentation](https://shenfun.readthedocs.io) for more details.

class shenfun.spectralbase.BoundaryConditions(bc, domain=None)[source]

Bases: dict

Boundary conditions for SpectralBase.

Parameters:

  • bc (str, dict or n-tuple) – The dictionary must have keys ‘left’ and ‘right’, to describe boundary conditions on the left and right boundaries, and another dictionary on left or right to describe the condition. For example, specify Dirichlet on both ends with:

    {'left': {'D': a}, 'right': {'D': b}}

    for some values a and b. Specify Neumann as:

    {'left': {'N': a}, 'right': {'N': b}}

    A mixture of 3 conditions:

    {'left': {'N': a}, 'right': {'D': b, 'N': c}}

    Etc. Any combination should be possible, and it should also be possible to use higher order derivatives with N2, N3 etc.

    If bc is an n-tuple, then we assume the basis function is:

    (None, a) - {'right': {'D': a}}
(a, None) - {'left': {'D': a}}
(a, b) - {'left': {'D': a}, 'right': {'D': b}}
(a, b, c, d) - {'left': {'D': a, 'N': b}, 'right': {'D': c, 'N': d}}
(a, b, c, d, e, f) - {'left': {'D': a, 'N': b, 'N2': c},
                      'right': {'D': d, 'N': e, 'N2': f}}
etc.

    If bc is a single string, then we assume the boundary conditions are described directly for generic function u, like:

    'u(-1)=0&&u(1)=0' - Dirichlet on boundaries x=-1 and x=1
'u'(-1)=0&&u'(1)=0' - Neumann on boundaries x=-1 and x=1
'u(-2)=a&&u(2)=b' - Dirichlet with values a and b on boundaries
    x=-2 and x=2.
etc.

  • domain (2-tuple of numbers, optional) – The domain, if different than (-1, 1).

class shenfun.spectralbase.MixedFunctionSpace(bases)[source]

Bases: object

Class for composite bases in 1D

Parameters:

bases (list) – List of instances of SpectralBase

dim()[source]

Return dimension of self (degrees of freedom)

dims()[source]

Return dimensions (degrees of freedom) for MixedFunctionSpace

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

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

Parameters:

  • padding_factor (number or tuple, optional) – Scale the number of quadrature points in self with this number, or tuple of numbers, one for each dimension.

  • dealias_direct (bool, optional) – Used only if padding_factor=1. Sets the 1/3 highest frequencies to zeros.

num_components()[source]

Return number of bases in mixed basis

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

L2_norm_sq(i)[source]

Return square of L2-norm

\[\| \phi_i \|^2_{\omega} = (\phi_i, \phi_i)_{\omega} = \int_{I} \phi_i \overline{\phi}_i \omega dx\]

where \(\phi_i\) is the i’th orthogonal basis function for the orthogonal basis in the given family.

Parameters:

i (int) – The number of the orthogonal basis function

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.

Return type:

array

backward(input_array=None, output_array=None, kind=None, mesh=None)[source]

Compute backward (inverse) transform

Parameters:

  • input_array (array, optional) – Expansion coefficients

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

  • kind (str, optional) –

    • ‘fast’ - Use fast transform on regular quadrature points

    • ‘recursive’ - Use low-memory implementation (only for polynomials)

    • ‘vandermonde’ - use Vandermonde on regular quadrature points

  • mesh (str or functionspace, optional) –

    • ‘quadrature’ - use quadrature mesh of self

    • ‘uniform’ - use uniform mesh

    • A function space with the same mesh distribution as self

Note

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

basis_function(i=0, x=x)[source]

Return basis function i

Parameters:

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

  • x (sympy Symbol, optional)

broadcast_to_ndims(x)[source]

Return 1D array x as an array of shape according to the dimensions() of the TensorProductSpace class that the instance of this class 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 FunctionSpace, TensorProductSpace, comm
>>> K0 = FunctionSpace(8, 'F', dtype='D')
>>> K1 = FunctionSpace(8, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (K0, K1), modify_spaces_inplace=True)
>>> x = np.arange(4)
>>> y = K0.broadcast_to_ndims(x)
>>> print(y.shape)
(4, 1)
cartesian_mesh(kind='quadrature')[source]

Return Cartesian mesh

Parameters:

kind (str, optional) –

  • ‘quadrature’ - Use quadrature mesh

  • ‘uniform’ - Use uniform mesh

Note

For Cartesian problems this function returns the same mesh as 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 space 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

Note

The domain factor is the length of the reference domain over the length of the true domain.

property dtype

Return datatype function space 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 function i at points x

Parameters:

  • x (float or array of floats)

  • i (int, optional) – Basis function 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 function 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 function 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

forward(input_array=None, output_array=None, kind=None)[source]

Compute forward transform

Parameters:

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

  • output_array (array, optional) – Expansion coefficients

  • kind (str, optional) –

    • ‘fast’ - use fast transform if implemented

    • ‘vandermonde’ - Use Vandermonde matrix

    • ‘recursive’ - Use low-memory implementation (only for polynomials)

Note

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

get_adaptive(fun=None, reltol=1e-12, abstol=1e-15)[source]

Return space (otherwise as self) with number of quadrature points determined by fitting fun

Returns:

A new space with adaptively found number of quadrature points

Return type:

SpectralBase

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

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

Parameters:

  • padding_factor (float, optional) – Create output array of shape padding_factor times non-padded shape

  • 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.

  • kwargs (keyword arguments) – Any other keyword arguments used in the creation of the bases.

Returns:

The space to be used for dealiasing

Return type:

SpectralBase

get_homogeneous(**kwargs)[source]

Return space (otherwise as self) with homogeneous boundary conditions

Parameters:

kwargs (keyword arguments) – Any keyword arguments used in the creation of the bases.

Returns:

A new space with homogeneous boundary conditions, otherwise as self.

Return type:

SpectralBase

get_measured_array(array)[source]

Return array times Jacobian determinant

Parameters:

array (array)

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, map_true_domain=False)[source]

Return weights times measure

Parameters:

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

  • measure (1 or sympy.Expr)

get_orthogonal(**kwargs)[source]

Return orthogonal space (otherwise as self)

Returns:

The orthogonal space in the same family, and otherwise as self.

Return type:

SpectralBase

get_refined(N, **kwargs)[source]

Return space (otherwise as self) with N quadrature points

Parameters:

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

  • kwargs (keyword arguments) – Any other keyword arguments used in the creation of the bases.

Returns:

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

Return type:

SpectralBase

get_testspace(kind='Galerkin', **kwargs)[source]

Return appropriate test space

Parameters:

  • kind (str, optional) –

    • ‘Galerkin’ - or ‘G’

    • ‘Petrov-Galerkin’ - or ‘PG’

    If kind = ‘Galerkin’ or ‘G’, then return self, corresponding to a regular Galerkin method.

    If kind = ‘Petrov-Galerkin’ or ‘PG’, then return a test space that makes use of the testbases Phi1, Phi2, Phi3 or Phi4, where the actual choice is made based on the number of boundary conditions in self. One boundary condition uses Phi1, two uses Phi2 etc.

    For Petrov-Galerkin the returned basis corresponds to \(\{\phi^{(k)}_n\}\), where

    \[\phi_n^{(k)} = \frac{(1-x^2)^k}{h^{(k)}_{n+k}}\frac{d^k}{dx^k}Q_{n+k}\]

  • kwargs (keyword arguments, optional) – Any other keyword arguments used in the creation of the test bases. Only used if kind=’Petrov-Galerkin’.

Note

If self is not a Jacobi polynomial basis, then return self, corresponding to a regular Galerkin method.

Return type:

Instance of SpectralBase

get_unplanned(**kwargs)[source]

Return unplanned space (otherwise as self)

Parameters:

kwargs (keyword arguments, optional) – Any keyword argument used in the creation of the unplanned space. Could be any one of

  • quad

  • domain

  • dtype

  • padding_factor

  • dealias_direct

  • coordinates

  • bcs

  • scaled

Not all will be applicable for all spaces.

Returns:

Space not planned for a TensorProductSpace

Return type:

SpectralBase

l2_norm_sq(i=None)[source]

Return square of l2-norm

\[\| u \|^2_{N,\omega} = (u, u)_{N,\omega} = \sun_{j=0}^{N-1} u(x_j) \overline{u}(x_j) \omega_j\]

where \(u=\{\phi_i\}_{i=0}^{N-1}\) and \(\phi_i\) is the i’th orthogonal basis function in the given family.

Parameters:

i (None or int) – If None then return the square of the l2-norm for all i=0, 1, …, N-1. Else, return for given i.

map_expression_true_domain(f, x=None)[source]

Return expression f mapped to true domain as a function of the reference coordinate

Parameters:

  • f (Sympy expression)

  • x (Sympy symbol, optional) – coordinate

map_reference_domain(x)[source]

Return true point x mapped to reference domain

Parameters:

x (coordinate or array of points)

map_true_domain(x)[source]

Return reference point x mapped to true domain

Parameters:

x (coordinate or array of points)

mesh(bcast=True, map_true_domain=True, kind='quadrature')[source]

Return the computational mesh

Parameters:

  • bcast (bool) – Whether or not to broadcast to dimensions() if an instance of this basis belongs to a TensorProductSpace. The returned mesh then has shape one in all ndims-1 dimensions apart from self.axis.

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

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

Note

For curvilinear problems this function returns the computational mesh. For the corresponding Cartesian mesh, use cartesian_mesh()

orthogonal_basis_function(i=0, x=x)[source]

Return the orthogonal basis function i

Parameters:

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

  • x (sympy Symbol, optional)

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.

Note

The weight of the space is included in the returned quadrature weights.

property rank

Return rank of function space

Note

This is 1 for MixedFunctionSpace and 0 otherwise

scalar_product(input_array=None, output_array=None, kind=None)[source]

Compute weighted scalar product

Parameters:

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

  • output_array (array, optional) – The result of the scalar product

  • kind (str, optional) –

    • ‘fast’ - use fast transform if implemented

    • ‘vandermonde’ - use Vandermonde matrix

    • ‘recursive’ - Use low-memory implementation (only for polynomials)

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 space

sympy_L2_norm_sq(i=i)[source]

Return sympy function for square of L2-norm

\[\| \phi_i \|^2_{N,\omega} = (\phi_i, \phi_i)_{\omega} = \int_{I} \phi_i \overline{\phi}_i \omega dx\]

where \(\phi_i\) is the i’th orthogonal basis function in the given family.

Parameters:

i (Sympy Symbol, optional)

Return type:

Sympy Function

sympy_l2_norm_sq(i=i)[source]

Return sympy function for square of l2-norm

\[\| \phi_i \|^2_{N,\omega} = (\phi_i, \phi_i)_{N,\omega} = \sun_{j=0}^{N-1} \phi_i(x_j) \overline{\phi}_i(x_j) \omega_j\]

where \(\phi_i\) is the i’th orthogonal basis function in the given family.

Parameters:

  • i (Sympy Symbol, optional)

  • implicit (bool, optional) – Whether to use an unevaluated Sympy function instead of the actual value of the l2-norm.

Return type:

Sympy Function

property tensor_rank

Return tensor rank of function space

Note

This is the number of free indices in the tensor, 0 for scalar, 1 for vector etc. It is None for a composite space that is not a tensor.

property tensorproductspace

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

Note

A 1D function space may be part of several :class:`.TensorProductSpace`s, 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 (orthogonal) basis of the family.

Evaluates basis function \(\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 orthogonal basis functions for either family. The true Vandermonde matrix of a basis is obtained through SpectralBase.evaluate_basis_all().

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

Return the wavenumbermesh

Parameters:

bcast (bool) – Whether or not to broadcast to dimensions() if an instance of this basis belongs to a TensorProductSpace

weight(x=None)[source]

Weight of inner product space

Parameters:

x (coordinate)

shenfun.spectralbase.getBCGeneric(CompositeBase)[source]

Dynamic class factory for boundary spaces

Parameters:

CompositeBase (SpectralBase) – Dynamic inheritance, using base class as either one of

  • chebyshev.bases.CompositeBase

  • chebyshevu.bases.CompositeBase

  • legendre.bases.CompositeBase

  • ultraspherical.bases.CompositeBase

  • jacobi.bases.CompositeBase

  • laguerre.bases.CompositeBase

Returns:

Return type:

Either one of the classes

shenfun.spectralbase.getCompositeBase(Orthogonal)[source]

Dynamic class factory for Composite base class

Parameters:

Orthogonal (SpectralBase) – Dynamic inheritance, using base class as either one of

Returns:

  • chebyshev.bases.CompositeBase

  • chebyshevu.bases.CompositeBase

  • legendre.bases.CompositeBase

  • ultraspherical.bases.CompositeBase

  • jacobi.bases.CompositeBase

  • laguerre.bases.CompositeBase

Return type:

Either one of the classes

shenfun.spectralbase.get_norm_sq(v, u, method)[source]

Return square of weighted norm

\[`(\phi_i, \phi_i)_w`\]

where \(\phi_i\) is the orthogonal basis for a spectral family.

Parameters:

  • v (instance of SpectralBase) – Orthogonal test function

  • N (instance of SpectralBase) – Orthogonal trial function

  • method (str) – Type of integration

    • ‘exact’

    • ‘quadrature’

shenfun.spectralbase.inner_product(test, trial, measure=1, assemble=None, kind=None, fixed_resolution=None)[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 function is differentiated. The test represents the matrix row

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

  • measure (function of coordinate, optional) – The measure is in physical coordinates, not in the reference domain.

  • assemble (None or str, optional) – Determines how to perform the integration

    • ‘quadrature’ (default)

    • ‘exact’

    • ‘adaptive’

  • kind (None or str, optional) – The kind of method used to do quadrature.

    • ‘implemented’

    • ‘stencil’

    • ‘vandermonde’

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 ShenDirichlet
>>> SD = ShenDirichlet(6)
>>> B = inner_product((SD, 0), (SD, 0))
>>> d = {-2: -np.pi/2,
...       0: np.array([1.5*np.pi, np.pi, np.pi, np.pi]),
...       2: -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.CompositeSpace(spaces)[source]

Bases: object

Class for composite tensorproductspaces.

The mixed spaces are Cartesian products of TensorProductSpace or other instances of CompositeSpace.

Parameters:

spaces (list) – List of spaces

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)

dims()[source]

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

dims_composite()[source]

Return dimensions per axis (degrees of freedom) of all spaces in self

eval(points, coefficients, output_array=None, method=1)[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. Method 0 is a low-memory cython version. Using method = 1 (default) leads to a faster cython implementation that, on the downside, uses more memory. The final, method = 2, is a python implementation used mainly for verification.

get_offsets()[source]

Return offsets for spaces

global_shape(forward_output=False)[source]

Return global shape for CompositeSpace

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.

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)

property rank

Return rank of space

Note

This is 1 for all composite spaces, since we use a flat storage of all composite arrays.

shape(forward_output=False)[source]

Return shape of arrays for CompositeSpace

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 CompositeSpace 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 CompositeSpace

slice()[source]

The slices of dofs for each dimension

property tensor_rank

Return rank of tensor

Note

This is the number of free indices in the tensor, 0 for scalar, 1 for vector etc. It is None for a composite space that is not a tensor.

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.TensorProductSpace(comm, bases, axes=None, dtype=None, slab=False, collapse_fourier=False, backward_from_pencil=False, coordinates=None, modify_spaces_inplace=False, **kw)[source]

Bases: PFFT

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 (SpectralBase)

  • 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. See Coordinates.

  • modify_spaces_inplace (bool, optional) – Whether or not a copy should be made of the input functionspaces. If True, then the input spaces will be modified inplace.

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

cartesian_mesh(map_true_domain=True, kind='quadrature')[source]

Return Cartesian mesh

Parameters:

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

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

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

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

property domains

Return domains for all directions

dtype(forward_output=False)[source]

The type of transformed arrays

Parameters:

forward_output (bool, optional) – If True then return dtype of an array that is the result of a forward transform. Otherwise, return the dtype of an array that is input to a forward transform.

eval(points, coefficients, output_array=None, method=1)[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 method 0 is a low-memory cython version. Using method = 1 (default) leads to a faster cython implementation that, on the downside, uses more memory. The final, method = 2, is a python implementation.

get_adaptive(fun=None, reltol=1e-12, abstol=1e-15)[source]

Return space (otherwise as self) with number of quadrature points determined by fitting fun

Returns:

A new space with adaptively found number of quadrature points

Return type:

TensorProductSpace

Note

Only spaces defined with N=0 quadrature points are adapted, the others are kept as is

Example

>>> from shenfun import FunctionSpace, Function, TensorProductSpace, comm
>>> import sympy as sp
>>> r, theta = psi = sp.symbols('x,y', real=True, positive=True)
>>> rv = (r*sp.cos(theta), r*sp.sin(theta))
>>> B0 = FunctionSpace(0, 'C', domain=(0, 1))
>>> F0 = FunctionSpace(0, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (B0, F0), coordinates=(psi, rv))
>>> u = Function(T, buffer=((1-r)*r)**4*(sp.cos(10*theta)))
>>> print(u.global_shape)
(9, 11)
get_dealiased(padding_factor=1.5, dealias_direct=False)[source]

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

Parameters:

  • padding_factor (number or tuple, optional) – Scale the number of quadrature points in self with this number, or tuple of numbers, one for each dimension.

  • dealias_direct (bool, optional) – Used only if padding_factor=1. Sets the 1/3 highest frequencies to zeros.

get_homogeneous(**kwargs)[source]

Return space with homogeneous boundary conditions, and otherwise as self.

Parameters:

kwargs (keyword arguments) – Any arguments used in the creation of the different bases.

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_orthogonal()[source]

Return tensor product space using the orthogonal basis for all directions

get_refined(N)[source]

Return space (otherwise as self) refined to new shape

Parameters:

N (Number or sequence of integers) – The global shape of the new space. If N is a number, then refine the global shape by multiplying self’s global shape by N. If N is an array, then N is treated as the new global shape of the returned space.

get_testspace(kind='Galerkin', **kwargs)[source]

Return appropriate test space

Parameters:

  • kind (str, optional) –

    • ‘Galerkin’ - or ‘G’

    • ‘Petrov-Galerkin’ - or ‘PG’

    If kind = ‘Galerkin’, then return self, corresponding to a regular Galerkin method.

    If kind = ‘Petrov-Galerkin’ or ‘PG’, then return a test space that makes use of the testbases Phi1, Phi2, Phi3 or Phi4, where the choice is made based on the number of boundary conditions in self. One boundary condition uses Phi1, two uses Phi2 etc.

  • kwargs (keyword arguments, optional) – Any other keyword arguments used in the creation of the test bases.

Note

Petrov-Galerkin is only possible for Jacobi polynomials. If mixed with, e.g., Fourier spaces, then regular Galerkin is used for the non-polynomial bases.

Return type:

Instance of TensorProductSpace

get_unplanned(tensorproductspace=False, **kwargs)[source]

Return unplanned bases otherwise as self. Or return a new TensorProductSpace from the unplanned bases.

Parameters:

  • tensorproductspace (boolean) – if True, then return a new TensorProductSpace if False, return only the bases

  • kwargs (keyword arguments) – Any arguments used in the creation of the different bases.

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.

local_cartesian_mesh(map_true_domain=True, kind='quadrature')[source]

Return local Cartesian mesh

Parameters:

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

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

local_mesh(bcast=False, map_true_domain=True, kind='quadrature')[source]

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

Parameters:

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

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

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

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

Return list of local wavenumbers

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(map_true_domain=True, kind='quadrature')[source]

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

Parameters:

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

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

num_components()[source]

Return number of scalar spaces in TensorProductSpace

property rank

Return tensor rank of TensorProductSpace

shape(forward_output=True)[source]

The local (to each processor) shape of data

Parameters:

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

size(forward_output=False)[source]

Return number of elements in TensorProductSpace

slice()[source]

The slices of dofs for each dimension

property use_fixed_gauge

Return True if TensorProductSpace contains only spaces with Neumann boundary conditions.

property volume

Return computational volume

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

Return list of wavenumbers

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.TensorSpace(space)[source]

Bases: VectorSpace

A special CompositeSpace for second rank tensors.

A TensorProductSpace created by a tensorproduct of 2 1D bases, will have tensors of length 4. A TensorProductSpace created from 3 1D bases will have tensors of length 9.

Parameters:

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

num_components()[source]

Return number of spaces in mixed space

property tensor_rank

Return rank of tensor

Note

This is the number of free indices in the tensor, 0 for scalar, 1 for vector etc. It is None for a composite space that is not a tensor.

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

Bases: CompositeSpace

Vector of TensorProductSpace instances.

This is simply a special CompositeSpace where the number of instances of the TensorProductSpace class equals the geometrical dimension of the problem.

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

Parameters:

  • space (TensorProductSpace or list of ndim instances of the)

  • :class:`.TensorProductSpace` class to create vector from.

num_components()[source]

Return number of spaces in mixed space

shape(forward_output=False)[source]

Return shape of arrays for VectorSpace

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 tensor_rank

Return rank of tensor

Note

This is the number of free indices in the tensor, 0 for scalar, 1 for vector etc. It is None for a composite space that is not a tensor.

shenfun.coordinates module

class shenfun.coordinates.Coordinates(psi, rv, assumptions=True, replace=(), measure=<function count_ops>)[source]

Bases: object

Class for handling curvilinear coordinates

Parameters:

  • psi (tuple or sp.Symbol) – The new coordinates

  • rv (tuple) – The position vector in terms of the new coordinates

  • assumptions (Sympy assumptions) – One or more Sympy assumptions

  • replace (sequence of two-tuples) – Use Sympy’s replace with these two-tuples

  • measure (Python function to replace Sympy’s count_ops.) –

    For example, to discourage the use of powers in an expression use:

    def discourage_powers(expr):
    POW = sp.Symbol('POW')
    count = sp.count_ops(expr, visual=True)
    count = count.replace(POW, 100)
    count = count.replace(sp.Symbol, type(sp.S.One))
    return count
get_cartesian_basis()[source]

Return Cartesian basis vectors

get_christoffel_second()[source]

Return Christoffel symbol of second kind

get_contravariant_basis()[source]

Return contravariant basisvectors

get_contravariant_metric_tensor()[source]

Return contravariant metric tensor

get_covariant_basis()[source]

Return covariant basisvectors

get_covariant_metric_tensor()[source]

Return covariant metric tensor

get_det_g(covariant=True)[source]

Return determinant of covariant metric tensor

get_normal_metric_tensor()[source]

Return normal metric tensor

get_scaling_factors()[source]

Return scaling factors

get_sqrt_det_g(covariant=True)[source]

Return square root of determinant of covariant metric tensor

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 or Numba.

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