diff --git a/pySDC/helpers/spectral_helper.py b/pySDC/helpers/spectral_helper.py index cd175bb738722db435eaf0649bec825ed073b5d0..d74aab8733f256272f0cb99a5ea6e1fd13fe3a34 100644 --- a/pySDC/helpers/spectral_helper.py +++ b/pySDC/helpers/spectral_helper.py @@ -194,8 +194,9 @@ class ChebychevHelper(SpectralHelper1D): 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 @@ class ChebychevHelper(SpectralHelper1D): 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 @@ class ChebychevHelper(SpectralHelper1D): (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 @@ class ChebychevHelper(SpectralHelper1D): 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 @@ class UltrasphericalHelper(ChebychevHelper): 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 @@ class UltrasphericalHelper(ChebychevHelper): 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): diff --git a/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py b/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py index 37a1ee7175e1a77946adcc5701c08b7c161314da..5052e5908c086b3e1db51506f06409b7723108cb 100644 --- a/pySDC/tests/test_helpers/test_spectral_helper_1d_chebychev.py +++ b/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') diff --git a/pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py b/pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py index d5ba10afbf18fab9eb031bd68c35597bb6601106..22d5fb5001d524918b7e2ebab8abd340fbc15080 100644 --- a/pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py +++ b/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()