Other functionality of shenfun

from shenfun import Dx
v = TestFunction(W)
du = Dx(v, 0, 1)

where the first argument to Dx is a basis function, the second (integer) is the axis to take the derivative over, and the third (integer) is the number of derivatives, e.g., $$ \begin{equation} \frac{\partial^2 v}{\partial y^2} = \text{Dx(v, 1, 2)}. \notag \tag{61} \end{equation} $$ The operator can be nested. To compute \( \frac{\partial^2 u}{\partial x \partial y} \) one may do

v = TestFunction(W)
du = Dx(Dx(v, 0, 1), 1, 1)

The operators work on TestFunctions, TrialFunctions or Functions, where only the last actually contain any data, because a Function is used to store the solution. Once a solution has been found, one may also manipulate it further using project in combination with operators on Functions. For example, to compute \( \partial u / \partial x \) of the solution to the biharmonic problem, one can do

K0 = Basis(N[0])
W0 = TensorProductSpace(comm, (K0, K1, K2))
du_hat = project(Dx(u, 0, 1), W0, uh_hat=u_hat)
du = Function(W0, False)
du = W0.backward(du_hat , du)

Note that we are here using a regular Chebyshev space (Basis) instead of the biharmonic, to avoid enforcing erroneous boundary conditions on the solution. We could in this case also, with advantage, have chosen a Dirichlet space, since the derivative of the biharmonic problem is known to be zero on the edges of the domain (at \( x=\pm 1 \)).

All problems considered thus far have been scalar valued. With shenfun there is also some functionality for working with vector equations. To this end, there is a class called VectorTensorProductSpace, and there is an additional operator, curl, that can only be used on vectors:

from shenfun.chebyshev.bases import ShenBiharmonicBasis
from shenfun.fourier.bases import R2CBasis, C2CBasis
from shenfun import Function, TensorProductSpace, VectorTensorProductSpace , curl
from shenfun import inner, curl, TestFunction
import numpy as np
from mpi4py import MPI
comm = MPI.COMM_WORLD

N = (32, 33, 34)
K0 = ShenBiharmonicBasis(N[0])
K1 = C2CBasis(N[1])
K2 = R2CBasis(N[2])

T = TensorProductSpace(comm, (K0, K1, K2))
Tk = VectorTensorProductSpace([T, T, T])

v = TestFunction(Tk)
u_ = Function(Tk, False)
u_[:] = np.random.random(u_.shape)
u_hat = Function(Tk)
u_hat = Tk.forward(u_, u_hat)
w_hat = inner(v, curl(u_), uh_hat=u_hat)

Vector equations have very similar form as scalar equations, but at the moment of writing the different equation components cannot be implicitly coupled.

• ← Prev
• Next →