shenfun.utilities package

Submodules

shenfun.utilities.findbasis module

get_bc_basis(bcs, family, alpha=None, beta=None, gn=1)[source]

Return boundary basis satisfying bcs.

Parameters:

  • bcs (dict) – The boundary conditions in dictionary form, see BoundaryConditions.

  • family (str) – Choose one of

    • Chebyshev

    • Chebyshevu

    • Legendre

    • Ultraspherical

    • Jacobi

    • Laguerre

  • alpha, beta (numbers, optional) – The Jacobi parameters, used only for Jacobi or Ultraspherical

  • gn (Scaling function for Jacobi polynomials)

get_stencil_matrix(bcs, family, alpha=None, beta=None, gn=1)[source]

Return stencil matrix.

Return the narrowest possible stencil matrix satisfying the homogeneous boundary conditions in bcs, and with leading coefficient equal to one.

Parameters:

  • bcs (dict) – The boundary conditions in dictionary form, see BoundaryConditions.

  • family (str) – Choose one of

    • Chebyshev

    • Chebyshevu

    • Legendre

    • Ultraspherical

    • Jacobi

    • Laguerre

  • alpha, beta (numbers, optional) – The Jacobi parameters, used only for Jacobi or Ultraspherical

  • gn (Scaling function for Jacobi polynomials)

shenfun.utilities.integrators module

Module for some integrators.

Integrators are set up to solve equations like

\[\frac{\partial u}{\partial t} = L u + N(u)\]

where \(u\) is the solution, \(L\) is a linear operator and \(N(u)\) is the nonlinear part of the right hand side.

There are two kinds of integrators, or time steppers. The first are ment to be subclassed and used one at the time. These are

  • IRK3: Implicit third order Runge-Kutta

  • RK4: Runge-Kutta fourth order

  • ETD: Exponential time differencing Euler method

  • ETDRK4: Exponential time differencing Runge-Kutta fourth order

See, e.g., H. Montanelli and N. Bootland “Solving periodic semilinear PDEs in 1D, 2D and 3D with exponential integrators”, https://arxiv.org/pdf/1604.08900.pdf

The second kind is ment to be used for systems of equations, one class instance for each equation. These are mainly IMEX Runge Kutta integrators:

  • IMEXRK3

  • IMEXRK111

  • IMEXRK222

  • IMEXRK443

See, e.g., https://github.com/spectralDNS/shenfun/blob/master/demo/ChannelFlow for an example of use.

Note

RK4, ETD and ETDRK4 can only be used with Fourier function spaces, as they assume all matrices are diagonal.

class BackwardEuler(T, L=None, N=None, update=None, **params)[source]

Bases: IntegratorBase

First order backward Euler

Parameters:

  • T (TensorProductSpace)

  • L (function of TrialFunction(T)) – To compute linear part of right hand side

  • N (function) – To compute nonlinear part of right hand side

  • update (function) – To be called at the end of a timestep

  • params (dictionary) – Any relevant keyword arguments

compute_rhs(u, u_hat, dU, dU1)[source]
setup(dt)[source]

Set up solver

solve(u, u_hat, dt, trange)[source]

Integrate forward in time

Parameters:

  • u (array) – The solution array in physical space

  • u_hat (array) – The solution array in spectral space

  • dt (float) – Timestep

  • trange (two-tuple) – Time and end time

class ETD(T, L=None, N=None, update=None, **params)[source]

Bases: IntegratorBase

Exponential time differencing Euler method

H. Montanelli and N. Bootland “Solving periodic semilinear PDEs in 1D, 2D and 3D with exponential integrators”, https://arxiv.org/pdf/1604.08900.pdf

Parameters:

  • T (TensorProductSpace)

  • L (function) – To compute linear part of right hand side

  • N (function) – To compute nonlinear part of right hand side

  • update (function) – To be called at the end of a timestep

  • params (dictionary) – Any relevant keyword arguments

setup(dt)[source]

Set up ETD ODE solver

solve(u, u_hat, dt, trange)[source]

Integrate forward in time

Parameters:

  • u (array) – The solution array in physical space

  • u_hat (array) – The solution array in spectral space

  • dt (float) – Timestep

  • trange (two-tuple) – Time and end time

class ETDRK4(T, L=None, N=None, update=None, **params)[source]

Bases: IntegratorBase

Exponential time differencing Runge-Kutta 4’th order method

H. Montanelli and N. Bootland “Solving periodic semilinear PDEs in 1D, 2D and 3D with exponential integrators”, https://arxiv.org/pdf/1604.08900.pdf

Parameters:

  • T (TensorProductSpace)

  • L (function) – To compute linear part of right hand side

  • N (function) – To compute nonlinear part of right hand side

  • update (function) – To be called at the end of a timestep

  • params (dictionary) – Any relevant keyword arguments

setup(dt)[source]

Set up ETDRK4 ODE solver

solve(u, u_hat, dt, trange)[source]

Integrate forward in time

Parameters:

  • u (array) – The solution array in physical space

  • u_hat (array) – The solution array in spectral space

  • dt (float) – Timestep

  • trange (two-tuple) – Time and end time

class IMEXRK011(v, u, L, N, dt, solver=None, name='U-equation', latex=None)[source]

Bases: PDEIMEXRK

classmethod steps()[source]
stages()[source]
class IMEXRK111(v, u, L, N, dt, solver=None, name='U-equation', latex=None)[source]

Bases: PDEIMEXRK

classmethod steps()[source]
stages()[source]
class IMEXRK222(v, u, L, N, dt, solver=None, name='U-equation', latex=None)[source]

Bases: PDEIMEXRK

classmethod steps()[source]
stages()[source]
class IMEXRK3(v, u, L, N, dt, solver=None, name='U-equation', latex=None)[source]

Bases: object

Solve partial differential equations of the form

\[\frac{\partial u}{\partial t} = N+Lu, \quad (1)\]

where \(N\) is a nonlinear term and \(L\) is a linear operator.

This is a third order accurate implicit Runge-Kutta solver

Parameters:

  • v (TestFunction)

  • u (Expr or Function) – Representing \(u\) in (1) If an Expr, then its basis must be a Function The Function will then hold the solution. If \(u\) is a Function, then this will hold the solution.

  • L (Linear operator) – Operates on \(u\)

  • N (Expr or sequence of Expr) – Nonlinear terms

  • dt (number) – Time step

  • solver (Linear solver, optional)

  • name (str, optional)

classmethod steps()[source]
assemble()[source]
compute_rhs(rk=0)[source]
solve_step(rk)[source]
stages()[source]
class IMEXRK443(v, u, L, N, dt, solver=None, name='U-equation', latex=None)[source]

Bases: PDEIMEXRK

classmethod steps()[source]
stages()[source]
class IRK3(T, L=None, N=None, update=None, **params)[source]

Bases: IntegratorBase

Third order implicit Runge Kutta

Parameters:

  • T (TensorProductSpace)

  • L (function of TrialFunction(T)) – To compute linear part of right hand side

  • N (function) – To compute nonlinear part of right hand side

  • update (function) – To be called at the end of a timestep

  • params (dictionary) – Any relevant keyword arguments

compute_rhs(u, u_hat, dU, dU1, rk)[source]
setup(dt)[source]

Set up solver

solve(u, u_hat, dt, trange)[source]

Integrate forward in time

Parameters:

  • u (array) – The solution array in physical space

  • u_hat (array) – The solution array in spectral space

  • dt (float) – Timestep

  • trange (two-tuple) – Time and end time

class RK4(T, L=None, N=None, update=None, **params)[source]

Bases: IntegratorBase

Regular 4’th order Runge-Kutta integrator

Parameters:

  • T (TensorProductSpace)

  • L (function) – To compute linear part of right hand side

  • N (function) – To compute nonlinear part of right hand side

  • update (function) – To be called at the end of a timestep

  • params (dictionary) – Any relevant keyword arguments

setup(dt)[source]

Set up RK4 ODE solver

solve(u, u_hat, dt, trange)[source]

Integrate forward in end_time

Parameters:

  • u (array) – The solution array in physical space

  • u_hat (array) – The solution array in spectral space

  • dt (float) – Timestep

  • trange (two-tuple) – Time and end time

shenfun.utilities.lagrangian_particles module

class LagrangianParticles(points, dt, u_hat)[source]

Bases: object

Class for tracking Lagrangian particles

Parameters:

  • points (array) – Initial location of particles. (D, N) array, with N particles in D dimensions

  • dt (float) – Time step

  • u_hat (Function) – Spectral Galerkin Function for the Eulerian velocity

rhs()[source]
step()[source]

Module contents

Module for implementing helper functions.

class CachedArrayDict[source]

Bases: MutableMapping

Dictionary for caching Numpy arrays (work arrays)

Example

>>> import numpy as np
>>> from shenfun.utilities import CachedArrayDict
>>> work = CachedArrayDict()
>>> a = np.ones((3, 4), dtype=int)
>>> w = work[(a, 0, True)] # create work array with shape as a
>>> print(w.shape)
(3, 4)
>>> print(w)
[[0 0 0 0]
 [0 0 0 0]
 [0 0 0 0]]
>>> w2 = work[(a, 1, True)] # Get different(note 1!) array of same shape/dtype
values() → an object providing a view on D's values[source]
Lambda(x)

Return

\[\Lambda(x) =\]

rac{Gamma(x+ rac{1}{2})}{Gamma(x+1)}

x : array of floats

cleanup(d, destroy_subcomm=True)[source]

Destroy MPI transfer objects in dictionary d

Parameters:

  • d (dictionary or tuple) – dictionary containing variables (e.g., locals(), vars()) or tuple containing objects with a destroy method

  • destroy_subcomm (bool, optional) – Whether to destroy communicators

clenshaw_curtis1D(u, quad='GC')[source]

Clenshaw-Curtis integration in 1D

cross(c, a, b)[source]

Cross product c = a x b

Parameters:

  • c (Array)

  • a (Array)

  • b (Array)

Returns:

c

Return type:

Array

dot(u, v, output_array=None, forward_output=True)[source]

Return dot product of u and v

Parameters:

  • u (Function or Array)

  • v (Function or Array)

  • output_array (Function or Array) – Must be of correct shape

  • forward_output (bool, optional) – Return the result as a Function if True, and an Array if False.

Note

This function uses the 3/2-rule for dealiasing all spaces if the input arrays are spectral Functions. The returned Function is projected to an orthogonal space, without boundary conditions.

dx(u, weighted=False)[source]

Compute integral of u over domain

\[\int_{\Omega} u dx\]
Parameters:

u (Array) – The Array to integrate

Note

This function assumes a standard reference domain. If the domain is not standard, then modify the result accordingly.

get_bc_basis(bcs, family, alpha=None, beta=None, gn=1)[source]

Return boundary basis satisfying bcs.

Parameters:

  • bcs (dict) – The boundary conditions in dictionary form, see BoundaryConditions.

  • family (str) – Choose one of

    • Chebyshev

    • Chebyshevu

    • Legendre

    • Ultraspherical

    • Jacobi

    • Laguerre

  • alpha, beta (numbers, optional) – The Jacobi parameters, used only for Jacobi or Ultraspherical

  • gn (Scaling function for Jacobi polynomials)

get_stencil_matrix(bcs, family, alpha=None, beta=None, gn=1)[source]

Return stencil matrix.

Return the narrowest possible stencil matrix satisfying the homogeneous boundary conditions in bcs, and with leading coefficient equal to one.

Parameters:

  • bcs (dict) – The boundary conditions in dictionary form, see BoundaryConditions.

  • family (str) – Choose one of

    • Chebyshev

    • Chebyshevu

    • Legendre

    • Ultraspherical

    • Jacobi

    • Laguerre

  • alpha, beta (numbers, optional) – The Jacobi parameters, used only for Jacobi or Ultraspherical

  • gn (Scaling function for Jacobi polynomials)

integrate_sympy(f, d)[source]

Exact definite integral using sympy

Try to convert expression f to a polynomial before integrating.

See sympy issue https://github.com/sympy/sympy/pull/18613 to why this is needed. Poly().integrate() is much faster than sympy.integrate() when applicable.

Parameters:

  • f (sympy expression)

  • d (3-tuple) – First item the symbol, next two the lower and upper integration limits

mayavi_show()[source]

Return show function that updates the mayavi figure in the background.

outer(a, b, c)[source]

Return outer product $c_{i,j} = a_i b_j$

Parameters:

  • a (Array of shape (N, …))

  • b (Array of shape (N, …))

  • c (Array of shape (N*N, …))

  • The outer product is taken over the first index of a and b,

  • for all remaining indices.

quiver3D(u, mesh=None, wrapaxes=(), slices=None, fig=None, kind='quadrature', **kw)[source]
Parameters:

  • u (Array)

  • mesh (list of Cartesian meshes [x, y, z], optional)

  • backend (str, optional) – plotly or mayavi

  • wrapaxes (sequence of integers, optional) – For domains that wrap around, extend mesh/u by one in given direction

  • fig (None, optional) – If None create a new figure

  • slices (None or sequence of slices, optional) – If only part of u should be plotted

  • kind (str, optional) –

    • ‘quadrature’ - Use quadrature mesh

    • ‘uniform’ - Use uniform mesh

  • kw (Keyword arguments, optional) – Arguments to mlab.quiver3d, for example ‘scale_factor’, ‘color’ or ‘’mode

reset_profile(prof)[source]

Reset profiler for kernprof

Parameters:

prof (The profiler)

scalar_product(v, f, output_array=None, assemble='exact')[source]

Return scalar product

\[(v, f)_w\]
Parameters:

  • v (TestFunction)

  • f (Sympy function)

  • output_array (Function)

  • assemble (str, optional) –

    • ‘exact’

    • ‘adaptive’

Note

The computed scalar product is not compensated for a non-standard domain size.

>>> from shenfun import inner, FunctionSpace, TestFunction
>>> import sympy as sp
>>> C = FunctionSpace(4, 'C')
>>> v = TestFunction(C)
>>> x = sp.Symbol('x', real=True)
>>> f = x**2
>>> inner(v, f, assemble='exact')
Function([1.57079633, 0.        , 0.78539816, 0.        ])
surf3D(u, mesh=None, backend='plotly', wrapaxes=(), slices=None, fig=None, **kw)[source]

Plot surface embedded in 3D

Parameters:

  • u (Function or Array)

  • mesh (list, str or None, optional) –

    • list of Cartesian meshes [x, y, z]

    • ‘quadrature’ - use quadrature mesh

    • ‘uniform’ - use uniform mesh

  • backend (str, optional) – plotly or mayavi

  • wrapaxes (sequence of integers, optional) – For domains that wrap around, extend mesh/u by one in given direction

  • slices (None or sequence of slices, optional) – If only part of u should be plotted

  • fig (Figure instance, optional)

  • kw (Keyword arguments, optional) – Used by plotly Surface. Possibly colorscale.

wrap_periodic(xs, axes=())[source]

Return arrays wrapped around periodically

Parameters:

  • xs (array or sequence of arrays)

  • axes (sequence of integers, optional) – Extend arrays in xs by one in direction given by wrap