Added linear transformation for non-standard domains in Chebychev and Ultraspherical methods (#536)

* Added linear transformation for non-standard domains in Chebychev and
Ultraspherical methods

Co-authored-by: tlunet

* Removed clutter in tests

Author: brownbaerchen

pySDC/helpers/spectral_helper.py

@@ -194,8 +194,9 @@ def __init__(self, *args, transform_type='fft', x0=-1, x1=1, **kwargs):
             x0 (float): Coordinate of left boundary. Note that only -1 is currently implented
             x1 (float): Coordinate of right boundary. Note that only +1 is currently implented
         """
-        assert x0 == -1
-        assert x1 == 1
+        # need linear transformation y = ax + b with a = (x1-x0)/2 and b = (x1+x0)/2
+        self.lin_trf_fac = (x1 - x0) / 2
+        self.lin_trf_off = (x1 + x0) / 2
         super().__init__(*args, x0=x0, x1=x1, **kwargs)
         self.transform_type = transform_type
 
@@ -214,7 +215,7 @@ def get_1dgrid(self):
         Returns:
             numpy.ndarray: 1D grid
         '''
-        return self.xp.cos(np.pi / self.N * (self.xp.arange(self.N) + 0.5))
+        return self.lin_trf_fac * self.xp.cos(np.pi / self.N * (self.xp.arange(self.N) + 0.5)) + self.lin_trf_off
 
     def get_wavenumbers(self):
         """Get the domain in spectral space"""
@@ -306,7 +307,7 @@ def get_integration_matrix(self, lbnd=0):
                 (n / (2 * (self.xp.arange(self.N) + 1)))[1::2]
                 * (-1) ** (self.xp.arange(self.N // 2))
                 / (np.append([1], self.xp.arange(self.N // 2 - 1) + 1))
-            )
+            ) * self.lin_trf_fac
         else:
             raise NotImplementedError(f'This function allows to integrate only from x=0, you attempted from x={lbnd}.')
         return S
@@ -327,7 +328,7 @@ def get_differentiation_matrix(self, p=1):
                 D[k, j] = 2 * j * ((j - k) % 2)
 
         D[0, :] /= 2
-        return self.sparse_lib.csc_matrix(self.xp.linalg.matrix_power(D, p))
+        return self.sparse_lib.csc_matrix(self.xp.linalg.matrix_power(D, p)) / self.lin_trf_fac**p
 
     def get_norm(self, N=None):
         '''
@@ -565,7 +566,7 @@ def get_differentiation_matrix(self, p=1):
         xp = self.xp
         N = self.N
         l = p
-        return 2 ** (l - 1) * factorial(l - 1) * sp.diags(xp.arange(N - l) + l, offsets=l)
+        return 2 ** (l - 1) * factorial(l - 1) * sp.diags(xp.arange(N - l) + l, offsets=l) / self.lin_trf_fac**p
 
     def get_S(self, lmbda):
         """
@@ -647,8 +648,10 @@ def get_integration_matrix(self):
         Returns:
             sparse integration matrix
         """
-        return self.sparse_lib.diags(1 / (self.xp.arange(self.N - 1) + 1), offsets=-1) @ self.get_basis_change_matrix(
-            p_out=1, p_in=0
+        return (
+            self.sparse_lib.diags(1 / (self.xp.arange(self.N - 1) + 1), offsets=-1)
+            @ self.get_basis_change_matrix(p_out=1, p_in=0)
+            * self.lin_trf_fac
         )
 
     def get_integration_constant(self, u_hat, axis):

pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py

@@ -79,7 +79,7 @@ def test_differentiation_matrix(N):
     from pySDC.helpers.spectral_helper import ChebychevHelper
 
     cheby = ChebychevHelper(N)
-    x = np.cos(np.pi / N * (np.arange(N) + 0.5))
+    x = cheby.get_1dgrid()
     coeffs = np.random.random(N)
     norm = cheby.get_norm()
 
@@ -91,6 +91,36 @@ def test_differentiation_matrix(N):
     assert np.allclose(exact, du)
 
 
+@pytest.mark.base
+@pytest.mark.parametrize('N', [4, 7, 32])
+@pytest.mark.parametrize('x0', [-1, 0])
+@pytest.mark.parametrize('x1', [0.789, 1])
+@pytest.mark.parametrize('p', [1, 2])
+def test_differentiation_non_standard_domain_size(N, x0, x1, p):
+    import numpy as np
+    import scipy
+    from pySDC.helpers.spectral_helper import ChebychevHelper
+
+    cheby = ChebychevHelper(N, x0=x0, x1=x1)
+    x = cheby.get_1dgrid()
+    assert all(x > x0)
+    assert all(x < x1)
+
+    coeffs = np.random.random(N)
+    u = np.polynomial.Chebyshev(coeffs)(x)
+    u_hat = cheby.transform(u)
+    du_exact = np.polynomial.Chebyshev(coeffs).deriv(p)(x)
+    du_hat_exact = cheby.transform(du_exact)
+
+    D = cheby.get_differentiation_matrix(p)
+
+    du_hat = D @ u_hat
+    du = cheby.itransform(du_hat)
+
+    assert np.allclose(du_hat_exact, du_hat), np.linalg.norm(du_hat_exact - du_hat)
+    assert np.allclose(du, du_exact), np.linalg.norm(du_exact - du)
+
+
 @pytest.mark.base
 @pytest.mark.parametrize('N', [4, 32])
 def test_integration_matrix(N):
@@ -121,7 +151,8 @@ def test_transform(N, d, transform_type):
     cheby = ChebychevHelper(N, transform_type=transform_type)
     u = np.random.random((d, N))
     norm = cheby.get_norm()
-    x = cheby.get_1dgrid()
+    x = (cheby.get_1dgrid() * cheby.lin_trf_fac + cheby.lin_trf_off) * cheby.lin_trf_fac + cheby.lin_trf_off
+    x = (cheby.get_1dgrid() - cheby.lin_trf_off) / cheby.lin_trf_fac
 
     itransform = cheby.itransform(u, axis=-1).real
 
@@ -370,15 +401,3 @@ def test_tau_method2D_diffusion(nz, nx, bc_val, plotting=False):
         assert np.allclose(
             polys[i](z), sol[0, i, :]
         ), f'Solution is incorrectly transformed back to real space at x={x[i]}'
-
-
-if __name__ == '__main__':
-    test_differentiation_matrix(4, 'T2U')
-    # test_transform(6, 1, 'fft')
-    # test_tau_method('T2U', -1.0, N=4, bc_val=3.0)
-    # test_tau_method2D('T2T', -1, nx=2**7, nz=2**6, bc_val=4.0, plotting=True)
-    # test_integration_matrix(5, 'T2U')
-    # test_integration_matrix2D(2**0, 2**2, 'T2U', 'z')
-    # test_differentiation_matrix2D(2**7, 2**7, 'T2U', 'mixed')
-    # test_integration_BC(6)
-    # test_filter(12, 2, 5, 'T2U')

pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py

@@ -21,6 +21,37 @@ def test_differentiation_matrix(N, p):
     assert np.allclose(exact, du)
 
 
+@pytest.mark.base
+@pytest.mark.parametrize('N', [4, 7, 32])
+@pytest.mark.parametrize('x0', [-1, 0])
+@pytest.mark.parametrize('x1', [0.789, 1])
+@pytest.mark.parametrize('p', [1, 2])
+def test_differentiation_non_standard_domain_size(N, x0, x1, p):
+    import numpy as np
+    import scipy
+    from pySDC.helpers.spectral_helper import UltrasphericalHelper
+
+    helper = UltrasphericalHelper(N, x0=x0, x1=x1)
+    x = helper.get_1dgrid()
+    assert all(x > x0)
+    assert all(x < x1)
+
+    coeffs = np.random.random(N)
+    u = np.polynomial.Chebyshev(coeffs)(x)
+    u_hat = helper.transform(u)
+    du_exact = np.polynomial.Chebyshev(coeffs).deriv(p)(x)
+    du_hat_exact = helper.transform(du_exact)
+
+    D = helper.get_differentiation_matrix(p)
+    Q = helper.get_basis_change_matrix(p_in=p, p_out=0)
+
+    du_hat = Q @ D @ u_hat
+    du = helper.itransform(du_hat)
+
+    assert np.allclose(du_hat_exact, du_hat), np.linalg.norm(du_hat_exact - du_hat)
+    assert np.allclose(du, du_exact), np.linalg.norm(du_exact - du)
+
+
 @pytest.mark.base
 @pytest.mark.parametrize('N', [4, 7, 32])
 def test_integration(N):
@@ -94,9 +125,3 @@ def test_poisson_problem(N, deg, Dirichlet_recombination):
 
     assert np.allclose(u_hat[deg + 3 :], 0)
     assert np.allclose(u_exact, u)
-
-
-if __name__ == '__main__':
-    # test_differentiation_matrix(6, 2)
-    test_poisson_problem(6, 1, True)
-    # test_integration()