Working on ultraspherical bases

README.rst

@@ -22,7 +22,7 @@ Description
 -----------
 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 <https://fenicsproject.org>`_, but applications are limited to multidimensional tensor product grids, using either Cartesian or curvilinear grids (e.g., but not limited to, polar, cylindrical, spherical or parabolic). The code is parallelized with MPI through the `mpi4py-fft <https://bitbucket.org/mpi4py/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 by using high-order *global* orthogonal basis functions (Fourier, Legendre, Chebyshev first and second kind, Jacobi, Laguerre and Hermite), as opposed to finite element codes that are using low-order *local* basis functions. Efficiency is ensured through vectorization (`Numpy <https://www.numpy.org/>`_), parallelization (`mpi4py <https://bitbucket.org/mpi4py/mpi4py>`_) and by moving critical routines to `Cython <https://cython.org/>`_ or `Numba <https://numba.pydata.org>`_. Shenfun has been used to run turbulence simulations (Direct Numerical Simulations) on thousands of processors on high-performance supercomputers, see the `spectralDNS <https://github.com/spectralDNS/spectralDNS>`_ repository.
+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 by using high-order *global* orthogonal basis functions (Fourier, Legendre, Chebyshev first and second kind, Ultraspherical, Jacobi, Laguerre and Hermite), as opposed to finite element codes that are using low-order *local* basis functions. Efficiency is ensured through vectorization (`Numpy <https://www.numpy.org/>`_), parallelization (`mpi4py <https://bitbucket.org/mpi4py/mpi4py>`_) and by moving critical routines to `Cython <https://cython.org/>`_ or `Numba <https://numba.pydata.org>`_. Shenfun has been used to run turbulence simulations (Direct Numerical Simulations) on thousands of processors on high-performance supercomputers, see the `spectralDNS <https://github.com/spectralDNS/spectralDNS>`_ repository.
 
 The demo folder contains several examples for the Poisson, Helmholtz and Biharmonic equations. For extended documentation and installation instructions see `ReadTheDocs <http://shenfun.readthedocs.org>`_. For interactive demos, see the `jupyter book <https://mikaem.github.io/shenfun-demos>`_. Note that shenfun currently comes with the possibility to use two non-periodic directions (see `biharmonic demo <https://github.com/spectralDNS/shenfun/blob/master/demo/biharmonic2D_2nonperiodic.py>`_), and equations may be solved coupled and implicit (see `MixedPoisson.py <https://github.com/spectralDNS/shenfun/blob/master/demo/MixedPoisson.py>`_).
 

docs/source/gettingstarted.rst

@@ -22,8 +22,9 @@ most important everyday tools are
 	* :func:`.project`
 	* :func:`.FunctionSpace`
 
-A good place to get started is by creating a :func:`.FunctionSpace`. There are seven families of
-function spaces: Fourier, Chebyshev (first and second kind), Legendre, Laguerre, Hermite and Jacobi.
+A good place to get started is by creating a :func:`.FunctionSpace`. There are eight families of
+function spaces: Fourier, Chebyshev (first and second kind), Legendre, Laguerre, Hermite, Ultraspherical
+and Jacobi.
 All spaces are defined on a one-dimensional reference
 domain, with their own basis functions and quadrature points. For example, we have
 the regular orthogonal Chebyshev space :math:`\text{span}\{T_k\}_{k=0}^{N-1}`, where :math:`T_k` is the

docs/source/shenfun.ultraspherical.rst

@@ -0,0 +1,31 @@
+shenfun.ultraspherical package
+==========================
+
+Submodules
+----------
+
+shenfun.ultraspherical.bases module
+-------------------------------
+
+.. automodule:: shenfun.ultraspherical.bases
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+shenfun.ultraspherical.matrices module
+----------------------------------
+
+.. automodule:: shenfun.ultraspherical.matrices
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: shenfun.ultraspherical
+    :members:
+    :undoc-members:
+    :show-inheritance:

setup.py

@@ -97,6 +97,7 @@ def version():
                     "shenfun.hermite",
                     "shenfun.chebyshev",
                     "shenfun.chebyshevu",
+                    "shenfun.ultraspherical",
                     "shenfun.fourier",
                     "shenfun.jacobi",
                     "shenfun.forms",

shenfun/__init__.py

@@ -29,7 +29,7 @@
 """
 #pylint: disable=wildcard-import,no-name-in-module
 
-__version__ = '3.3.0'
+__version__ = '4.0.0'
 __author__ = 'Mikael Mortensen'
 
 import numpy as np
@@ -42,6 +42,7 @@
 from . import hermite
 from . import fourier
 from . import jacobi
+from . import ultraspherical
 from . import matrixbase
 from . import la
 from .coordinates import Coordinates

shenfun/chebyshevu/bases.py

@@ -362,6 +362,8 @@ class Phi1(CompositeBase):
         Boundary conditions at, respectively, x=(-1, 1).
     domain : 2-tuple of floats, optional
         The computational domain
+    scaled : boolean, optional
+        Whether or not to scale basis function n by 1/(n+2)
     dtype : data-type, optional
         Type of input data in real physical space. Will be overloaded when
         basis is part of a :class:`.TensorProductSpace`.
@@ -397,7 +399,7 @@ def stencil_matrix(self, N=None):
         d0, d2 = np.zeros(N), np.zeros(N-2)
         d0[:-2] = sp.lambdify(n, self.b0n)(k[:N-2])
         d2[:] = sp.lambdify(n, self.b2n)(k[:N-2])
-        if self.is_scaled:
+        if self.is_scaled():
             return SparseMatrix({0: d0/(k+2), 2: d2/(k[:-2]+2)}, (N, N))
         return SparseMatrix({0: d0, 2: d2}, (N, N))
 
@@ -674,6 +676,12 @@ class CompactDirichlet(CompositeBase):
     different from 0. In one dimension :math:`\hat{u}_{N-2}=a` and
     :math:`\hat{u}_{N-1}=b`.
 
+    If the parameter `scaled=True`, then the first N-2 basis functions are
+
+    .. math::
+
+        \phi_k &= \frac{U_k}{k+1} - \frac{U_{k+2}}{k+3}, \, k=0, 1, \ldots, N-3, \\
+
     Parameters
     ----------
     N : int, optional
@@ -690,6 +698,8 @@ class CompactDirichlet(CompositeBase):
     dtype : data-type, optional
         Type of input data in real physical space. Will be overloaded when
         basis is part of a :class:`.TensorProductSpace`.
+    scaled : boolean, optional
+        Whether or not to use scaled basis function.
     padding_factor : float, optional
         Factor for padding backward transforms.
     dealias_direct : bool, optional
@@ -703,10 +713,11 @@ class CompactDirichlet(CompositeBase):
 
     """
     def __init__(self, N, quad="GU", bc=(0, 0), domain=(-1, 1), dtype=float,
-                 padding_factor=1, dealias_direct=False, coordinates=None, **kw):
+                 padding_factor=1, dealias_direct=False, coordinates=None,
+                 scaled= False, **kw):
         CompositeBase.__init__(self, N, quad=quad, domain=domain, dtype=dtype, bc=bc,
                                padding_factor=padding_factor, dealias_direct=dealias_direct,
-                               coordinates=coordinates)
+                               coordinates=coordinates, scaled=scaled)
 
     @staticmethod
     def boundary_condition():
@@ -720,8 +731,12 @@ def stencil_matrix(self, N=None):
         N = self.N if N is None else N
         k = np.arange(N)
         d0, d2 = np.zeros(N), np.zeros(N-2)
-        d0[:-2] = 1
-        d2[:] = -(k[:-2]+1)/(k[:-2]+3)
+        if not self.is_scaled():
+            d0[:-2] = 1
+            d2[:] = -(k[:-2]+1)/(k[:-2]+3)
+        else:
+            d0[:-2] = (k[:-2]+3)/(k[:-2]+1)
+            d2[:] = -1
         return SparseMatrix({0: d0, 2: d2}, (N, N))
 
     def slice(self):

shenfun/chebyshevu/matrices.py

@@ -125,7 +125,7 @@ def __init__(self, test, trial, scale=1, measure=1):
         K.shape = (N, N+2)
         K = K.diags('csr')
         B2 = Lmat(2, 0, 2, N, N+2, -half, -half, cn) # B^{(2)_{(2)}}
-        if not test[0].is_scaled:
+        if not test[0].is_scaled():
             k = np.arange(N+2)
             B2 = SparseMatrix({0: (k[:-2]+2)}, (N, N)).diags('csr')*B2
         M = B2 * K.T
@@ -139,7 +139,7 @@ class AP1SDmat(SpectralMatrix):
 
         a_{kj}=(\psi''_j, \phi_k)_w,
 
-    where the test function :math:`\phi_k \in` :class:`.chebyshevu.bases.Phi2`, the trial
+    where the test function :math:`\phi_k \in` :class:`.chebyshevu.bases.Phi1`, the trial
     function :math:`\psi_j \in` :class:`.chebyshev.bases.ShenDirichlet`, and test and
     trial spaces have dimensions of M and N, respectively.
 
@@ -148,7 +148,7 @@ def __init__(self, test, trial, scale=1, measure=1):
         assert isinstance(test[0], P1)
         assert isinstance(trial[0], SD)
         d = {0: -1, 2: 1}
-        if not test[0].is_scaled:
+        if not test[0].is_scaled():
             k = np.arange(test[0].N-2)
             d = {0: -(k+2), 2: k[:-2]+2}
         SpectralMatrix.__init__(self, d, test, trial, scale=scale, measure=measure)
@@ -177,7 +177,7 @@ def __init__(self, test, trial, scale=1, measure=1):
         K.shape = (N, N+2)
         K = K.diags('csr')
         B2 = Lmat(2, 0, 2, N, N+2, -half, -half, cn) # B^{(2)_{(2)}}
-        if not test[0].is_scaled:
+        if not test[0].is_scaled():
             k = np.arange(test[0].N)
             B2 = SparseMatrix({0: (k[:-2]+2)}, (N, N)).diags('csr')*B2
         M = B2 * K.T
@@ -201,7 +201,7 @@ def __init__(self, test, trial, scale=1, measure=1):
         assert isinstance(trial[0], SN)
         k = np.arange(test[0].N-2)
         d = {0: -(k/(k+2))**2, 2: 1}
-        if not test[0].is_scaled:
+        if not test[0].is_scaled():
             d = {0: -k**2/(k+2), 2: k[:-2]+2}
         SpectralMatrix.__init__(self, d, test, trial, scale=scale, measure=measure)
 

shenfun/forms/arguments.py

@@ -258,6 +258,46 @@ def FunctionSpace(N, family='Fourier', bc=None, dtype='d', quad=None,
 
         return B(N, **par)
 
+    elif family.lower() in ('ultraspherical', 'q'):
+        from shenfun import ultraspherical
+        if quad is not None:
+            assert quad == 'QG'
+            par['quad'] = quad
+
+        if scaled is not None:
+            assert isinstance(scaled, bool)
+            par['scaled'] = scaled
+
+        bases = defaultdict(lambda: ultraspherical.bases.Generic,
+                            {
+                                '': ultraspherical.bases.Orthogonal,
+                                'LDRD': ultraspherical.bases.CompactDirichlet,
+                                'LNRN': ultraspherical.bases.CompactNeumann,
+                                'LDLNRDRN': ultraspherical.bases.Phi2,
+                                'LDLNLN2RDRNRN2': ultraspherical.bases.Phi3,
+                                'LDLNLN2N3RDRNRN2N3': ultraspherical.bases.Phi4
+                            })
+
+        if basis is not None:
+            assert isinstance(basis, str)
+            B = getattr(ultraspherical.bases, basis)
+            if isinstance(bc, tuple):
+                par['bc'] = bc
+        else:
+            if isinstance(bc, (str, tuple, dict)):
+                domain = (-1, 1) if domain is None else domain
+                bcs = BoundaryConditions(bc, domain=domain)
+                key = bcs.orderednames()
+                par['bc'] = bcs
+
+            elif bc is None:
+                key = ''
+            else:
+                raise NotImplementedError
+            B = bases[''.join(key)]
+
+        return B(N, **par)
+
     elif family.lower() in ('laguerre', 'la'):
         from shenfun import laguerre
         if quad is not None:

shenfun/jacobi/bases.py

@@ -34,7 +34,6 @@
 from shenfun.spectralbase import SpectralBase, Transform, islicedict, \
     slicedict, getCompositeBase, BoundaryConditions
 from shenfun.matrixbase import SparseMatrix
-from .recursions import n
 
 try:
     import quadpy
@@ -45,9 +44,6 @@
     has_quadpy = False
     mp = None
 
-mode = config['bases']['jacobi']['mode']
-mode = mode if has_quadpy else 'numpy'
-
 xp = sp.Symbol('x', real=True)
 m, n, k = sp.symbols('m,n,k', real=True, integer=True)
 
@@ -64,7 +60,6 @@
            'Generic',
            'BCBase',
            'BCGeneric',
-           'mode',
            'has_quadpy',
            'mp']
 
@@ -136,6 +131,7 @@ def points_and_weights(self, N=None, map_true_domain=False, weighted=True, **kw)
         return points, weights
 
     def mpmath_points_and_weights(self, N=None, map_true_domain=False, weighted=True, **kw):
+        mode = config['bases']['jacobi']['mode']
         if mode == 'numpy' or not has_quadpy:
             return self.points_and_weights(N=N, map_true_domain=map_true_domain, weighted=weighted, **kw)
         if N is None:
@@ -148,6 +144,7 @@ def mpmath_points_and_weights(self, N=None, map_true_domain=False, weighted=True
         return points, weights
 
     def jacobi(self, x, alpha, beta, N):
+        mode = config['bases']['jacobi']['mode']
         V = np.zeros((x.shape[0], N))
         if mode == 'numpy':
             for n in range(N):
@@ -219,6 +216,7 @@ def bnd_values(k=0, alpha=0, beta=0):
         return bnd_values(alpha, beta, k=k)
 
     def evaluate_basis(self, x, i=0, output_array=None):
+        mode = config['bases']['jacobi']['mode']
         x = np.atleast_1d(x)
         if output_array is None:
             output_array = np.zeros(x.shape)
@@ -230,9 +228,10 @@ def evaluate_basis(self, x, i=0, output_array=None):
         return output_array
 
     def evaluate_basis_derivative(self, x=None, i=0, k=0, output_array=None):
+        mode = config['bases']['jacobi']['mode']
         if x is None:
             x = self.points_and_weights(mode=mode)[0]
-        #x = np.atleast_1d(x)
+        x = np.atleast_1d(x)
         if output_array is None:
             output_array = np.zeros(x.shape, dtype=self.dtype)
 
@@ -245,6 +244,7 @@ def evaluate_basis_derivative(self, x=None, i=0, k=0, output_array=None):
         return output_array
 
     def evaluate_basis_derivative_all(self, x=None, k=0, argument=0):
+        mode = config['bases']['jacobi']['mode']
         if x is None:
             x = self.mpmath_points_and_weights(mode=mode)[0]
         if mode == 'numpy':
@@ -380,7 +380,7 @@ class Phi2(CompositeBase):
     The basis functions :math:`\phi_k` for :math:`k=0, 1, \ldots, N-5` are
 
     .. math::
-        \phi_k &= \frac{(1-x^2)^2}{h^{(2,\alpha,\beta)}_{k+2}} \frac{d^$P^{(\alpha,\beta)}_{k+4}}{dx},
+        \phi_k &= \frac{(1-x^2)^2}{h^{(2,\alpha,\beta)}_{k+2}} \frac{d^2P^{(\alpha,\beta)}_{k+2}}{dx^2},
 
     where
 
@@ -440,6 +440,7 @@ def __init__(self, N, quad="JG", bc=(0, 0, 0, 0), domain=(-1., 1.), dtype=float,
         self.b0n = 2**(-a - b + 1)*sp.gamma(n + 1)*sp.gamma(a + b + n + 3)/((a + b + 2*n + 2)*(a + b + 2*n + 3)*(a + b + 2*n + 4)*sp.gamma(a + n + 1)*sp.gamma(b + n + 1))
         self.b2n = 2**(-a - b + 1)*(a + b + 2*n + 5)*(a**2 - 4*a*b - 2*a*n - 5*a + b**2 - 2*b*n - 5*b - 2*n**2 - 10*n - 12)*sp.gamma(n + 3)*sp.gamma(a + b + n + 3)/((a + b + 2*n + 3)*(a + b + 2*n + 4)*(a + b + 2*n + 6)*(a + b + 2*n + 7)*sp.gamma(a + n + 3)*sp.gamma(b + n + 3))
         self.b4n = 2**(-a - b + 1)*sp.gamma(n + 5)*sp.gamma(a + b + n + 3)/((a + b + 2*n + 6)*(a + b + 2*n + 7)*(a + b + 2*n + 8)*sp.gamma(a + n + 3)*sp.gamma(b + n + 3))
+
         if self.alpha != self.beta:
             self.b1n = 2**(-a - b + 2)*(a - b)*sp.gamma(n + 2)*sp.gamma(a + b + n + 3)/((a + b + 2*n + 2)*(a + b + 2*n + 4)*(a + b + 2*n + 6)*sp.gamma(a + n + 2)*sp.gamma(b + n + 2))
             self.b3n = -2**(-a - b + 2)*(a - b)*sp.gamma(n + 4)*sp.gamma(a + b + n + 3)/((a + b + 2*n + 4)*(a + b + 2*n + 6)*(a + b + 2*n + 8)*sp.gamma(a + n + 3)*sp.gamma(b + n + 3))
@@ -458,13 +459,13 @@ def stencil_matrix(self, N=None):
             return self._stencil_matrix[N]
         k = np.arange(N)
         d0, d2, d4 = np.zeros(N), np.zeros(N-2), np.zeros(N-4)
-        d0[:-4] = sp.lambdify(n, self.b0n)(k[:N-4])
-        d2[:-2] = sp.lambdify(n, self.b2n)(k[:N-4])
-        d4[:] = sp.lambdify(n, self.b4n)(k[:N-4])
+        d0[:-4] = sp.lambdify(n, sp.simplify(self.b0n))(k[:N-4])
+        d2[:-2] = sp.lambdify(n, sp.simplify(self.b2n))(k[:N-4])
+        d4[:] = sp.lambdify(n, sp.simplify(self.b4n))(k[:N-4])
         if self.alpha != self.beta:
             d1, d3 = np.zeros(N-1), np.zeros(N-3)
-            d1[:-3] = sp.lambdify(n, self.b1n)(k[:N-4])
-            d3[:-1] = sp.lambdify(n, self.b3n)(k[:N-4])
+            d1[:-3] = sp.lambdify(n, sp.simplify(self.b1n))(k[:N-4])
+            d3[:-1] = sp.lambdify(n, sp.simplify(self.b3n))(k[:N-4])
             self._stencil_matrix[N] = SparseMatrix({0: d0, 1: d1, 2: d2, 3: d3, 4: d4}, (N, N))
         else:
             self._stencil_matrix[N] = SparseMatrix({0: d0, 2: d2, 4: d4}, (N, N))
@@ -584,7 +585,7 @@ class Phi4(CompositeBase):
     The basis functions :math:`\phi_k` for :math:`k=0, 1, \ldots, N-9` are
 
     .. math::
-        \phi_k &= \frac{(1-x^2)^4}{h^{(4,\alpha,\beta)}_{k+4}} \frac{d^4P^{(\alpha,\beta)}_{k+4}}{dx}, \\
+        \phi_k &= \frac{(1-x^2)^4}{h^{(4,\alpha,\beta)}_{k+4}} \frac{d^4P^{(\alpha,\beta)}_{k+4}}{dx^4}, \\
 
     where