shenfun.forms package

Submodules

shenfun.forms.arguments module

class shenfun.forms.arguments.Expr(basis, terms=None, scales=None, indices=None)[source]

Bases: object

Class for spectral Galerkin forms

An Expression that is linear in TestFunction, TrialFunction or Function. Exprs are used as input to inner() or project().

Parameters
  • basis (BasisFunction) – TestFunction, TrialFunction or Function

  • terms (Numpy array of ndim = 3) – Describes operations performed in Expr

    • Index 0: Vector component. If Expr is rank = 0, then terms[0] = 1. For vectors it equals ndim

    • Index 1: One for each term in the form. For example div(grad(u)) has three terms in 3D:

    \[\partial^2u/\partial x^2 + \partial^2u/\partial y^2 + \partial^2u/\partial z^2\]
    • Index 2: The operations stored as an array of length = dim

    The Expr div(grad(u)), where u is a scalar, is as such represented as an array of shape (1, 3, 3), 1 meaning it’s a scalar, the first 3 because the Expr consists of the sum of three terms, and the last 3 because it is 3D. The entire representation is:

    array([[[2, 0, 0],
            [0, 2, 0],
            [0, 0, 2]]])
    

    where the first [2, 0, 0] term has two derivatives in first direction and none in the others, the second [0, 2, 0] has two derivatives in second direction, etc.

  • scales (Numpy array of shape == terms.shape[:2]) – Representing a scalar multiply of each inner product

  • indices (Numpy array of shape == terms.shape[:2]) – Index into MixedTensorProductSpace. Only used when basis of form has rank > 0

Examples

>>> from shenfun import *
>>> from mpi4py import MPI
>>> comm = MPI.COMM_WORLD
>>> C0 = Basis(16, 'F', dtype='D')
>>> C1 = Basis(16, 'F', dtype='D')
>>> R0 = Basis(16, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (C0, C1, R0))
>>> v = TestFunction(T)
>>> e = div(grad(v))
>>> e.terms()
array([[[2, 0, 0],
        [0, 2, 0],
        [0, 0, 2]]])
>>> e2 = grad(v)
>>> e2.terms()
array([[[1, 0, 0]],

       [[0, 1, 0]],

       [[0, 0, 1]]])

Note that e2 in the example has shape (3, 1, 3). The first 3 because it is a vector, the 1 because each vector item contains one term, and the final 3 since it is a 3-dimensional tensor product space.

property argument

Return argument of Expr’s basis

property base

Return base BasisFunction used in Expr

basis()[source]

Return basis of Expr

property dimensions

Return ndim of Expr

eval(x, output_array=None)[source]

Return expression evaluated on x

Parameters

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

expr_rank()[source]

Return rank of Expr

function_space()[source]

Return function space of basis in Expr

index()[source]
indices()[source]

Return indices of Expr

num_components()[source]

Return number of components in Expr

num_terms()[source]

Return number of terms in Expr

property rank

Return rank of Expr’s basis

scales()[source]

Return scales of Expr

terms()[source]

Return terms of Expr

class shenfun.forms.arguments.BasisFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Bases: object

Base class for arguments to shenfun’s Exprs

Parameters
property argument

Return argument of basis

property base

Return base

property basespace

Return base space

property dimensions

Return dimensions of function space

expr_rank()[source]

Return rank of expression involving basis

function_space()[source]

Return function space of BasisFunction

index()[source]

Return index into base space

num_components()[source]

Return number of components in basis

offset()[source]

Return offset of this basis

The offset is the number of scalar TensorProductSpace`es ahead of this space in a :class:.MixedTensorProductSpace`.

property rank

Return rank of basis

class shenfun.forms.arguments.TestFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Bases: shenfun.forms.arguments.BasisFunction

Test function - BasisFunction with argument = 0

Parameters
  • space (TensorProductSpace)

  • index (int, optional) – Component of basis with rank > 0

  • basespace (The base MixedTensorProductSpace if space is a) – subspace.

  • offset (int) – The number of scalar spaces (i.e., :class:`.TensorProductSpace`es) ahead of this space

  • base (The base TestFunction)

property argument

Return argument of basis

class shenfun.forms.arguments.TrialFunction(space, index=0, basespace=None, offset=0, base=None)[source]

Bases: shenfun.forms.arguments.BasisFunction

Trial function - BasisFunction with argument = 1

Parameters
  • space (TensorProductSpace)

  • index (int, optional) – Component of basis with rank > 0

  • basespace (The base MixedTensorProductSpace if space is a) – subspace.

  • offset (int) – The number of scalar spaces (i.e., :class:`.TensorProductSpace`es) ahead of this space

  • base (The base TrialFunction)

property argument

Return argument of basis

class shenfun.forms.arguments.Function(space, val=0, buffer=None)[source]

Bases: shenfun.forms.arguments.ShenfunBaseArray, shenfun.forms.arguments.BasisFunction

Spectral Galerkin function for given TensorProductSpace or Basis()

The Function is the product of all 1D basis expansions, that for each dimension is defined like

\[u(x) = \sum_{k \in \mathcal{K}} \hat{u}_k \psi_k(x),\]

where \(\psi_k(x)\) are the trial functions and \(\{\hat{u}_k\}_{k\in\mathcal{K}}\) are the expansion coefficients. Here an index set \(\mathcal{K}=0, 1, \ldots, N\) is used to simplify notation.

For an M+1-dimensional TensorProductSpace with coordinates \(x_0, x_1, \ldots, x_M\) we get

\[u(x_{0}, x_{1}, \ldots, x_{M}) = \sum_{k_0 \in \mathcal{K}_0}\sum_{k_1 \in \mathcal{K}_1} \ldots \sum_{k_M \in \mathcal{K}_M} \hat{u}_{k_0, k_1, \ldots k_M} \psi_{k_0}(x_0) \psi_{k_1}(x_1) \ldots \psi_{k_M}(x_M),\]

where \(\mathcal{K}_j\) is the index set for the wavenumber mesh along axis \(j\).

Note that for a Cartesian mesh in 3D it would be natural to use coordinates \((x, y, z) = (x_0, x_1, x_2)\) and the expansion would be the simpler and somewhat more intuitive

\[u(x, y, z) = \sum_{l \in \mathcal{K}_0}\sum_{m \in \mathcal{K}_1} \sum_{n \in \mathcal{K}_2} \hat{u}_{l, m, n} \psi_{l}(x) \psi_{m}(y) \psi_{n}(z).\]

The Function’s values (the Numpy array) represent the \(\hat{u}\) array. The trial functions \(\psi\) may differ in the different directions. They are chosen when creating the TensorProductSpace.

Parameters
  • space (TensorProductSpace)

  • val (int or float) – Value used to initialize array

  • buffer (Numpy array, Function or sympy Expr) – If array it must be of correct shape. A sympy expression is evaluated on the quadrature mesh and forward transformed to create the buffer array.

Note

For more information, see numpy.ndarray

Examples

>>> from mpi4py import MPI
>>> from shenfun import Basis, TensorProductSpace, Function
>>> K0 = Basis(8, 'F', dtype='D')
>>> K1 = Basis(8, 'F', dtype='d')
>>> T = TensorProductSpace(MPI.COMM_WORLD, [K0, K1])
>>> u = Function(T)
>>> K2 = Basis(8, 'C', bc=(0, 0))
>>> T2 = TensorProductSpace(MPI.COMM_WORLD, [K0, K1, K2])
>>> v = Function(T2)
assign(u_hat)[source]

Assign self to u_hat of possibly different size

Parameters

u_hat (Function) – Function of possibly different shape than self. Must have the same function_space

backward(output_array=None, uniform=False)[source]

Return Function evaluated on quadrature mesh

eval(x, output_array=None)[source]

Evaluate Function at points

Parameters
  • points (float or array of floats)

  • coefficients (array) – Expansion coefficients

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

Examples

>>> import sympy as sp
>>> K0 = Basis(9, 'F', dtype='D')
>>> K1 = Basis(8, 'F', dtype='d')
>>> T = TensorProductSpace(MPI.COMM_WORLD, [K0, K1], axes=(0, 1))
>>> X = T.local_mesh()
>>> x, y = sp.symbols("x,y")
>>> ue = sp.sin(2*x) + sp.cos(3*y)
>>> ul = sp.lambdify((x, y), ue, 'numpy')
>>> ua = Array(T, buffer=ul(*X))
>>> points = np.random.random((2, 4))
>>> u = ua.forward()
>>> u0 = u.eval(points).real
>>> assert np.allclose(u0, ul(*points))
property forward_output

Return whether self is the result of a forward transform

mask_nyquist(mask=None)[source]

Set self to have zeros in Nyquist coefficients

refine(N, output_array=None)[source]

Return self with new number of quadrature points

Parameters

N (number or sequence of numbers) – The new number of quadrature points

Note

If N is smaller than for self, then a truncated array is returned. If N is greater than before, then the returned array is padded with zeros.

to_ortho(output_array=None)[source]

Project Function to orthogonal basis

class shenfun.forms.arguments.Array[source]

Bases: shenfun.forms.arguments.ShenfunBaseArray

Numpy array for TensorProductSpace

The Array is the result of a Function evaluated on its quadrature mesh.

The Function is the product of all 1D basis expansions, that for each dimension is defined like

\[u(x) = \sum_{k \in \mathcal{K}} \hat{u}_k \psi_k(x),\]

where \(\psi_k(x)\) are the trial functions and \(\{\hat{u}_k\}_{k\in\mathcal{K}}\) are the expansion coefficients. Here an index set \(\mathcal{K}=0, 1, \ldots, N\) is used to simplify notation.

For an M+1-dimensional TensorProductSpace with coordinates \(x_0, x_1, \ldots, x_M\) we get

\[u(x_{0}, x_{1}, \ldots, x_{M}) = \sum_{k_0 \in \mathcal{K}_0}\sum_{k_1 \in \mathcal{K}_1} \ldots \sum_{k_M \in \mathcal{K}_M} \hat{u}_{k_0, k_1, \ldots k_M} \psi_{k_0}(x_0) \psi_{k_1}(x_1) \ldots \psi_{k_M}(x_M),\]

where \(\mathcal{K}_j\) is the index set for the wavenumber mesh along axis \(j\).

Note that for a Cartesian mesh in 3D it would be natural to use coordinates \((x, y, z) = (x_0, x_1, x_2)\) and the expansion would be the simpler and somewhat more intuitive

\[u(x, y, z) = \sum_{l \in \mathcal{K}_0}\sum_{m \in \mathcal{K}_1} \sum_{n \in \mathcal{K}_2} \hat{u}_{l, m, n} \psi_{l}(x) \psi_{m}(y) \psi_{n}(z).\]

The Array’s values (the Numpy array) represent the left hand side, evaluated on the Cartesian quadrature mesh. With this we mean the \(u(x_i, y_j, z_k)\) array, where \(\{x_i\}_{i=0}^{N_0}\), \(\{y_j\}_{j=0}^{N_1}\) and \(\{z_k\}_{k=0}^{N_2}\) represent the mesh along the three directions. The quadrature mesh is then

\[(x_i, y_j, z_k) \quad \forall \, (i, j, k) \in [0, 1, \ldots, N_0] \times [0, 1, \ldots, N_1] \times [0, 1, \ldots, N_2]\]

The entire spectral Galerkin function can be obtained using the Function class.

Parameters
  • space (TensorProductSpace or SpectralBase)

  • val (int or float) – Value used to initialize array

  • buffer (Numpy array, Function or sympy Expr) – If array it must be of correct shape. A sympy expression is evaluated on the quadrature mesh and the result is used as buffer.

  • .. note:: For more information, see `numpy.ndarray <https (//docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html>`_)

Examples

>>> from mpi4py import MPI
>>> from shenfun import Basis, TensorProductSpace, Function
>>> K0 = Basis(8, 'F', dtype='D')
>>> K1 = Basis(8, 'F', dtype='d')
>>> FFT = TensorProductSpace(MPI.COMM_WORLD, [K0, K1])
>>> u = Array(FFT)
forward(output_array=None)[source]

Return Function used to evaluate Array

property forward_output

Return whether self is the result of a forward transform

offset()[source]

Return offset of this basis

The offset is the number of scalar TensorProductSpace`es ahead of this Arrays space in a :class:.MixedTensorProductSpace`.

shenfun.forms.arguments.Basis(N, family='Fourier', bc=None, dtype='d', quad=None, domain=None, scaled=None, padding_factor=1.0, dealias_direct=False, **kw)[source]

Return basis for one dimension

Parameters
  • N (int) – Number of quadrature points

  • family (str, optional) – Choose one of

    • Chebyshev or C,

    • Legendre or L,

    • Fourier or F,

    • Laguerre or La,

    • Hermite or H

  • bc (str or two-tuple, optional) – Choose one of

    • two-tuple (a, b) - Dirichlet boundary condition with \(v(-1)=a\) and \(v(1)=b\). For solving Poisson equation.

    • Dirichlet - Homogeneous Dirichlet

    • Neumann - Homogeneous Neumann

    • Biharmonic - Homogeneous Dirichlet and Neumann at both ends

  • dtype (str or np.dtype, optional) – The datatype of physical space (input to forward transforms)

  • quad (str, optional) – Type of quadrature

    • For family=Chebyshev:

      • GL - Chebyshev-Gauss-Lobatto

      • GC - Chebyshev-Gauss

    • For family=Legendre:

      • LG - Legendre-Gauss

      • GL - Legendre-Gauss-Lobatto

    • For family=Laguerre:

      • LG - Laguerre-Gauss

    • For family=Hermite:

      • HG - Hermite-Gauss

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

  • scaled (bool) – Whether to use scaled basis (only Legendre)

  • padding_factor (float, optional) – For padding backward transform (for dealiasing, and only for Fourier)

  • dealias_direct (bool, optional) – Use 2/3-rule dealiasing (only Fourier)

Examples

>>> from shenfun import Basis
>>> F0 = Basis(16, 'F')
>>> C1 = Basis(32, 'C', quad='GC')

shenfun.forms.inner module

This module contains the inner function that computes the weighted inner product.

shenfun.forms.inner.inner(expr0, expr1, output_array=None, level=0)[source]

Return weighted discrete inner product of linear or bilinear form

\[(f, g)_w^N = \sum_{i\in\mathcal{I}}f(x_i) \overline{g}(x_i) w_i \approx \int_{\Omega} g\, \overline{f}\, w\, dx\]

where \(\mathcal{I}=0, 1, \ldots, N, N \in \mathbb{Z}^+\), \(f\) is an expression linear in a TestFunction, and \(g\) is an expression that is linear in TrialFunction or Function, or it is simply an Array (a solution interpolated on the quadrature mesh in physical space). \(w\) is a weight associated with chosen basis, and \(w_i\) are quadrature weights.

If the expressions are created in a multidimensional TensorProductSpace, then the sum above is over all dimensions. In 2D it becomes:

\[(f, g)_w^N = \sum_{i\in\mathcal{I}}\sum_{j\in\mathcal{J}} f(x_i, y_j) \overline{g}(x_i, y_j) w_j w_i\]

where \(\mathcal{J}=0, 1, \ldots, M, M \in \mathbb{Z}^+\).

Parameters
  • expr0, expr1 (Expr, BasisFunction, Array) – or number. Either one can be an expression involving a BasisFunction (TestFunction, TrialFunction or Function) an Array or a number. With expressions (Expr) on a BasisFunction we typically mean terms like div(u) or grad(u), where u is any one of the different types of BasisFunction. One of expr0 or expr1 need to be an expression on a TestFunction. If the second then involves a TrialFunction, a matrix is returned. If one of expr0/expr1 involves a TestFunction and the other one is an expression on a Function, or a plain Array, then a linear form is assembled and a Function is returned. If either expr0 or expr1 is a number (typically 1) or a tuple of numbers for a vector, then the inner product represents a non-weighted integral over the domain. If a single number, then the other expression must be a scalar Array or a Function. If a tuple of numbers, then the Array/Function must be a vector.

  • output_array (Function) – Optional return array for linear form.

  • level (int) – The level of postprocessing for assembled matrices. Applies only to bilinear forms

    • 0 Full postprocessing - diagonal matrices to scale arrays and add equal matrices

    • 1 Diagonal matrices to scale arrays, but don’t add equal matrices

    • 2 No postprocessing, return all assembled matrices

Returns

Function for linear forms.

SparseMatrix for bilinear 1D forms.

TPMatrix or list of TPMatrix for bilinear multidimensional forms.

Number, for non-weighted integral where either one of the arguments is a number.

Return type

Depending on dimensionality and the arguments to the forms

See also

project()

Example

Compute mass matrix of Shen’s Chebyshev Dirichlet basis:

>>> from shenfun import Basis
>>> from shenfun import TestFunction, TrialFunction
>>> SD = Basis(6, 'Chebyshev', bc=(0, 0))
>>> u = TrialFunction(SD)
>>> v = TestFunction(SD)
>>> B = inner(v, u)
>>> 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(abs(B[k]-v) < 1e-7) for k, v in d.items()]
[True, True, True]

shenfun.forms.operators module

This module contains the implementation of operators acting on arguments.

shenfun.forms.operators.div(test)[source]

Return div(test)

Parameters

test (Expr or BasisFunction) – Must be rank > 0 (cannot take divergence of scalar)

shenfun.forms.operators.grad(test)[source]

Return grad(test)

Parameters

test (Expr or BasisFunction)

shenfun.forms.operators.Dx(test, x, k=1)[source]

Return k’th order partial derivative in direction x

Parameters
  • test (Expr or BasisFunction)

  • x (int) – axis to take derivative over

  • k (int) – Number of derivatives

shenfun.forms.operators.curl(test)[source]

Return curl of test

Parameters

test (Expr or BasisFunction)

shenfun.forms.project module

shenfun.forms.project.project(uh, T, output_array=None, fill=True, use_to_ortho=True, use_assign=True)[source]

Project uh to tensor product space T

Find \(u \in T\), such that

\[(u - u_h, v)_w^N = 0 \quad \forall v \in T\]
Parameters
  • uh (Instance of either one of) –

  • T (TensorProductSpace or MixedTensorProductSpace)

  • output_array (Function, optional) – Return array

  • fill (bool, optional) – Whether to fill the output_array with zeros before projection

  • use_to_ortho (bool, optional) – Whether to use fast to_ortho method for projection of Functions to orthogonal space.

  • use_assign (bool, optional) – Whether to use fast assign method for projection of Function to a denser space of the same kind.

Returns

The projection of uh in T

Return type

Function

See also

inner()

Example

>>> import numpy as np
>>> from mpi4py import MPI
>>> from shenfun import Basis, project, TensorProductSpace, Array, \
...     Function, Dx
>>> N = 16
>>> comm = MPI.COMM_WORLD
>>> T0 = Basis(N, 'C')
>>> K0 = Basis(N, 'F', dtype='d')
>>> T = TensorProductSpace(comm, (T0, K0))
>>> uj = Array(T)
>>> uj[:] = np.random.random(uj.shape)
>>> u = Function(T)
>>> u = project(uj, T, output_array=u) # Same as u = T.forward(uj, u)
>>> du = project(Dx(u, 0, 1), T)

Module contents

shenfun.forms.extract_bc_matrices(mats)[source]

Extract boundary matrices from list of mats

Parameters

mats (list of list of :class:`.TPMatrix`es)

Returns

list of boundary matrices.

Return type

list

Note

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