diff --git a/pySDC/helpers/spectral_helper.py b/pySDC/helpers/spectral_helper.py
index 49079d0791453ffc5650d246b366833d62966777..2cf912fbf4abd59d1946a7e7924c5bc8336c85a9 100644
--- a/pySDC/helpers/spectral_helper.py
+++ b/pySDC/helpers/spectral_helper.py
@@ -1,11 +1,24 @@
import numpy as np
import scipy
from pySDC.implementations.datatype_classes.mesh import mesh
+from scipy.special import factorial
class SpectralHelper1D:
"""
- Base class for 1D spectral discretizations.
+ Abstract base class for 1D spectral discretizations. Defines a common interface with parameters and functions that
+ all bases need to have.
+
+ When implementing new bases, please take care to use the modules that are supplied as class attributes to enable
+ the code for GPUs.
+
+ Attributes:
+ N (int): Resolution
+ x0 (float): Coordinate of left boundary
+ x1 (float): Coordinate of right boundary
+ L (float): Length of the domain
+ useGPU (bool): Whether to use GPUs
+
"""
fft_lib = scipy.fft
@@ -13,7 +26,16 @@ class SpectralHelper1D:
linalg = scipy.sparse.linalg
xp = np
- def __init__(self, N, x0=None, x1=None, useGPU=False, **kwargs):
+ def __init__(self, N, x0=None, x1=None, useGPU=False):
+ """
+ Constructor
+
+ Args:
+ N (int): Resolution
+ x0 (float): Coordinate of left boundary
+ x1 (float): Coordinate of right boundary
+ useGPU (bool): Whether to use GPUs
+ """
self.N = N
self.x0 = x0
self.x1 = x1
@@ -38,9 +60,21 @@ class SpectralHelper1D:
cls.fft_lib = fft_lib
def get_Id(self):
+ """
+ Get identity matrix
+
+ Returns:
+ sparse diagonal identity matrix
+ """
return self.sparse_lib.eye(self.N)
def get_zero(self):
+ """
+ Get a matrix with all zeros of the correct size.
+
+ Returns:
+ sparse matrix with zeros everywhere
+ """
return 0 * self.get_Id()
def get_differentiation_matrix(self):
@@ -70,18 +104,31 @@ class SpectralHelper1D:
def get_basis_change_matrix(self, *args, **kwargs):
"""
- Some spectral discretization change the basis during differentiation. This method can be used to transfer between the various bases.
+ Some spectral discretization change the basis during differentiation. This method can be used to transfer
+ between the various bases.
+
+ This method accepts arbitrary arguments that may not be used in order to provide an easy interface for multi-
+ dimensional bases. For instance, you may combine an FFT discretization with an ultraspherical discretization.
+ The FFT discretization will always be in the same base, but the ultraspherical discretization uses a different
+ base for every derivative. You can then ask all bases for transfer matrices from one ultraspherical derivative
+ base to the next. The FFT discretization will ignore this and return an identity while the ultraspherical
+ discretization will return the desired matrix. After a Kronecker product, you get the 2D version of the matrix
+ you want. This is what the `SpectralHelper` does when you call the method of the same name on it.
+
+ Returns:
+ sparse bases change matrix
"""
return self.sparse_lib.eye(self.N)
- def get_BC(self, kind, **kwargs):
+ def get_BC(self, kind):
"""
To facilitate boundary conditions (BCs) we use either a basis where all functions satisfy the BCs automatically,
as is the case in FFT basis for periodic BCs, or boundary bordering. In boundary bordering, specific lines in
the matrix are replaced by the boundary conditions as obtained by this method.
Args:
- kind (str): The type of BC you want to implement
+ kind (str): The type of BC you want to implement please refer to the implementations of this method in the
+ individual 1D bases for what is implemented
Returns:
self.xp.array: Boundary condition
@@ -136,6 +183,17 @@ class ChebychevHelper(SpectralHelper1D):
"""
def __init__(self, *args, transform_type='fft', x0=-1, x1=1, **kwargs):
+ """
+ Constructor.
+ Please refer to the parent class for additional arguments. Notably, you have to supply a resolution `N` and you
+ may choose to run on GPUs via the `useGPU` argument.
+
+ Args:
+ transform_type ('fft' or 'dct'): Either use DCT functions directly implemented in the transform library or
+ use the FFT from the library to compute the DCT
+ 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
super().__init__(*args, x0=x0, x1=x1, **kwargs)
@@ -159,6 +217,7 @@ class ChebychevHelper(SpectralHelper1D):
return self.xp.cos(np.pi / self.N * (self.xp.arange(self.N) + 0.5))
def get_wavenumbers(self):
+ """Get the domain in spectral space"""
return self.xp.arange(self.N)
def get_conv(self, name, N=None):
@@ -215,9 +274,30 @@ class ChebychevHelper(SpectralHelper1D):
return mat
def get_basis_change_matrix(self, conv='T2T', **kwargs):
+ """
+ As the differentiation matrix in Chebychev-T base is dense but is sparse when simultaneously changing base to
+ Chebychev-U, you may need a basis change matrix to transfer the other matrices as well. This function returns a
+ conversion matrix from `ChebychevHelper.get_conv`. Not that `**kwargs` are used to absorb arguments for other
+ bases, see documentation of `SpectralHelper1D.get_basis_change_matrix`.
+
+ Args:
+ conv (str): Conversion code, i.e. T2U
+
+ Returns:
+ Sparse conversion matrix
+ """
return self.get_conv(conv)
def get_integration_matrix(self, lbnd=0):
+ """
+ Get matrix for integration
+
+ Args:
+ lbnd (float): Lower bound for integration, only 0 is currently implemented
+
+ Returns:
+ Sparse integration matrix
+ """
S = self.sparse_lib.diags(1 / (self.xp.arange(self.N - 1) + 1), offsets=-1) @ self.get_conv('T2U')
n = self.xp.arange(self.N)
if lbnd == 0:
@@ -228,7 +308,7 @@ class ChebychevHelper(SpectralHelper1D):
/ (np.append([1], self.xp.arange(self.N // 2 - 1) + 1))
)
else:
- raise NotImplementedError
+ raise NotImplementedError(f'This function allows to integrate only from x=0, you attempted from x={lbnd}.')
return S
def get_differentiation_matrix(self, p=1):
@@ -250,7 +330,15 @@ class ChebychevHelper(SpectralHelper1D):
return self.sparse_lib.csc_matrix(self.xp.linalg.matrix_power(D, p))
def get_norm(self, N=None):
- '''get normalization for converting Chebychev coefficients and DCT'''
+ '''
+ Get normalization for converting Chebychev coefficients and DCT
+
+ Args:
+ N (int, optional): Resolution
+
+ Returns:
+ self.xp.array: Normalization
+ '''
N = self.N if N is None else N
norm = self.xp.ones(N) / N
norm[0] /= 2
@@ -321,9 +409,19 @@ class ChebychevHelper(SpectralHelper1D):
return self.fft_utils
def transform(self, u, axis=-1, **kwargs):
- if self.transform_type == 'dct':
- return self.fft_lib.dct(u, axis=axis) * self.norm
- elif self.transform_type == 'fft':
+ """
+ 1D DCT along axis. `kwargs` will be passed on to the FFT library.
+
+ Args:
+ u: Data you want to transform
+ axis (int): Axis you want to transform along
+
+ Returns:
+ Data in spectral space
+ """
+ if self.transform_type.lower() == 'dct':
+ return self.fft_lib.dct(u, axis=axis, **kwargs) * self.norm
+ elif self.transform_type.lower() == 'fft':
result = u.copy()
shuffle = [slice(0, s, 1) for s in u.shape]
@@ -341,9 +439,19 @@ class ChebychevHelper(SpectralHelper1D):
result.real[...] = V.real[...]
return result
else:
- raise NotImplementedError
+ raise NotImplementedError(f'Please choose a transform type from fft and dct, not {self.transform_type=}')
def itransform(self, u, axis=-1):
+ """
+ 1D inverse DCT along axis.
+
+ Args:
+ u: Data you want to transform
+ axis (int): Axis you want to transform along
+
+ Returns:
+ Data in physical space
+ """
assert self.norm.shape[0] == u.shape[axis]
if self.transform_type == 'dct':
@@ -367,16 +475,21 @@ class ChebychevHelper(SpectralHelper1D):
def get_BC(self, kind, **kwargs):
"""
- Get boundary condition row for boundary bordering.
+ Get boundary condition row for boundary bordering. `kwargs` will be passed on to implementations of the BC of
+ the kind you choose. Specifically, `x` for `'dirichlet'` boundary condition, which is the coordinate at which to
+ set the BC.
+
+ Args:
+ kind ('integral' or 'dirichlet'): Kind of boundary condition you want
"""
if kind.lower() == 'integral':
return self.get_integ_BC_row(**kwargs)
elif kind.lower() == 'dirichlet':
return self.get_Dirichlet_BC_row(**kwargs)
else:
- return super().get_BC(kind, **kwargs)
+ return super().get_BC(kind)
- def get_integ_BC_row(self, **kwargs):
+ def get_integ_BC_row(self):
"""
Get a row for generating integral BCs with T polynomials.
It returns the values of the integrals of T polynomials over the entire interval.
@@ -390,7 +503,7 @@ class ChebychevHelper(SpectralHelper1D):
me[0] = 2.0
return me
- def get_Dirichlet_BC_row(self, x, **kwargs):
+ def get_Dirichlet_BC_row(self, x):
"""
Get a row for generating Dirichlet BCs at x with T polynomials.
It returns the values of the T polynomials at x.
@@ -427,7 +540,7 @@ class ChebychevHelper(SpectralHelper1D):
return sp.eye(N) - sp.diags(xp.ones(N - 2), offsets=+2)
-class Ultraspherical(ChebychevHelper):
+class UltrasphericalHelper(ChebychevHelper):
"""
This implementation follows https://doi.org/10.1137/120865458.
The ultraspherical method works in Chebychev polynomials as well, but also uses various Gegenbauer polynomials.
@@ -452,13 +565,12 @@ class Ultraspherical(ChebychevHelper):
xp = self.xp
N = self.N
l = p
- from scipy.special import factorial
-
return 2 ** (l - 1) * factorial(l - 1) * sp.diags(xp.arange(N - l) + l, offsets=l)
def get_S(self, lmbda):
"""
- Get matrix for bumping the derivative base by one
+ Get matrix for bumping the derivative base by one from lmbda to lmbda + 1. This is the same language as in
+ https://doi.org/10.1137/120865458.
Args:
lmbda (int): Ingoing derivative base
@@ -481,7 +593,7 @@ class Ultraspherical(ChebychevHelper):
return self.sparse_lib.csc_matrix(mat)
- def get_conv(self, p_out, p_in=0):
+ def get_basis_change_matrix(self, p_in=0, p_out=0, **kwargs):
"""
Get a conversion matrix from derivative base `p_in` to `p_out`.
@@ -497,6 +609,7 @@ class Ultraspherical(ChebychevHelper):
return mat_fwd
else:
+ # We have to invert the matrix on CPU because the GPU equivalent is not implemented in CuPy at the time of writing.
import scipy.sparse as sp
if self.useGPU:
@@ -506,13 +619,49 @@ class Ultraspherical(ChebychevHelper):
return self.sparse_lib.csc_matrix(mat_bck)
- def get_basis_change_matrix(self, p_in=0, p_out=0, **kwargs):
- return self.get_conv(p_in=p_in, p_out=p_out)
-
def get_integration_matrix(self):
- return self.sparse_lib.diags(1 / (self.xp.arange(self.N - 1) + 1), offsets=-1) @ self.get_conv(p_out=1, p_in=0)
+ """
+ Get an integration matrix. Please use `UltrasphericalHelper.get_integration_constant` afterwards to compute the
+ integration constant such that integration starts from x=-1.
+
+ Example:
+
+ .. code-block:: python
+
+ import numpy as np
+ from pySDC.helpers.spectral_helper import UltrasphericalHelper
+
+ N = 4
+ helper = UltrasphericalHelper(N)
+ coeffs = np.random.random(N)
+ coeffs[-1] = 0
+
+ poly = np.polynomial.Chebyshev(coeffs)
+
+ S = helper.get_integration_matrix()
+ U_hat = S @ coeffs
+ U_hat[0] = helper.get_integration_constant(U_hat, axis=-1)
+
+ assert np.allclose(poly.integ(lbnd=-1).coef[:-1], U_hat)
+
+ 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
+ )
def get_integration_constant(self, u_hat, axis):
+ """
+ Get integration constant for lower bound of -1. See documentation of `UltrasphericalHelper.get_integration_matrix` for details.
+
+ Args:
+ u_hat: Solution in spectral space
+ axis: Axis you want to integrate over
+
+ Returns:
+ Integration constant, has one less dimension than `u_hat`
+ """
slices = [
None,
] * u_hat.ndim
@@ -522,6 +671,17 @@ class Ultraspherical(ChebychevHelper):
class FFTHelper(SpectralHelper1D):
def __init__(self, *args, x0=0, x1=2 * np.pi, **kwargs):
+ """
+ Constructor.
+ Please refer to the parent class for additional arguments. Notably, you have to supply a resolution `N` and you
+ may choose to run on GPUs via the `useGPU` argument.
+
+ Args:
+ transform_type ('fft' or 'dct'): Either use DCT functions directly implemented in the transform library or
+ use the FFT from the library to compute the DCT
+ x0 (float, optional): Coordinate of left boundary
+ x1 (float, optional): Coordinate of right boundary
+ """
super().__init__(*args, x0=x0, x1=x1, **kwargs)
def get_1dgrid(self):
@@ -548,28 +708,70 @@ class FFTHelper(SpectralHelper1D):
sparse differentiation matrix
"""
k = self.get_wavenumbers()
- import scipy.sparse as sp
if self.useGPU:
+ # Have to raise the matrix to power p on CPU because the GPU equivalent is not implemented in CuPy at the time of writing.
+ import scipy.sparse as sp
+
D = self.sparse_lib.diags(1j * k).get()
return self.sparse_lib.csc_matrix(sp.linalg.matrix_power(D, p))
- elif self.sparse_lib == sp:
+ else:
return self.linalg.matrix_power(self.sparse_lib.diags(1j * k), p)
def get_integration_matrix(self, p=1):
+ """
+ Get integration matrix to compute `p`-th integral over the entire domain.
+
+ Args:
+ p (int): Order of integral you want to compute
+
+ Returns:
+ sparse integration matrix
+ """
k = self.xp.array(self.get_wavenumbers(), dtype='complex128')
k[0] = 1j * self.L
return self.linalg.matrix_power(self.sparse_lib.diags(1 / (1j * k)), p)
def transform(self, u, axis=-1, **kwargs):
+ """
+ 1D FFT along axis. `kwargs` are passed on to the FFT library.
+
+ Args:
+ u: Data you want to transform
+ axis (int): Axis you want to transform along
+
+ Returns:
+ transformed data
+ """
return self.fft_lib.fft(u, axis=axis, **kwargs)
def itransform(self, u, axis=-1):
+ """
+ Inverse 1D FFT.
+
+ Args:
+ u: Data you want to transform
+ axis (int): Axis you want to transform along
+
+ Returns:
+ transformed data
+ """
return self.fft_lib.ifft(u, axis=axis)
- def get_BC(self, kind, **kwargs):
+ def get_BC(self, kind):
+ """
+ Get a sort of boundary condition. You can use `kind=integral`, to fix the integral, or you can use `kind=Nyquist`.
+ The latter is not really a boundary condition, but is used to set the Nyquist mode to some value, preferably zero.
+ You should set the Nyquist mode zero when the solution in physical space is real and the resolution is even.
+
+ Args:
+ kind ('integral' or 'nyquist'): Kind of BC
+
+ Returns:
+ self.xp.ndarray: Boundary condition row
+ """
if kind.lower() == 'integral':
- return self.get_integ_BC_row(**kwargs)
+ return self.get_integ_BC_row()
elif kind.lower() == 'nyquist':
assert (
self.N % 2 == 0
@@ -578,14 +780,23 @@ class FFTHelper(SpectralHelper1D):
BC[self.get_Nyquist_mode_index()] = 1
return BC
else:
- return super().get_BC(kind, **kwargs)
+ return super().get_BC(kind)
def get_Nyquist_mode_index(self):
+ """
+ Compute the index of the Nyquist mode, i.e. the mode with the lowest wavenumber, which doesn't have a positive
+ counterpart for even resolution. This means real waves of this wave number cannot be properly resolved and you
+ are best advised to set this mode zero if representing real functions on even-resolution grids is what you're
+ after.
+
+ Returns:
+ int: Index of the Nyquist mode
+ """
k = self.get_wavenumbers()
Nyquist_mode = min(k)
return self.xp.where(k == Nyquist_mode)[0][0]
- def get_integ_BC_row(self, **kwargs):
+ def get_integ_BC_row(self):
"""
Only the 0-mode has non-zero integral with FFT basis in periodic BCs
"""
@@ -600,6 +811,26 @@ class SpectralHelper:
- Easily assemble matrices containing multiple equations
- Direct product of 1D bases to solve problems in more dimensions
- Distribute the FFTs to facilitate concurrency.
+
+ Attributes:
+ comm (mpi4py.Intracomm): MPI communicator
+ debug (bool): Perform additional checks at extra computational cost
+ useGPU (bool): Whether to use GPUs
+ axes (list): List of 1D bases
+ components (list): List of strings of the names of components in the equations
+ full_BCs (list): List of Dictionaries containing all information about the boundary conditions
+ BC_mat (list): List of lists of sparse matrices to put BCs into and eventually assemble the BC matrix from
+ BCs (sparse matrix): Matrix containing only the BCs
+ fft_cache (dict): Cache FFTs of various shapes here to facilitate padding and so on
+ BC_rhs_mask (self.xp.ndarray): Mask values that contain boundary conditions in the right hand side
+ BC_zero_index (self.xp.ndarray): Indeces of rows in the matrix that are replaced by BCs
+ BC_line_zero_matrix (sparse matrix): Matrix that zeros rows where we can then add the BCs in using `BCs`
+ rhs_BCs_hat (self.xp.ndarray): Boundary conditions in spectral space
+ global_shape (tuple): Global shape of the solution as in `mpi4py-fft`
+ local_slice (slice): Local slice of the solution as in `mpi4py-fft`
+ fft_obj: When using distributed FFTs, this will be a parallel transform object from `mpi4py-fft`
+ init (tuple): This is the same `init` that is used throughout the problem classes
+ init_forward (tuple): This is the equivalent of `init` in spectral space
"""
xp = np
@@ -628,6 +859,14 @@ class SpectralHelper:
cls.dtype = cupy_mesh
def __init__(self, comm=None, useGPU=False, debug=False):
+ """
+ Constructor
+
+ Args:
+ comm (mpi4py.Intracomm): MPI communicator
+ useGPU (bool): Whether to use GPUs
+ debug (bool): Perform additional checks at extra computational cost
+ """
self.comm = comm
self.debug = debug
self.useGPU = useGPU
@@ -683,6 +922,8 @@ class SpectralHelper:
def add_axis(self, base, *args, **kwargs):
"""
Add an axis to the domain by deciding on suitable 1D base.
+ Arguments to the bases are forwarded using `*args` and `**kwargs`. Please refer to the documentation of the 1D
+ bases for possible arguments.
Args:
base (str): 1D spectral method
@@ -695,7 +936,7 @@ class SpectralHelper:
elif base.lower() in ['fft', 'fourier', 'ffthelper']:
self.axes.append(FFTHelper(*args, **kwargs))
elif base.lower() in ['ultraspherical', 'gegenbauer']:
- self.axes.append(Ultraspherical(*args, **kwargs))
+ self.axes.append(UltrasphericalHelper(*args, **kwargs))
else:
raise NotImplementedError(f'{base=!r} is not implemented!')
self.axes[-1].xp = self.xp
@@ -751,6 +992,7 @@ class SpectralHelper:
Use this method for boundary bordering. It gets the respective matrix row and embeds it into a matrix.
Pay attention that if you have multiple BCs in a single equation, you need to put them in different lines.
Typically, the last line that does not contain a BC is the best choice.
+ Forward arguments for the boundary conditions using `kwargs`. Refer to documentation of 1D bases for details.
Args:
axis (int): Axis you want to add the BC to
@@ -798,6 +1040,7 @@ class SpectralHelper:
"""
Remove a BC from the matrix. This is useful e.g. when you add a non-scalar BC and then need to selectively
remove single BCs again, as in incompressible Navier-Stokes, for instance.
+ Forward arguments for the boundary conditions using `kwargs`. Refer to documentation of 1D bases for details.
Args:
component (str): Name of the component the BC should act on
@@ -831,6 +1074,7 @@ class SpectralHelper:
"""
Add a BC to the matrix. Note that you need to convert the list of lists of BCs that this method generates to a
single sparse matrix by calling `setup_BCs` after adding/removing all BCs.
+ Forward arguments for the boundary conditions using `kwargs`. Refer to documentation of 1D bases for details.
Args:
component (str): Name of the component the BC should act on
@@ -910,10 +1154,16 @@ class SpectralHelper:
if len(BCs) > 0:
u_hat = self.transform(u, axes=(axis - self.ndim,))
for BC in BCs:
+ kwargs = {
+ key: value
+ for key, value in BC.items()
+ if key not in ['component', 'equation', 'axis', 'v', 'line', 'scalar']
+ }
+
if axis == 0:
- get = self.axes[axis].get_BC(**BC) @ u_hat[self.index(BC['component'])]
+ get = self.axes[axis].get_BC(**kwargs) @ u_hat[self.index(BC['component'])]
elif axis == 1:
- get = u_hat[self.index(BC['component'])] @ self.axes[axis].get_BC(**BC)
+ get = u_hat[self.index(BC['component'])] @ self.axes[axis].get_BC(**kwargs)
want = BC['v']
assert self.xp.allclose(
get, want
@@ -996,7 +1246,32 @@ class SpectralHelper:
"""
Add the left hand part (that you want to solve implicitly) of an equation to a list of lists of sparse matrices
that you will convert to an operator later.
- For instance, the `relations` for diffusion in component `u` would look like `relations = {'u': -Dxx}`.
+
+ Example:
+ Setup linear operator `L` for 1D heat equation using Chebychev method in first order form and T-to-U
+ preconditioning:
+
+ .. code-block:: python
+ helper = SpectralHelper()
+
+ helper.add_axis(base='chebychev', N=8)
+ helper.add_component(['u', 'ux'])
+ helper.setup_fft()
+
+ I = helper.get_Id()
+ Dx = helper.get_differentiation_matrix(axes=(0,))
+ T2U = helper.get_basis_change_matrix('T2U')
+
+ L_lhs = {
+ 'ux': {'u': -T2U @ Dx, 'ux': T2U @ I},
+ 'u': {'ux': -(T2U @ Dx)},
+ }
+
+ operator = helper.get_empty_operator_matrix()
+ for line, equation in L_lhs.items():
+ helper.add_equation_lhs(operator, line, equation)
+
+ L = helper.convert_operator_matrix_to_operator(operator)
Args:
A (list of lists of sparse matrices): The operator to be
@@ -1009,6 +1284,7 @@ class SpectralHelper:
def convert_operator_matrix_to_operator(self, M):
"""
Promote the list of lists of sparse matrices to a single sparse matrix that can be used as linear operator.
+ See documentation of `SpectralHelper.add_equation_lhs` for an example.
Args:
M (list of lists of sparse matrices): The operator to be
@@ -1035,7 +1311,7 @@ class SpectralHelper:
grids = [self.axes[i].get_1dgrid()[self.local_slice[i]] for i in range(len(self.axes))][::-1]
return self.xp.meshgrid(*grids)
- def get_fft(self, axes=None, direction='object', padding=None, shape=None, comm=None):
+ def get_fft(self, axes=None, direction='object', padding=None, shape=None):
"""
When using MPI, we use `PFFT` objects generated by mpi4py-fft
@@ -1044,7 +1320,6 @@ class SpectralHelper:
direction (str): use "forward" or "backward" to get functions for performing the transforms or "object" to get the PFFT object
padding (tuple): Padding for dealiasing
shape (tuple): Shape of the transform
- comm (mpi4py.Intracomm): MPI communicator
Returns:
transform
@@ -1075,7 +1350,7 @@ class SpectralHelper:
from mpi4py_fft import PFFT
_fft = PFFT(
- comm=self.comm if comm is None else comm,
+ comm=self.comm,
shape=shape,
axes=sorted(axes),
dtype='D',
@@ -1134,10 +1409,6 @@ class SpectralHelper:
)
self.BC_mat = self.get_empty_operator_matrix()
- self.BC_mask = self.xp.zeros(
- shape=self.init[0],
- dtype=bool,
- )
self.BC_rhs_mask = self.xp.zeros(
shape=self.init[0],
dtype=bool,
@@ -1154,6 +1425,7 @@ class SpectralHelper:
Returns:
transformed solution
"""
+ # TODO: clean up and try putting more of this in the 1D bases
fft = self.get_fft(axes, 'forward', **kwargs)
return fft(u, axes=axes)
@@ -1171,6 +1443,7 @@ class SpectralHelper:
Returns:
transformed solution
'''
+ # TODO: clean up and try putting more of this in the 1D bases
if self.debug:
assert self.xp.allclose(u.imag, 0), 'This function can only handle real input.'
@@ -1216,10 +1489,19 @@ class SpectralHelper:
def transform_single_component(self, u, axes=None, padding=None):
"""
Transform a single component of the solution
+
+ Args:
+ u data to transform:
+ axes (tuple): Axes over which to transform
+ padding (list): Padding factor for transform. E.g. a padding factor of 2 will discard the upper half of modes after transforming
+
+ Returns:
+ Transformed data
"""
+ # TODO: clean up and try putting more of this in the 1D bases
trfs = {
ChebychevHelper: self._transform_dct,
- Ultraspherical: self._transform_dct,
+ UltrasphericalHelper: self._transform_dct,
FFTHelper: self._transform_fft,
}
@@ -1279,7 +1561,15 @@ class SpectralHelper:
def transform(self, u, axes=None, padding=None):
"""
- Transform all components of the solution
+ Transform all components from physical space to spectral space
+
+ Args:
+ u data to transform:
+ axes (tuple): Axes over which to transform
+ padding (list): Padding factor for transform. E.g. a padding factor of 2 will discard the upper half of modes after transforming
+
+ Returns:
+ Transformed data
"""
axes = tuple(-i - 1 for i in range(self.ndim)[::-1]) if axes is None else axes
padding = (
@@ -1302,6 +1592,7 @@ class SpectralHelper:
return self.xp.stack(result)
def _transform_ifft(self, u, axes, **kwargs):
+ # TODO: clean up and try putting more of this in the 1D bases
ifft = self.get_fft(axes, 'backward', **kwargs)
return ifft(u, axes=axes)
@@ -1309,6 +1600,7 @@ class SpectralHelper:
'''
This will only ever return real values!
'''
+ # TODO: clean up and try putting more of this in the 1D bases
if self.debug:
assert self.xp.allclose(u.imag, 0), 'This function can only handle real input.'
@@ -1357,10 +1649,22 @@ class SpectralHelper:
return v.real
def itransform_single_component(self, u, axes=None, padding=None):
+ """
+ Inverse transform over single component of the solution
+
+ Args:
+ u data to transform:
+ axes (tuple): Axes over which to transform
+ padding (list): Padding factor for transform. E.g. a padding factor of 2 will add as many zeros as there were modes before before transforming
+
+ Returns:
+ Transformed data
+ """
+ # TODO: clean up and try putting more of this in the 1D bases
trfs = {
FFTHelper: self._transform_ifft,
ChebychevHelper: self._transform_idct,
- Ultraspherical: self._transform_idct,
+ UltrasphericalHelper: self._transform_idct,
}
axes = tuple(-i - 1 for i in range(self.ndim)[::-1]) if axes is None else axes
@@ -1417,7 +1721,9 @@ class SpectralHelper:
def get_aligned(self, u, axis_in, axis_out, fft=None, forward=False, **kwargs):
"""
- Realign the data along the axis when using distributed FFTs
+ Realign the data along the axis when using distributed FFTs. `kwargs` will be used to get the correct PFFT
+ object from `mpi4py-fft`, which has suitable transfer classes for the shape of data. Hence, they should include
+ shape especially, if applicable.
Args:
u: The solution
@@ -1479,6 +1785,17 @@ class SpectralHelper:
return _in.redistribute(axis_out)
def itransform(self, u, axes=None, padding=None):
+ """
+ Inverse transform over all components of the solution
+
+ Args:
+ u data to transform:
+ axes (tuple): Axes over which to transform
+ padding (list): Padding factor for transform. E.g. a padding factor of 2 will add as many zeros as there were modes before before transforming
+
+ Returns:
+ Transformed data
+ """
axes = tuple(-i - 1 for i in range(self.ndim)[::-1]) if axes is None else axes
padding = (
[
@@ -1500,11 +1817,27 @@ class SpectralHelper:
return self.xp.stack(result)
def get_local_slice_of_1D_matrix(self, M, axis):
+ """
+ Get the local version of a 1D matrix. When using distributed FFTs, each rank will carry only a subset of modes,
+ which you can sort out via the `SpectralHelper.local_slice` attribute. When constructing a 1D matrix, you can
+ use this method to get the part corresponding to the modes carried by this rank.
+
+ Args:
+ M (sparse matrix): Global 1D matrix you want to get the local version of
+ axis (int): Direction in which you want the local version. You will get the global matrix in other directions. This means slab decomposition only.
+
+ Returns:
+ sparse local matrix
+ """
return M.tocsc()[self.local_slice[axis], self.local_slice[axis]]
def get_filter_matrix(self, axis, **kwargs):
"""
- Get bandpass filter along `axis`
+ Get bandpass filter along `axis`. See the documentation `get_filter_matrix` in the 1D bases for what kwargs are
+ admissible.
+
+ Returns:
+ sparse bandpass matrix
"""
if self.ndim == 1:
return self.axes[0].get_filter_matrix(**kwargs)
@@ -1515,7 +1848,7 @@ class SpectralHelper:
def get_differentiation_matrix(self, axes, **kwargs):
"""
- Get differentiation matrix along specified axis.
+ Get differentiation matrix along specified axis. `kwargs` are forwarded to the 1D base implementation.
Args:
axes (tuple): Axes along which to differentiate.
@@ -1619,7 +1952,7 @@ class SpectralHelper:
return I
- def get_Dirichlet_recombination_matrix(self, axis=-1, **kwargs):
+ def get_Dirichlet_recombination_matrix(self, axis=-1):
"""
Get Dirichlet recombination matrix along axis. Not that it only makes sense in directions discretized with variations of Chebychev bases.
@@ -1633,10 +1966,10 @@ class SpectralHelper:
ndim = len(self.axes)
if ndim == 1:
- C = self.axes[0].get_Dirichlet_recombination_matrix(**kwargs)
+ C = self.axes[0].get_Dirichlet_recombination_matrix()
elif ndim == 2:
axis2 = (axis + 1) % ndim
- C1D = self.axes[axis].get_Dirichlet_recombination_matrix(**kwargs)
+ C1D = self.axes[axis].get_Dirichlet_recombination_matrix()
I1D = self.axes[axis2].get_Id()
@@ -1653,7 +1986,7 @@ class SpectralHelper:
def get_basis_change_matrix(self, axes=None, **kwargs):
"""
Some spectral bases do a change between bases while differentiating. This method returns matrices that changes the basis to whatever you want.
- Refer to the methods of the same name of the 1D bases to learn what parameters you need to pass here.
+ Refer to the methods of the same name of the 1D bases to learn what parameters you need to pass here as `kwargs`.
Args:
axes (tuple): Axes along which to change basis.
diff --git a/pySDC/implementations/problem_classes/Burgers.py b/pySDC/implementations/problem_classes/Burgers.py
index 285d1421616afb02a9c9d86e14e03309f59791ce..ef957f54d0dc33851b0688a73c1399c95ee6d96c 100644
--- a/pySDC/implementations/problem_classes/Burgers.py
+++ b/pySDC/implementations/problem_classes/Burgers.py
@@ -23,7 +23,18 @@ class Burgers1D(GenericSpectralLinear):
dtype_f = imex_mesh
def __init__(self, N=64, epsilon=0.1, BCl=1, BCr=-1, f=0, mode='T2U', **kwargs):
- self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ 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
+ """
+ self._makeAttributeAndRegister('N', 'epsilon', 'BCl', 'BCr', 'f', 'mode', localVars=locals(), readOnly=True)
bases = [{'base': 'cheby', 'N': N}]
components = ['u', 'ux']
@@ -96,7 +107,7 @@ class Burgers1D(GenericSpectralLinear):
def get_fig(self): # pragma: no cover
"""
- Get a figure suitable to plot the solution of this problem
+ Get a figure suitable to plot the solution of this problem.
Returns
-------
@@ -110,7 +121,7 @@ class Burgers1D(GenericSpectralLinear):
def plot(self, u, t=None, fig=None, comp='u'): # pragma: no cover
r"""
- Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
+ Plot the solution.
Parameters
----------
@@ -118,8 +129,8 @@ class Burgers1D(GenericSpectralLinear):
Solution to be plotted
t : float
Time to display at the top of the figure
- fig : matplotlib.pyplot.figure.Figure
- Figure with the correct structure
+ fig : matplotlib.pyplot.figure.Figure, optional
+ Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.
Returns
-------
@@ -157,7 +168,18 @@ class Burgers2D(GenericSpectralLinear):
dtype_f = imex_mesh
def __init__(self, nx=64, nz=64, epsilon=0.1, fux=2, fuz=1, mode='T2U', **kwargs):
- self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ nx (int): Spatial resolution in x direction
+ nz (int): Spatial resolution in z direction
+ epsilon (float): viscosity
+ 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
+ """
+ self._makeAttributeAndRegister('nx', 'nz', 'epsilon', 'fux', 'fuz', 'mode', localVars=locals(), readOnly=True)
bases = [
{'base': 'fft', 'N': nx},
diff --git a/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py b/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
index 7f808a7efe70a2e2c0b228f5e360782dbf11d350..7a3423ec9792dd242d01494711b911879eecf1b2 100644
--- a/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
+++ b/pySDC/implementations/problem_classes/HeatEquation_Chebychev.py
@@ -7,11 +7,27 @@ from pySDC.implementations.problem_classes.generic_spectral import GenericSpectr
class Heat1DChebychev(GenericSpectralLinear):
+ """
+ 1D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1) using a Chebychev spectral method.
+ """
+
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)
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ nvars (int): Resolution
+ a (float): Left BC value
+ b (float): Right BC value
+ f (int): Frequency of the solution
+ nu (float): Diffusion parameter
+ mode ('T2T' or 'T2U'): Mode for Chebychev method.
+
+ """
+ self._makeAttributeAndRegister('nvars', 'a', 'b', 'f', 'nu', 'mode', localVars=locals(), readOnly=True)
bases = [{'base': 'chebychev', 'N': nvars}]
components = ['u', 'ux']
@@ -58,7 +74,18 @@ class Heat1DChebychev(GenericSpectralLinear):
f[iu] = me[iu]
return f
- def u_exact(self, t, noise=0):
+ def u_exact(self, t, noise=0, seed=666):
+ """
+ Get exact solution at time `t`
+
+ Args:
+ t (float): When you want the exact solution
+ noise (float): Add noise of this level
+ seed (int): Random seed for the noise
+
+ Returns:
+ Heat1DChebychev.dtype_u: Exact solution
+ """
xp = self.xp
iu, iux = self.index(self.components)
u = self.spectral.u_init
@@ -76,7 +103,7 @@ class Heat1DChebychev(GenericSpectralLinear):
if noise > 0:
assert t == 0
_noise = self.u_init
- rng = self.xp.random.default_rng(seed=666)
+ rng = self.xp.random.default_rng(seed=seed)
_noise[iu] = rng.normal(size=u[iu].shape)
noise_hat = self.transform(_noise)
low_pass = self.get_filter_matrix(axis=0, kmax=self.nvars - 2)
@@ -95,11 +122,25 @@ class Heat1DChebychev(GenericSpectralLinear):
class Heat1DUltraspherical(GenericSpectralLinear):
+ """
+ 1D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1) using an ultraspherical spectral method.
+ """
+
dtype_u = mesh
dtype_f = mesh
def __init__(self, nvars=128, a=0, b=0, f=1, nu=1.0, **kwargs):
- self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ nvars (int): Resolution
+ a (float): Left BC value
+ b (float): Right BC value
+ f (int): Frequency of the solution
+ nu (float): Diffusion parameter
+ """
+ self._makeAttributeAndRegister('nvars', 'a', 'b', 'f', 'nu', localVars=locals(), readOnly=True)
bases = [{'base': 'ultraspherical', 'N': nvars}]
components = ['u']
@@ -147,7 +188,18 @@ class Heat1DUltraspherical(GenericSpectralLinear):
f[iu][...] = me[iu]
return f
- def u_exact(self, t, noise=0):
+ def u_exact(self, t, noise=0, seed=666):
+ """
+ Get exact solution at time `t`
+
+ Args:
+ t (float): When you want the exact solution
+ noise (float): Add noise of this level
+ seed (int): Random seed for the noise
+
+ Returns:
+ Heat1DUltraspherical.dtype_u: Exact solution
+ """
xp = self.xp
iu = self.index('u')
u = self.spectral.u_init
@@ -161,7 +213,7 @@ class Heat1DUltraspherical(GenericSpectralLinear):
if noise > 0:
assert t == 0
_noise = self.u_init
- rng = self.xp.random.default_rng(seed=666)
+ rng = self.xp.random.default_rng(seed=seed)
_noise[iu] = rng.normal(size=u[iu].shape)
noise_hat = self.transform(_noise)
low_pass = self.get_filter_matrix(axis=0, kmax=self.nvars - 2)
@@ -178,13 +230,34 @@ class Heat1DUltraspherical(GenericSpectralLinear):
class Heat2DChebychev(GenericSpectralLinear):
+ """
+ 2D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1)x(-1,1) using spectral methods based on FFT and Chebychev.
+ """
+
dtype_u = mesh
dtype_f = mesh
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):
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ nx (int): Resolution in x-direction
+ ny (int): Resolution in y-direction
+ base_x (str): Spectral base in x-direction
+ base_y (str): Spectral base in y-direction
+ a (float): BC at y=0 and x=-1
+ b (float): BC at y=0 and x=+1
+ c (float): BC at y=1 and x = -1
+ fx (int): Horizontal frequency of initial conditions
+ fy (int): Vertical frequency of initial conditions
+ nu (float): Diffusion parameter
+ """
assert nx % 2 == 1 or base_x == 'fft'
assert ny % 2 == 1 or base_y == 'fft'
- self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
+ self._makeAttributeAndRegister(
+ 'nx', 'ny', 'base_x', 'base_y', 'a', 'b', 'c', 'fx', 'fy', 'nu', localVars=locals(), readOnly=True
+ )
bases = [{'base': base_x, 'N': nx}, {'base': base_y, 'N': ny}]
components = ['u', 'ux', 'uy']
@@ -265,13 +338,34 @@ class Heat2DChebychev(GenericSpectralLinear):
class Heat2DUltraspherical(GenericSpectralLinear):
+ """
+ 2D Heat equation with Dirichlet Boundary conditions discretized on (-1, 1)x(-1,1) using spectral methods based on FFT and Gegenbauer.
+ """
+
dtype_u = mesh
dtype_f = mesh
def __init__(
self, nx=128, ny=128, base_x='fft', base_y='ultraspherical', a=0, b=0, c=0, fx=1, fy=1, nu=1.0, **kwargs
):
- self._makeAttributeAndRegister(*locals().keys(), localVars=locals(), readOnly=True)
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ nx (int): Resolution in x-direction
+ ny (int): Resolution in y-direction
+ base_x (str): Spectral base in x-direction
+ base_y (str): Spectral base in y-direction
+ a (float): BC at y=0 and x=-1
+ b (float): BC at y=0 and x=+1
+ c (float): BC at y=1 and x = -1
+ fx (int): Horizontal frequency of initial conditions
+ fy (int): Vertical frequency of initial conditions
+ nu (float): Diffusion parameter
+ """
+ self._makeAttributeAndRegister(
+ 'nx', 'ny', 'base_x', 'base_y', 'a', 'b', 'c', 'fx', 'fy', 'nu', localVars=locals(), readOnly=True
+ )
bases = [{'base': base_x, 'N': nx}, {'base': base_y, 'N': ny}]
components = ['u']
diff --git a/pySDC/implementations/problem_classes/RayleighBenard.py b/pySDC/implementations/problem_classes/RayleighBenard.py
index b521d01d8d45a2a28615087aa32bd5132f1a1b76..3a16f2ad063ccf5b968c653602e71ffaaf4d2734 100644
--- a/pySDC/implementations/problem_classes/RayleighBenard.py
+++ b/pySDC/implementations/problem_classes/RayleighBenard.py
@@ -10,7 +10,30 @@ from pySDC.implementations.convergence_controller_classes.check_convergence impo
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.
+ Rayleigh-Benard Convection is a variation of incompressible Navier-Stokes.
+
+ The equations we solve are
+
+ u_x + v_z = 0
+ T_t - kappa (T_xx + T_zz) = -uT_x - vT_z
+ u_t - nu (u_xx + u_zz) + p_x = -uu_x - vu_z
+ v_t - nu (v_xx + v_zz) + p_z - T = -uv_x - vv_z
+
+ with u the horizontal velocity, v the vertical velocity (in z-direction), T the temperature, p the pressure, indices
+ denoting derivatives, kappa=(Rayleigh * Prandl)**(-1/2) and nu = (Rayleigh / Prandl)**(-1/2). Everything on the left
+ hand side, that is the viscous part, the pressure gradient and the buoyancy due to temperature are treated
+ implicitly, while the non-linear convection part on the right hand side is integrated explicitly.
+
+ The domain, vertical boundary conditions and pressure gauge are
+
+ Omega = [0, 8) x (-1, 1)
+ T(z=+1) = 0
+ T(z=-1) = 2
+ u(z=+-1) = v(z=+-1) = 0
+ integral over p = 0
+
+ The spectral discretization uses FFT horizontally, implying periodic BCs, and an ultraspherical method vertically to
+ facilitate the Dirichlet BCs.
Parameters:
Prandl (float): Prandl number
@@ -29,13 +52,25 @@ class RayleighBenard(GenericSpectralLinear):
self,
Prandl=1,
Rayleigh=2e6,
- nx=257,
+ nx=256,
nz=64,
BCs=None,
dealiasing=3 / 2,
comm=None,
**kwargs,
):
+ """
+ Constructor. `kwargs` are forwarded to parent class constructor.
+
+ Args:
+ Prandl (float): Prandtl number
+ Rayleigh (float): Rayleigh number
+ nx (int): Resolution in x-direction
+ nz (int): Resolution in z direction
+ BCs (dict): Vertical boundary conditions
+ dealiasing (float): Dealiasing for evaluating the non-linear part in real space
+ comm (mpi4py.Intracomm): Space communicator
+ """
BCs = {} if BCs is None else BCs
BCs = {
'T_top': 0,
@@ -108,8 +143,10 @@ class RayleighBenard(GenericSpectralLinear):
M_lhs = {i: {i: U02 @ Id} for i in ['u', 'v', 'T']}
self.setup_M(M_lhs)
+ # Prepare going from second (first for divergence free equation) derivative basis back to Chebychov-T
self.base_change = self._setup_operator({**{comp: {comp: S2} for comp in ['u', 'v', 'T']}, 'p': {'p': S1}})
+ # BCs
self.add_BC(
component='p', equation='p', axis=1, v=self.BCs['p_integral'], kind='integral', line=-1, scalar=True
)
@@ -238,7 +275,7 @@ class RayleighBenard(GenericSpectralLinear):
def plot(self, u, t=None, fig=None, quantity='T'): # pragma: no cover
r"""
- Plot the solution. Please supply a figure with the same structure as returned by ``self.get_fig``.
+ Plot the solution.
Parameters
----------
@@ -247,7 +284,9 @@ class RayleighBenard(GenericSpectralLinear):
t : float
Time to display at the top of the figure
fig : matplotlib.pyplot.figure.Figure
- Figure with the correct structure
+ Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.
+ quantity : (str)
+ quantity you want to plot
Returns
-------
diff --git a/pySDC/implementations/problem_classes/generic_spectral.py b/pySDC/implementations/problem_classes/generic_spectral.py
index 3b523b2b15730dc573497553125ca55fcf4ebbd4..04f88823fc4f2889c5fe21400721720f4d934b8c 100644
--- a/pySDC/implementations/problem_classes/generic_spectral.py
+++ b/pySDC/implementations/problem_classes/generic_spectral.py
@@ -6,9 +6,26 @@ from pySDC.core.errors import ParameterError
class GenericSpectralLinear(Problem):
"""
- Generic class to solve problems of the form M u_t + L u = y, with mass matrix M, linear operator L and some right hand side y using spectral methods.
+ Generic class to solve problems of the form M u_t + L u = y, with mass matrix M, linear operator L and some right
+ hand side y using spectral methods.
L may contain algebraic conditions, as long as (M + dt L) is invertible.
+ Note that the `__getattr__` method is overloaded to pass requests on to the spectral helper if they are not
+ attributes of this class itself. For instance, you can add a BC by calling `self.spectral.add_BC` or equivalently
+ `self.add_BC`.
+
+ You can port problems derived from this more or less seamlessly to GPU by using the numerical libraries that are
+ class attributes of the spectral helper. This class will automatically switch the datatype using the `setup_GPU` class method.
+
+ Attributes:
+ spectral (pySDC.helpers.spectral_helper.SpectralHelper): Spectral helper
+ work_counters (dict): Dictionary for counting work
+ cached_factorizations (dict): Dictionary of cached matrix factorizations for solving
+ L (sparse matrix): Linear operator
+ M (sparse matrix): Mass matrix
+ diff_mask (list): Mask for separating differential and algebraic terms
+ Pl (sparse matrix): Left preconditioner
+ Pr (sparse matrix): Right preconditioner
"""
@classmethod
@@ -40,8 +57,7 @@ class GenericSpectralLinear(Problem):
max_cached_factorizations=12,
spectral_space=True,
real_spectral_coefficients=False,
- *args,
- **kwargs,
+ debug=False,
):
"""
Base class for problems discretized with spectral methods.
@@ -58,6 +74,7 @@ class GenericSpectralLinear(Problem):
max_cached_factorizations (int): Number of matrix decompositions to cache before starting eviction
spectral_space (bool): If yes, the solution will not be transformed back after solving and evaluating the RHS, and is expected as input in spectral space to these functions
real_spectral_coefficients (bool): If yes, allow only real values in spectral space, otherwise, allow complex.
+ debug (bool): Make additional tests at extra computational cost
"""
solver_args = {} if solver_args is None else solver_args
self._makeAttributeAndRegister(
@@ -70,9 +87,10 @@ class GenericSpectralLinear(Problem):
'comm',
'spectral_space',
'real_spectral_coefficients',
+ 'debug',
localVars=locals(),
)
- self.spectral = SpectralHelper(comm=comm, useGPU=useGPU)
+ self.spectral = SpectralHelper(comm=comm, useGPU=useGPU, debug=debug)
if useGPU:
self.setup_GPU()
@@ -93,6 +111,15 @@ class GenericSpectralLinear(Problem):
self.cached_factorizations = {}
def __getattr__(self, name):
+ """
+ Pass requests on to the helper if they are not directly attributes of this class for convenience.
+
+ Args:
+ name (str): Name of the attribute you want
+
+ Returns:
+ request
+ """
return getattr(self.spectral, name)
def _setup_operator(self, LHS):
@@ -117,7 +144,7 @@ class GenericSpectralLinear(Problem):
The argument is meant to be a dictionary with the line you want to write the equation in as the key and the relationship between components as another dictionary. For instance, you can add an algebraic condition capturing a first derivative relationship between u and ux as follows:
```
- Dx = self.get_self.get_differentiation_matrix(axes=(0,))
+ Dx = self.get_differentiation_matrix(axes=(0,))
I = self.get_Id()
LHS = {'ux': {'u': Dx, 'ux': -I}}
self.setup_L(LHS)
@@ -134,16 +161,16 @@ class GenericSpectralLinear(Problem):
'''
Setup mass matrix, see documentation of ``GenericSpectralLinear.setup_L``.
'''
- self.diff_index = list(LHS.keys())
- self.diff_mask = [me in self.diff_index for me in self.components]
+ diff_index = list(LHS.keys())
+ self.diff_mask = [me in diff_index for me in self.components]
self.M = self._setup_operator(LHS)
def setup_preconditioner(self, Dirichlet_recombination=True, left_preconditioner=True):
"""
- Get left and right precondioners. A right preconditioner of D2T will result in Dirichlet recombination. 10/10 would recommend!
+ Get left and right preconditioners.
Args:
- Dirichlet_recombination (bool): Basis conversion for right precondioner
+ Dirichlet_recombination (bool): Basis conversion for right preconditioner. Useful for Chebychev and Ultraspherical methods. 10/10 would recommend.
left_preconditioner (bool): If True, it will interleave the variables and reverse the Kronecker product
"""
sp = self.spectral.sparse_lib
@@ -177,9 +204,14 @@ class GenericSpectralLinear(Problem):
def solve_system(self, rhs, dt, u0=None, *args, skip_itransform=False, **kwargs):
"""
- Solve (M + dt*L)u=rhs. This requires that you setup the operators before using the functions ``GenericSpectralLinear.setup_L`` and ``GenericSpectralLinear.setup_M``. Note that the mass matrix need not be invertible, as long as (M + dt*L) is. This allows to solve some differential algebraic equations.
+ Do an implicit Euler step to solve M u_t + Lu = rhs, with M the mass matrix and L the linear operator as setup by
+ ``GenericSpectralLinear.setup_L`` and ``GenericSpectralLinear.setup_M``.
+
+ The implicit Euler step is (M - dt L) u = M rhs. Note that M need not be invertible as long as (M + dt*L) is.
+ This means solving with dt=0 to mimic explicit methods does not work for all problems, in particular simple DAEs.
- Note that in implicit Euler, the right hand side will be composed of the initial conditions. We don't want that in lines that don't depend on time. Therefore, we multiply the right hand side by the mass matrix. This means you can only do algebraic conditions that add up to zero. But you can easily overload this function with something more generic if needed.
+ Note that by putting M rhs on the right hand side, this function can only solve algebraic conditions equal to
+ zero. If you want something else, it should be easy to overload this function.
"""
sp = self.spectral.sparse_lib
diff --git a/pySDC/tests/test_helpers/test_spectral_helper.py b/pySDC/tests/test_helpers/test_spectral_helper.py
index 433fb03c8de2e5845271afda8426639672fcf198..627632b14b98132b28e1318eaf7465b3dafdada8 100644
--- a/pySDC/tests/test_helpers/test_spectral_helper.py
+++ b/pySDC/tests/test_helpers/test_spectral_helper.py
@@ -19,7 +19,7 @@ def test_integration_matrix2D(nx, nz, variant, axes, useMPI=False, **kwargs):
helper = SpectralHelper(comm=comm, debug=True)
helper.add_axis(base='fft', N=nx)
- helper.add_axis(base='cheby', N=nz, mode=variant)
+ helper.add_axis(base='cheby', N=nz)
helper.setup_fft()
Z, X = helper.get_grid()
@@ -65,7 +65,7 @@ def test_differentiation_matrix2D(nx, nz, variant, axes, bx, bz, useMPI=False, *
helper = SpectralHelper(comm=comm, debug=True)
helper.add_axis(base=bx, N=nx)
- helper.add_axis(base=bz, N=nz, mode=variant)
+ helper.add_axis(base=bz, N=nz)
helper.setup_fft()
Z, X = helper.get_grid()
@@ -120,7 +120,7 @@ def test_identity_matrix2D(nx, nz, variant, bx, useMPI=False, **kwargs):
helper = SpectralHelper(comm=comm, debug=True)
helper.add_axis(base=bx, N=nx)
- helper.add_axis(base='cheby', N=nz, mode=variant)
+ helper.add_axis(base='cheby', N=nz)
helper.setup_fft()
Z, X = helper.get_grid()
@@ -226,8 +226,6 @@ def _test_transform_dealias(
u2_hat_expect[0][xp.logical_and(xp.abs(Kx) == 0, Kz == 0)] += 2 / nx
u_expect[0] += np.cos(f * X) * 2 / nx
u_expect_pad[0] += np.cos(f * X_pad) * 2 / nx
- # u2_expect = u_expect**2
- # u2_expect_pad = u_expect_pad**2
elif axis == -1:
f = nz // 2 + 1
u_hat[0][xp.logical_and(xp.abs(Kz) == f, Kx == 0)] += 1
@@ -238,8 +236,6 @@ def _test_transform_dealias(
coef[f] = 1 / nx
u_expect[0] = np.polynomial.Chebyshev(coef)(Z)
u_expect_pad[0] = np.polynomial.Chebyshev(coef)(Z_pad)
- # u2_expect = u_expect**2
- # u2_expect_pad = u_expect_pad**2
elif axis in [(-1, -2), (-2, -1)]:
fx = nx // 3
fz = nz // 2 + 1
@@ -257,9 +253,7 @@ def _test_transform_dealias(
coef[fz] = 1 / nx
u_expect[0] = np.cos(fx * X) * 2 / nx + np.polynomial.Chebyshev(coef)(Z)
- # u2_expect = u_expect**2
u_expect_pad[0] = np.cos(fx * X_pad) * 2 / nx + np.polynomial.Chebyshev(coef)(Z_pad)
- # u2_expect_pad = u_expect_pad**2
else:
raise NotImplementedError
@@ -280,60 +274,6 @@ def _test_transform_dealias(
assert not np.allclose(u2_hat_expect, helper.transform(u2)), 'Test is too boring, no dealiasing needed'
-# @pytest.mark.base
-# @pytest.mark.parametrize('bx', ['fft', 'cheby'])
-# @pytest.mark.parametrize('bz', ['fft', 'cheby'])
-# @pytest.mark.parametrize('axis', [0, 1])
-# def test_transform_dealias_back_and_forth(
-# bx,
-# bz,
-# axis,
-# nx=2**1 + 1,
-# nz=2**1 + 1,
-# padding=3 / 2,
-# axes=(
-# -2,
-# -1,
-# ),
-# useMPI=False,
-# **kwargs,
-# ):
-# import numpy as np
-# from pySDC.helpers.spectral_helper import SpectralHelper
-#
-# if useMPI:
-# from mpi4py import MPI
-#
-# comm = MPI.COMM_WORLD
-# rank = comm.rank
-# else:
-# comm = None
-#
-# helper = SpectralHelper(comm=comm, debug=True)
-# helper.add_axis(base=bx, N=nx)
-# helper.add_axis(base=bz, N=nz)
-# helper.setup_fft()
-# xp = helper.xp
-#
-# shape = helper.shape
-# helper_padded = helper.get_zero_padded_version(axis=axis, padding=padding)
-#
-# u = helper.u_init
-# u[:] = xp.random.rand(*shape)
-#
-# u_hat = helper.transform(u)
-# u_hat_padded = helper_padded.get_padded(u_hat)
-#
-# u_padded = helper_padded.itransform(u_hat_padded)
-#
-# u_2_padded_hat = helper_padded.transform(u_padded)
-# u_2_hat = helper_padded.get_unpadded(u_2_padded_hat)
-#
-# assert xp.allclose(u_2_padded_hat.real, u_hat_padded.real)
-# assert xp.allclose(u_2_padded_hat.imag, u_hat_padded.imag)
-# assert xp.allclose(u_2_hat, u_hat)
-
-
@pytest.mark.base
@pytest.mark.parametrize('nx', [3, 8])
@pytest.mark.parametrize('nz', [3, 8])
@@ -493,12 +433,6 @@ def test_tau_method(bc, N, bc_val, kind='Dirichlet'):
rhs_hat = helper.transform(rhs, axes=(-1,))
sol_hat = sp.linalg.spsolve(A, rhs_hat.flatten())
- # print(_sol_hat)
- # sol_hat = helper.remove_tau_terms(_sol_hat)
- # print(A.toarray())
- # import matplotlib.pyplot as plt
- # plt.imshow(A.real/np.abs(A))
- # plt.show()
sol_poly = np.polynomial.Chebyshev(sol_hat)
d_sol_poly = sol_poly.deriv(1)
@@ -538,7 +472,7 @@ def test_tau_method2D(variant, nz, nx, bc_val, bc=-1, useMPI=False, plotting=Fal
helper = SpectralHelper(comm=comm, debug=True)
helper.add_axis('fft', N=nx)
- helper.add_axis('cheby', N=nz, mode=variant)
+ helper.add_axis('cheby', N=nz)
helper.add_component(['u'])
helper.setup_fft()
@@ -563,8 +497,6 @@ def test_tau_method2D(variant, nz, nx, bc_val, bc=-1, useMPI=False, plotting=Fal
# prepare system to solve
A = helper.put_BCs_in_matrix(A)
- # rhs = helper.put_BCs_in_rhs(helper.u_init)
- # rhs_hat = helper.transform(rhs, axes=(-1, -2))
rhs_hat = helper.put_BCs_in_rhs_hat(helper.u_init_forward)
# solve the system
@@ -584,7 +516,6 @@ def test_tau_method2D(variant, nz, nx, bc_val, bc=-1, useMPI=False, plotting=Fal
import matplotlib.pyplot as plt
im = plt.pcolormesh(X, Z, sol[0])
- # im = plt.pcolormesh(X, Z, error)
plt.colorbar(im)
plt.xlabel('x')
plt.ylabel('t')
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 bd7bd2eb3489f31d8e47e3f8d51926eebc5524ce..d5ba10afbf18fab9eb031bd68c35597bb6601106 100644
--- a/pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py
+++ b/pySDC/tests/test_helpers/test_spectral_helper_1d_ultraspherical.py
@@ -6,9 +6,9 @@ import pytest
@pytest.mark.parametrize('p', [1, 2, 3, 4])
def test_differentiation_matrix(N, p):
import numpy as np
- from pySDC.helpers.spectral_helper import Ultraspherical
+ from pySDC.helpers.spectral_helper import UltrasphericalHelper
- helper = Ultraspherical(N)
+ helper = UltrasphericalHelper(N)
x = helper.get_1dgrid()
coeffs = np.random.random(N)
@@ -25,9 +25,9 @@ def test_differentiation_matrix(N, p):
@pytest.mark.parametrize('N', [4, 7, 32])
def test_integration(N):
import numpy as np
- from pySDC.helpers.spectral_helper import Ultraspherical
+ from pySDC.helpers.spectral_helper import UltrasphericalHelper
- helper = Ultraspherical(N)
+ helper = UltrasphericalHelper(N)
coeffs = np.random.random(N)
coeffs[-1] = 0
@@ -47,12 +47,12 @@ def test_integration(N):
def test_poisson_problem(N, deg, Dirichlet_recombination):
import numpy as np
import scipy.sparse as sp
- from pySDC.helpers.spectral_helper import Ultraspherical
+ from pySDC.helpers.spectral_helper import UltrasphericalHelper
a = 0
b = 4
- helper = Ultraspherical(N)
+ helper = UltrasphericalHelper(N)
x = helper.get_1dgrid()
f = x**deg * (deg + 1) * (deg + 2) * (a - b) / 2
diff --git a/pySDC/tests/test_problems/test_Burgers.py b/pySDC/tests/test_problems/test_Burgers.py
index 12bcedc3a13323961f4a9efca6b620e69fa28cff..ca540e0ab5ef350ef9912f3c90d33f0093a28906 100644
--- a/pySDC/tests/test_problems/test_Burgers.py
+++ b/pySDC/tests/test_problems/test_Burgers.py
@@ -144,7 +144,6 @@ def test_Burgers2D_solver(mode, nx=2**6, nz=2**6, plotting=False):
mode=mode,
comm=comm,
left_preconditioner=True,
- right_preconditioning='D2T',
)
iu, iv, iux, iuz, ivx, ivz = P.index(P.components)
diff --git a/pySDC/tests/test_problems/test_heat_chebychev.py b/pySDC/tests/test_problems/test_heat_chebychev.py
index acc4b4169f4fafd546ebe237b6ac4e68f268d29f..8a5f150d7ba3a83ec28b437b9b4e4bbb68145103 100644
--- a/pySDC/tests/test_problems/test_heat_chebychev.py
+++ b/pySDC/tests/test_problems/test_heat_chebychev.py
@@ -23,7 +23,6 @@ def test_heat1d_chebychev(a, b, f, noise, use_ultraspherical, spectral_space, nv
f=f,
nu=1e-2,
left_preconditioner=False,
- right_preconditioning='T2T',
debug=True,
spectral_space=spectral_space,
)