Skip to content

Commit

Permalink
Added some problem documentation and renamed the guy
Browse files Browse the repository at this point in the history
  • Loading branch information
brownbaerchen committed Sep 6, 2024
1 parent 821829d commit 6fbfa72
Show file tree
Hide file tree
Showing 7 changed files with 98 additions and 68 deletions.
34 changes: 17 additions & 17 deletions pySDC/helpers/spectral_helper.py
Original file line number Diff line number Diff line change
Expand Up @@ -117,12 +117,12 @@ def get_1dgrid(self):
raise NotImplementedError


class ChebychovHelper(SpectralHelper1D):
class ChebychevHelper(SpectralHelper1D):
"""
The Chebychov base consists of special kinds of polynomials, with the main advantage that you can easily transform
The Chebychev base consists of special kinds of polynomials, with the main advantage that you can easily transform
between physical and spectral space by discrete cosine transform.
The differentiation in the Chebychov T base is dense, but can be preconditioned to yield a differentiation operator
that moves to Chebychov U basis during differentiation, which is sparse. When using this technique, problems need to
The differentiation in the Chebychev T base is dense, but can be preconditioned to yield a differentiation operator
that moves to Chebychev U basis during differentiation, which is sparse. When using this technique, problems need to
be formulated in first order formulation.
This implementation is largely based on the Dedalus paper (arXiv:1905.10388).
Expand All @@ -142,9 +142,9 @@ def __init__(self, *args, transform_type='fft', x0=-1, x1=1, **kwargs):

def get_1dgrid(self):
'''
Generates a 1D grid with Chebychov points. These are clustered at the boundary. You need this kind of grid to
use discrete cosine transformation (DCT) to get the Chebychov representation. If you want a different grid, you
need to do an affine transformation before any Chebychov business.
Generates a 1D grid with Chebychev points. These are clustered at the boundary. You need this kind of grid to
use discrete cosine transformation (DCT) to get the Chebychev representation. If you want a different grid, you
need to do an affine transformation before any Chebychev business.
Returns:
numpy.ndarray: 1D grid
Expand All @@ -157,8 +157,8 @@ def get_wavenumbers(self):
def get_conv(self, name, N=None):
'''
Get conversion matrix between different kinds of polynomials. The supported kinds are
- T: Chebychov polynomials of first kind
- U: Chebychov polynomials of second kind
- T: Chebychev polynomials of first kind
- U: Chebychev polynomials of second kind
- D: Dirichlet recombination.
You get the desired matrix by choosing a name as ``A2B``. I.e. ``T2U`` for the conversion matrix from T to U.
Expand Down Expand Up @@ -420,11 +420,11 @@ def get_Dirichlet_recombination_matrix(self):
return sp.eye(N) - sp.diags(xp.ones(N - 2), offsets=+2)


class Ultraspherical(ChebychovHelper):
class Ultraspherical(ChebychevHelper):
"""
This implementation follows https://doi.org/10.1137/120865458.
The ultraspherical method works in Chebychov polynomials as well, but also uses various Gegenbauer polynomials.
The idea is that for every derivative of Chebychov T polynomials, there is a basis of Gegenbauer polynomials where the differentiation matrix is a single off-diagonal.
The ultraspherical method works in Chebychev polynomials as well, but also uses various Gegenbauer polynomials.
The idea is that for every derivative of Chebychev T polynomials, there is a basis of Gegenbauer polynomials where the differentiation matrix is a single off-diagonal.
There are also conversion operators from one derivative basis to the next that are sparse.
This basis is used like this: For every equation that you have, look for the highest derivative and bump all matrices to the correct basis. If your highest derivative is 2 and you have an identity, it needs to get bumped from 0 to 1 and from 1 to 2. If you have a first derivative as well, it needs to be bumped from 1 to 2.
Expand Down Expand Up @@ -672,7 +672,7 @@ def add_axis(self, base, *args, **kwargs):

if base.lower() in ['chebychov', 'chebychev', 'cheby', 'chebychovhelper']:
kwargs['transform_type'] = kwargs.get('transform_type', 'fft')
self.axes.append(ChebychovHelper(*args, **kwargs))
self.axes.append(ChebychevHelper(*args, **kwargs))
elif base.lower() in ['fft', 'fourier', 'ffthelper']:
self.axes.append(FFTHelper(*args, **kwargs))
elif base.lower() in ['ultraspherical', 'gegenbauer']:
Expand Down Expand Up @@ -1146,7 +1146,7 @@ def transform_single_component(self, u, axes=None, padding=None):
Transform a single component of the solution
"""
trfs = {
ChebychovHelper: self._transform_dct,
ChebychevHelper: self._transform_dct,
Ultraspherical: self._transform_dct,
FFTHelper: self._transform_fft,
}
Expand Down Expand Up @@ -1283,7 +1283,7 @@ def _transform_idct(self, u, axes, padding=None, **kwargs):
def itransform_single_component(self, u, axes=None, padding=None):
trfs = {
FFTHelper: self._transform_ifft,
ChebychovHelper: self._transform_idct,
ChebychevHelper: self._transform_idct,
Ultraspherical: self._transform_idct,
}

Expand Down Expand Up @@ -1504,10 +1504,10 @@ def get_Id(self):

def get_Dirichlet_recombination_matrix(self, axis=-1, **kwargs):
"""
Get Dirichlet recombination matrix along axis. Not that it only makes sense in directions discretized with variations of Chebychov bases.
Get Dirichlet recombination matrix along axis. Not that it only makes sense in directions discretized with variations of Chebychev bases.
Args:
axis (int): Axis you discretized with Chebychov
axis (int): Axis you discretized with Chebychev
Returns:
sparse matrix
Expand Down
25 changes: 22 additions & 3 deletions pySDC/implementations/problem_classes/Burgers.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,17 @@

class Burgers1D(GenericSpectralLinear):
"""
See https://en.wikipedia.org/wiki/Burgers'_equation
Discretization is done with a Chebychov method, which requires a first order derivative formulation.
See https://en.wikipedia.org/wiki/Burgers'_equation for the equation that is solved.
Discretization is done with a Chebychev method, which requires a first order derivative formulation.
Feel free to do a more efficient implementation using an ultraspherical method to avoid the first order business.
Parameters:
N (int): Spatial resolution
epsilon (float): viscosity
BCl (float): Value at left boundary
BCr (float): Value at right boundary
f (int): Frequency of the initial conditions
mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
"""

dtype_u = mesh
Expand Down Expand Up @@ -131,7 +139,18 @@ def plot(self, u, t=None, fig=None, comp='u'): # pragma: no cover

class Burgers2D(GenericSpectralLinear):
"""
See documentation of `Burgers1D`.
See https://en.wikipedia.org/wiki/Burgers'_equation for the equation that is solved.
This implementation is discretized with FFTs in x and Chebychev in z.
Parameters:
nx (int): Spatial resolution in x direction
nz (int): Spatial resolution in z direction
epsilon (float): viscosity
BCl (float): Value at left boundary
BCr (float): Value at right boundary
fux (int): Frequency of the initial conditions in x-direction
fuz (int): Frequency of the initial conditions in z-direction
mode (str): 'T2U' or 'T2T'. Use 'T2U' to get sparse differentiation matrices
"""

dtype_u = mesh
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,14 +6,14 @@
from pySDC.implementations.problem_classes.generic_spectral import GenericSpectralLinear


class Heat1DChebychov(GenericSpectralLinear):
class Heat1DChebychev(GenericSpectralLinear):
dtype_u = mesh
dtype_f = mesh

def __init__(self, nvars=128, a=0, b=0, f=1, nu=1.0, mode='T2U', **kwargs):
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)

bases = [{'base': 'chebychov', 'N': nvars}]
bases = [{'base': 'chebychev', 'N': nvars}]
components = ['u', 'ux']

super().__init__(bases, components, **kwargs)
Expand Down Expand Up @@ -153,11 +153,11 @@ def u_exact(self, t, noise=0):
return u


class Heat2DChebychov(GenericSpectralLinear):
class Heat2DChebychev(GenericSpectralLinear):
dtype_u = mesh
dtype_f = mesh

def __init__(self, nx=128, ny=128, base_x='fft', base_y='chebychov', a=0, b=0, c=0, fx=1, fy=1, nu=1.0, **kwargs):
def __init__(self, nx=128, ny=128, base_x='fft', base_y='chebychev', a=0, b=0, c=0, fx=1, fy=1, nu=1.0, **kwargs):
assert nx % 2 == 1 or base_x == 'fft'
assert ny % 2 == 1 or base_y == 'fft'
self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
Expand All @@ -184,20 +184,20 @@ def __init__(self, nx=128, ny=128, base_x='fft', base_y='chebychov', a=0, b=0, c
self.setup_M(M_lhs)

for base in [base_x, base_y]:
assert base in ['chebychov', 'fft']
assert base in ['chebychev', 'fft']

alpha = (self.b - self.a) / 2.0
beta = (self.c - self.b) / 2.0
gamma = (self.c + self.a) / 2.0

if base_x == 'chebychov':
if base_x == 'chebychev':
y = self.Y[0, :]
if self.base_y == 'fft':
self.add_BC(component='u', equation='u', axis=0, x=-1, v=beta * y - alpha + gamma, kind='Dirichlet')
self.add_BC(component='ux', equation='ux', axis=0, v=2 * alpha, kind='integral')
else:
assert a == b, f'Need periodic boundary conditions in x for {base_x} method!'
if base_y == 'chebychov':
if base_y == 'chebychev':
x = self.X[:, 0]
self.add_BC(component='u', equation='u', axis=1, x=-1, v=alpha * x - beta + gamma, kind='Dirichlet')
self.add_BC(component='uy', equation='uy', axis=1, v=2 * beta, kind='integral')
Expand Down
15 changes: 13 additions & 2 deletions pySDC/implementations/problem_classes/RayleighBenard.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,19 @@


class RayleighBenard(GenericSpectralLinear):
"""
Rayleigh-Benard Convection is a variation of incompressible Navier-Stokes. See, for instance https://doi.org/10.1007/s00791-020-00332-3.
Parameters:
Prandl (float): Prandl number
Rayleigh (float): Rayleigh number
nx (int): Horizontal resolution
nz (int): Vertical resolution
BCs (dict): Can specify boundary conditions here
dealiasing (float): Dealiasing factor for evaluating the non-linear part
comm (mpi4py.Intracomm): Space communicator
"""

dtype_u = mesh
dtype_f = imex_mesh

Expand All @@ -21,7 +34,6 @@ def __init__(
BCs=None,
dealiasing=3 / 2,
comm=None,
debug=False,
**kwargs,
):
BCs = {} if BCs is None else BCs
Expand Down Expand Up @@ -50,7 +62,6 @@ def __init__(
'BCs',
'dealiasing',
'comm',
'debug',
localVars=locals(),
readOnly=True,
)
Expand Down
2 changes: 1 addition & 1 deletion pySDC/implementations/problem_classes/generic_spectral.py
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ def setup_preconditioner(self, Dirichlet_recombination=True, left_preconditioner

self.Pl = self.spectral.sparse_lib.csc_matrix(R)

if Dirichlet_recombination and type(self.axes[-1]).__name__ in ['ChebychovHelper, Ultraspherical']:
if Dirichlet_recombination and type(self.axes[-1]).__name__ in ['ChebychevHelper, Ultraspherical']:
_Pr = self.spectral.get_Dirichlet_recombination_matrix(axis=-1)
else:
_Pr = Id
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,9 @@
@pytest.mark.parametrize('N', [4, 32])
def test_D2T_conversion_matrices(N):
import numpy as np
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)

x = np.linspace(-1, 1, N)
D2T = cheby.get_conv('D2T')
Expand All @@ -31,9 +31,9 @@ def test_D2T_conversion_matrices(N):
def test_T_U_conversion(N):
import numpy as np
from scipy.special import chebyt, chebyu
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)

T2U = cheby.get_conv('T2U')
U2T = cheby.get_conv('U2T')
Expand Down Expand Up @@ -61,11 +61,11 @@ def eval_poly(poly, coeffs, x):
@pytest.mark.base
@pytest.mark.parametrize('name', ['T2U', 'T2D', 'T2T'])
def test_conversion_inverses(name):
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper
import numpy as np

N = 8
cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
P = cheby.get_conv(name)
Pinv = cheby.get_conv(name[::-1])
assert np.allclose((P @ Pinv).toarray(), np.diag(np.ones(N)))
Expand All @@ -76,9 +76,9 @@ def test_conversion_inverses(name):
def test_differentiation_matrix(N):
import numpy as np
import scipy
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
x = np.cos(np.pi / N * (np.arange(N) + 0.5))
coeffs = np.random.random(N)
norm = cheby.get_norm()
Expand All @@ -95,9 +95,9 @@ def test_differentiation_matrix(N):
@pytest.mark.parametrize('N', [4, 32])
def test_integration_matrix(N):
import numpy as np
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
coeffs = np.random.random(N)
coeffs[-1] = 0

Expand All @@ -116,9 +116,9 @@ def test_integration_matrix(N):
def test_transform(N, d, transform_type):
import scipy
import numpy as np
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper

cheby = ChebychovHelper(N, transform_type=transform_type)
cheby = ChebychevHelper(N, transform_type=transform_type)
u = np.random.random((d, N))
norm = cheby.get_norm()
x = cheby.get_1dgrid()
Expand All @@ -137,10 +137,10 @@ def test_transform(N, d, transform_type):
@pytest.mark.base
@pytest.mark.parametrize('N', [8, 32])
def test_integration_BC(N):
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper
import numpy as np

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
coef = np.random.random(N)

BC = cheby.get_integ_BC_row()
Expand All @@ -155,11 +155,11 @@ def test_integration_BC(N):
@pytest.mark.base
@pytest.mark.parametrize('N', [4, 32])
def test_norm(N):
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper
import numpy as np
import scipy

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
coeffs = np.random.random(N)
x = cheby.get_1dgrid()
norm = cheby.get_norm()
Expand All @@ -185,10 +185,10 @@ def test_tau_method(bc, N, bc_val):
The test verifies that the solution satisfies the perturbed equation and the boundary condition.
'''
from pySDC.helpers.spectral_helper import ChebychovHelper
from pySDC.helpers.spectral_helper import ChebychevHelper
import numpy as np

cheby = ChebychovHelper(N)
cheby = ChebychevHelper(N)
x = cheby.get_1dgrid()

coef = np.append(np.zeros(N - 1), [1])
Expand Down Expand Up @@ -221,13 +221,13 @@ def test_tau_method(bc, N, bc_val):
def test_tau_method2D(bc, nz, nx, bc_val, plotting=False):
'''
solve u_z - 0.1u_xx -u_x + tau P = 0, u(bc) = sin(bc_val*x) -> space-time discretization of advection-diffusion
problem. We do FFT in x-direction and Chebychov in z-direction.
problem. We do FFT in x-direction and Chebychev in z-direction.
'''
from pySDC.helpers.spectral_helper import ChebychovHelper, FFTHelper
from pySDC.helpers.spectral_helper import ChebychevHelper, FFTHelper
import numpy as np
import scipy.sparse as sp

cheby = ChebychovHelper(nz)
cheby = ChebychevHelper(nz)
fft = FFTHelper(nx)

# generate grid
Expand All @@ -239,7 +239,7 @@ def test_tau_method2D(bc, nz, nx, bc_val, plotting=False):
bcs = np.sin(bc_val * x)
rhs = np.zeros_like(X)
rhs[:, -1] = bcs
rhs_hat = fft.transform(rhs, axis=-2) # the rhs is already in Chebychov spectral space
rhs_hat = fft.transform(rhs, axis=-2) # the rhs is already in Chebychev spectral space

# generate matrices
Dx = fft.get_differentiation_matrix(p=2) * 1e-1 + fft.get_differentiation_matrix()
Expand Down Expand Up @@ -299,11 +299,11 @@ def test_tau_method2D_diffusion(nz, nx, bc_val, plotting=False):
'''
Solve a Poisson problem with funny Dirichlet BCs in z-direction and periodic in x-direction.
'''
from pySDC.helpers.spectral_helper import ChebychovHelper, FFTHelper
from pySDC.helpers.spectral_helper import ChebychevHelper, FFTHelper
import numpy as np
import scipy.sparse as sp

cheby = ChebychovHelper(nz)
cheby = ChebychevHelper(nz)
fft = FFTHelper(nx)

# generate grid
Expand All @@ -315,7 +315,7 @@ def test_tau_method2D_diffusion(nz, nx, bc_val, plotting=False):
rhs = np.zeros((2, nx, nz)) # components u and u_x
rhs[0, :, -1] = np.sin(bc_val * x) + 1
rhs[1, :, -1] = 3 * np.exp(-((x - 3.6) ** 2)) + np.cos(x)
rhs_hat = fft.transform(rhs, axis=-2) # the rhs is already in Chebychov spectral space
rhs_hat = fft.transform(rhs, axis=-2) # the rhs is already in Chebychev spectral space

# generate 1D matrices
Dx = fft.get_differentiation_matrix()
Expand Down
Loading

0 comments on commit 6fbfa72

Please sign in to comment.