Skip to content

Chebyshev Filters

ChebyshevFilter1D

Bases: Module

1D Chebyshev spectral filter for smoothing and numerical stabilisation.

Mathematical Formulation

Filters are applied as a multiplicative mask in Chebyshev-coefficient space:

ãₖ = F(k) · aₖ

where aₖ are Chebyshev coefficients of u(x) = Σₖ aₖ Tₖ(x/L).

Exponential filter F(k) = exp(−α (k/kₘₐₓ)ᵖ) # kₘₐₓ = N for GL, N−1 for Gauss

Hyperviscosity filter F(k) = exp(−ν_h kᵖ Δt)

Attributes

grid : ChebyshevGrid1D Underlying 1D Chebyshev grid (provides the forward/inverse transform).

Examples

Smooth a noisy field with the default machine-epsilon exponential filter:

import jax.numpy as jnp grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0) flt = ChebyshevFilter1D(grid=grid) u = jnp.sin(jnp.pi * grid.x) + 1e-2 * jnp.cos(31 * jnp.pi * grid.x) u_smooth = flt.exponential_filter(u) # high mode is suppressed u_damped = flt.hyperviscosity(u, nu_hyper=1e-4, dt=0.01, power=4)

Source code in spectraldiffx/_src/chebyshev/filters.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
class ChebyshevFilter1D(eqx.Module):
    """1D Chebyshev spectral filter for smoothing and numerical stabilisation.

    Mathematical Formulation
    ------------------------
    Filters are applied as a multiplicative mask in Chebyshev-coefficient space:

        ãₖ = F(k) · aₖ

    where aₖ are Chebyshev coefficients of u(x) = Σₖ aₖ Tₖ(x/L).

    Exponential filter
        F(k) = exp(−α (k/kₘₐₓ)ᵖ)          # kₘₐₓ = N for GL, N−1 for Gauss

    Hyperviscosity filter
        F(k) = exp(−ν_h kᵖ Δt)

    Attributes
    ----------
    grid : ChebyshevGrid1D
        Underlying 1D Chebyshev grid (provides the forward/inverse transform).

    Examples
    --------
    Smooth a noisy field with the default machine-epsilon exponential filter:

    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid1D.from_N_L(N=32, L=1.0)
    >>> flt = ChebyshevFilter1D(grid=grid)
    >>> u = jnp.sin(jnp.pi * grid.x) + 1e-2 * jnp.cos(31 * jnp.pi * grid.x)
    >>> u_smooth = flt.exponential_filter(u)  # high mode is suppressed
    >>> u_damped = flt.hyperviscosity(u, nu_hyper=1e-4, dt=0.01, power=4)
    """

    grid: ChebyshevGrid1D

    def exponential_filter(
        self,
        u: Num[Array, "Npts"],
        alpha: float = 36.0,
        power: int = 16,
        spectral: bool = False,
    ) -> Num[Array, "Npts"]:
        """Apply an exponential cut-off filter in Chebyshev mode space.

            F(k) = exp(−α (k/kₘₐₓ)ᵖ)

        Near unity for low k and falls off sharply near kₘₐₓ.  The default
        α = 36 gives F(kₘₐₓ) ≈ exp(−36) ≈ 2·10⁻¹⁶ (≈ double-precision ε).

        Parameters
        ----------
        u : Num[Array, "Npts"]
            Physical-space field, or Chebyshev coefficients if ``spectral=True``.
        alpha : float
            Damping strength. Default 36.0.
        power : int
            Sharpening exponent (even integer). Default 16.  Must be > 0.
        spectral : bool
            If ``True``, ``u`` is treated as Chebyshev coefficients.

        Returns
        -------
        Num[Array, "Npts"]
            Filtered field (physical if ``spectral=False``, else coefficients).
        """
        if alpha < 0:
            raise ValueError(f"alpha must be >= 0, got {alpha}")
        if power <= 0:
            raise ValueError(f"power must be > 0, got {power}")
        a = u if spectral else self.grid.transform(u)
        N = self.grid.N
        n_modes = N + 1 if self.grid.node_type == "gauss-lobatto" else N
        # Normalise by the highest mode index so F(kₘₐₓ) = exp(−α).
        k_max = max(1, n_modes - 1)
        real_dtype = _coeff_dtype(a)
        k = jnp.arange(n_modes, dtype=real_dtype)
        filter_mask = jnp.exp(-alpha * (k / k_max) ** power)
        a_f = a * filter_mask
        return a_f if spectral else self.grid.transform(a_f, inverse=True)

    def hyperviscosity(
        self,
        u: Num[Array, "Npts"],
        nu_hyper: float,
        dt: float,
        power: int = 4,
        spectral: bool = False,
    ) -> Num[Array, "Npts"]:
        """Apply hyperviscous damping in Chebyshev mode space.

            F(k) = exp(−ν_h kᵖ Δt)

        Simulates high-order diffusion: ∂u/∂t = (−1)^(p/2) ν_h ∂ᵖu/∂xᵖ.

        Parameters
        ----------
        u : Num[Array, "Npts"]
            Physical-space field or Chebyshev coefficients.
        nu_hyper : float
            Hyperviscosity coefficient (≥ 0).
        dt : float
            Time step for the damping (≥ 0).
        power : int
            Diffusion order.  Default 4 (biharmonic).  Must be > 0.
        spectral : bool
            If ``True``, ``u`` is treated as Chebyshev coefficients.

        Returns
        -------
        Num[Array, "Npts"]
            Damped field or coefficients.
        """
        if nu_hyper < 0:
            raise ValueError(f"nu_hyper must be >= 0, got {nu_hyper}")
        if power <= 0:
            raise ValueError(f"power must be > 0, got {power}")
        a = u if spectral else self.grid.transform(u)
        N = self.grid.N
        n_modes = N + 1 if self.grid.node_type == "gauss-lobatto" else N
        real_dtype = _coeff_dtype(a)
        k = jnp.arange(n_modes, dtype=real_dtype)
        filter_mask = jnp.exp(-nu_hyper * k**power * dt)
        a_f = a * filter_mask
        return a_f if spectral else self.grid.transform(a_f, inverse=True)

Functions

exponential_filter(u, alpha=36.0, power=16, spectral=False)

Apply an exponential cut-off filter in Chebyshev mode space.

F(k) = exp(−α (k/kₘₐₓ)ᵖ)

Near unity for low k and falls off sharply near kₘₐₓ. The default α = 36 gives F(kₘₐₓ) ≈ exp(−36) ≈ 2·10⁻¹⁶ (≈ double-precision ε).

Parameters

u : Num[Array, "Npts"] Physical-space field, or Chebyshev coefficients if spectral=True. alpha : float Damping strength. Default 36.0. power : int Sharpening exponent (even integer). Default 16. Must be > 0. spectral : bool If True, u is treated as Chebyshev coefficients.

Returns

Num[Array, "Npts"] Filtered field (physical if spectral=False, else coefficients).

Source code in spectraldiffx/_src/chebyshev/filters.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def exponential_filter(
    self,
    u: Num[Array, "Npts"],
    alpha: float = 36.0,
    power: int = 16,
    spectral: bool = False,
) -> Num[Array, "Npts"]:
    """Apply an exponential cut-off filter in Chebyshev mode space.

        F(k) = exp(−α (k/kₘₐₓ)ᵖ)

    Near unity for low k and falls off sharply near kₘₐₓ.  The default
    α = 36 gives F(kₘₐₓ) ≈ exp(−36) ≈ 2·10⁻¹⁶ (≈ double-precision ε).

    Parameters
    ----------
    u : Num[Array, "Npts"]
        Physical-space field, or Chebyshev coefficients if ``spectral=True``.
    alpha : float
        Damping strength. Default 36.0.
    power : int
        Sharpening exponent (even integer). Default 16.  Must be > 0.
    spectral : bool
        If ``True``, ``u`` is treated as Chebyshev coefficients.

    Returns
    -------
    Num[Array, "Npts"]
        Filtered field (physical if ``spectral=False``, else coefficients).
    """
    if alpha < 0:
        raise ValueError(f"alpha must be >= 0, got {alpha}")
    if power <= 0:
        raise ValueError(f"power must be > 0, got {power}")
    a = u if spectral else self.grid.transform(u)
    N = self.grid.N
    n_modes = N + 1 if self.grid.node_type == "gauss-lobatto" else N
    # Normalise by the highest mode index so F(kₘₐₓ) = exp(−α).
    k_max = max(1, n_modes - 1)
    real_dtype = _coeff_dtype(a)
    k = jnp.arange(n_modes, dtype=real_dtype)
    filter_mask = jnp.exp(-alpha * (k / k_max) ** power)
    a_f = a * filter_mask
    return a_f if spectral else self.grid.transform(a_f, inverse=True)

hyperviscosity(u, nu_hyper, dt, power=4, spectral=False)

Apply hyperviscous damping in Chebyshev mode space.

F(k) = exp(−ν_h kᵖ Δt)

Simulates high-order diffusion: ∂u/∂t = (−1)^(p/2) ν_h ∂ᵖu/∂xᵖ.

Parameters

u : Num[Array, "Npts"] Physical-space field or Chebyshev coefficients. nu_hyper : float Hyperviscosity coefficient (≥ 0). dt : float Time step for the damping (≥ 0). power : int Diffusion order. Default 4 (biharmonic). Must be > 0. spectral : bool If True, u is treated as Chebyshev coefficients.

Returns

Num[Array, "Npts"] Damped field or coefficients.

Source code in spectraldiffx/_src/chebyshev/filters.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def hyperviscosity(
    self,
    u: Num[Array, "Npts"],
    nu_hyper: float,
    dt: float,
    power: int = 4,
    spectral: bool = False,
) -> Num[Array, "Npts"]:
    """Apply hyperviscous damping in Chebyshev mode space.

        F(k) = exp(−ν_h kᵖ Δt)

    Simulates high-order diffusion: ∂u/∂t = (−1)^(p/2) ν_h ∂ᵖu/∂xᵖ.

    Parameters
    ----------
    u : Num[Array, "Npts"]
        Physical-space field or Chebyshev coefficients.
    nu_hyper : float
        Hyperviscosity coefficient (≥ 0).
    dt : float
        Time step for the damping (≥ 0).
    power : int
        Diffusion order.  Default 4 (biharmonic).  Must be > 0.
    spectral : bool
        If ``True``, ``u`` is treated as Chebyshev coefficients.

    Returns
    -------
    Num[Array, "Npts"]
        Damped field or coefficients.
    """
    if nu_hyper < 0:
        raise ValueError(f"nu_hyper must be >= 0, got {nu_hyper}")
    if power <= 0:
        raise ValueError(f"power must be > 0, got {power}")
    a = u if spectral else self.grid.transform(u)
    N = self.grid.N
    n_modes = N + 1 if self.grid.node_type == "gauss-lobatto" else N
    real_dtype = _coeff_dtype(a)
    k = jnp.arange(n_modes, dtype=real_dtype)
    filter_mask = jnp.exp(-nu_hyper * k**power * dt)
    a_f = a * filter_mask
    return a_f if spectral else self.grid.transform(a_f, inverse=True)

ChebyshevFilter2D

Bases: Module

2D Chebyshev spectral filter on [−Lx, Lx] × [−Ly, Ly].

Applies separable 1D exponential or hyperviscosity kernels:

F(kₓ, kᵧ) = Fₓ(kₓ) · Fᵧ(kᵧ)

Attributes

grid : ChebyshevGrid2D Underlying 2D Chebyshev grid.

Examples

import jax.numpy as jnp grid = ChebyshevGrid2D.from_N_L(Nx=24, Ny=24, Lx=1.0, Ly=1.0) flt = ChebyshevFilter2D(grid=grid) X, Y = grid.X u = jnp.sin(2 * jnp.pi * X) * jnp.cos(3 * jnp.pi * Y) u_smooth = flt.exponential_filter(u, alpha=20.0, power=8)

Source code in spectraldiffx/_src/chebyshev/filters.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
class ChebyshevFilter2D(eqx.Module):
    """2D Chebyshev spectral filter on [−Lx, Lx] × [−Ly, Ly].

    Applies separable 1D exponential or hyperviscosity kernels:

        F(kₓ, kᵧ) = Fₓ(kₓ) · Fᵧ(kᵧ)

    Attributes
    ----------
    grid : ChebyshevGrid2D
        Underlying 2D Chebyshev grid.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = ChebyshevGrid2D.from_N_L(Nx=24, Ny=24, Lx=1.0, Ly=1.0)
    >>> flt = ChebyshevFilter2D(grid=grid)
    >>> X, Y = grid.X
    >>> u = jnp.sin(2 * jnp.pi * X) * jnp.cos(3 * jnp.pi * Y)
    >>> u_smooth = flt.exponential_filter(u, alpha=20.0, power=8)
    """

    grid: ChebyshevGrid2D

    def exponential_filter(
        self,
        u: Num[Array, "Nypts Nxpts"],
        alpha: float = 36.0,
        power: int = 16,
        spectral: bool = False,
    ) -> Num[Array, "Nypts Nxpts"]:
        """Separable 2D exponential filter.

            F(kₓ, kᵧ) = exp(−α (kₓ/kₓ_max)ᵖ) · exp(−α (kᵧ/kᵧ_max)ᵖ)

        Parameters
        ----------
        u : Num[Array, "Nypts Nxpts"]
            Physical-space field or Chebyshev coefficients.
        alpha : float
            Damping strength (≥ 0).  Default 36.0.
        power : int
            Sharpening exponent (> 0).  Default 16.
        spectral : bool
            If ``True``, ``u`` is treated as spectral coefficients.

        Returns
        -------
        Num[Array, "Nypts Nxpts"]
        """
        if alpha < 0:
            raise ValueError(f"alpha must be >= 0, got {alpha}")
        if power <= 0:
            raise ValueError(f"power must be > 0, got {power}")
        a = u if spectral else self.grid.transform(u)
        Nx, Ny = self.grid.Nx, self.grid.Ny
        nx_modes = Nx + 1 if self.grid.node_type == "gauss-lobatto" else Nx
        ny_modes = Ny + 1 if self.grid.node_type == "gauss-lobatto" else Ny
        # Normalise by the highest mode index; derive dtype from coefficients.
        kx_max = max(1, nx_modes - 1)
        ky_max = max(1, ny_modes - 1)
        real_dtype = _coeff_dtype(a)
        kx = jnp.arange(nx_modes, dtype=real_dtype)
        ky = jnp.arange(ny_modes, dtype=real_dtype)
        Fx = jnp.exp(-alpha * (kx / kx_max) ** power)
        Fy = jnp.exp(-alpha * (ky / ky_max) ** power)
        # Separable filter: F = Fᵧ[:, None] · Fₓ[None, :]
        a_f = a * Fy[:, None] * Fx[None, :]
        return a_f if spectral else self.grid.transform(a_f, inverse=True)

    def hyperviscosity(
        self,
        u: Num[Array, "Nypts Nxpts"],
        nu_hyper: float,
        dt: float,
        power: int = 4,
        spectral: bool = False,
    ) -> Num[Array, "Nypts Nxpts"]:
        """Separable 2D hyperviscosity filter.

            F(kₓ, kᵧ) = exp(−ν_h kₓᵖ Δt) · exp(−ν_h kᵧᵖ Δt)

        Parameters
        ----------
        u : Num[Array, "Nypts Nxpts"]
            Physical-space field or Chebyshev coefficients.
        nu_hyper : float
            Hyperviscosity coefficient (≥ 0).
        dt : float
            Time step for the damping (≥ 0).
        power : int
            Diffusion order (> 0).  Default 4 (biharmonic).
        spectral : bool
            If ``True``, ``u`` is treated as spectral coefficients.

        Returns
        -------
        Num[Array, "Nypts Nxpts"]
        """
        if nu_hyper < 0:
            raise ValueError(f"nu_hyper must be >= 0, got {nu_hyper}")
        if power <= 0:
            raise ValueError(f"power must be > 0, got {power}")
        a = u if spectral else self.grid.transform(u)
        Nx, Ny = self.grid.Nx, self.grid.Ny
        nx_modes = Nx + 1 if self.grid.node_type == "gauss-lobatto" else Nx
        ny_modes = Ny + 1 if self.grid.node_type == "gauss-lobatto" else Ny
        real_dtype = _coeff_dtype(a)
        kx = jnp.arange(nx_modes, dtype=real_dtype)
        ky = jnp.arange(ny_modes, dtype=real_dtype)
        # Separable: exp(−ν kₓᵖ Δt) · exp(−ν kᵧᵖ Δt)
        Fx = jnp.exp(-nu_hyper * kx**power * dt)
        Fy = jnp.exp(-nu_hyper * ky**power * dt)
        a_f = a * Fy[:, None] * Fx[None, :]
        return a_f if spectral else self.grid.transform(a_f, inverse=True)

Functions

exponential_filter(u, alpha=36.0, power=16, spectral=False)

Separable 2D exponential filter.

F(kₓ, kᵧ) = exp(−α (kₓ/kₓ_max)ᵖ) · exp(−α (kᵧ/kᵧ_max)ᵖ)
Parameters

u : Num[Array, "Nypts Nxpts"] Physical-space field or Chebyshev coefficients. alpha : float Damping strength (≥ 0). Default 36.0. power : int Sharpening exponent (> 0). Default 16. spectral : bool If True, u is treated as spectral coefficients.

Returns

Num[Array, "Nypts Nxpts"]

Source code in spectraldiffx/_src/chebyshev/filters.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def exponential_filter(
    self,
    u: Num[Array, "Nypts Nxpts"],
    alpha: float = 36.0,
    power: int = 16,
    spectral: bool = False,
) -> Num[Array, "Nypts Nxpts"]:
    """Separable 2D exponential filter.

        F(kₓ, kᵧ) = exp(−α (kₓ/kₓ_max)ᵖ) · exp(−α (kᵧ/kᵧ_max)ᵖ)

    Parameters
    ----------
    u : Num[Array, "Nypts Nxpts"]
        Physical-space field or Chebyshev coefficients.
    alpha : float
        Damping strength (≥ 0).  Default 36.0.
    power : int
        Sharpening exponent (> 0).  Default 16.
    spectral : bool
        If ``True``, ``u`` is treated as spectral coefficients.

    Returns
    -------
    Num[Array, "Nypts Nxpts"]
    """
    if alpha < 0:
        raise ValueError(f"alpha must be >= 0, got {alpha}")
    if power <= 0:
        raise ValueError(f"power must be > 0, got {power}")
    a = u if spectral else self.grid.transform(u)
    Nx, Ny = self.grid.Nx, self.grid.Ny
    nx_modes = Nx + 1 if self.grid.node_type == "gauss-lobatto" else Nx
    ny_modes = Ny + 1 if self.grid.node_type == "gauss-lobatto" else Ny
    # Normalise by the highest mode index; derive dtype from coefficients.
    kx_max = max(1, nx_modes - 1)
    ky_max = max(1, ny_modes - 1)
    real_dtype = _coeff_dtype(a)
    kx = jnp.arange(nx_modes, dtype=real_dtype)
    ky = jnp.arange(ny_modes, dtype=real_dtype)
    Fx = jnp.exp(-alpha * (kx / kx_max) ** power)
    Fy = jnp.exp(-alpha * (ky / ky_max) ** power)
    # Separable filter: F = Fᵧ[:, None] · Fₓ[None, :]
    a_f = a * Fy[:, None] * Fx[None, :]
    return a_f if spectral else self.grid.transform(a_f, inverse=True)

hyperviscosity(u, nu_hyper, dt, power=4, spectral=False)

Separable 2D hyperviscosity filter.

F(kₓ, kᵧ) = exp(−ν_h kₓᵖ Δt) · exp(−ν_h kᵧᵖ Δt)
Parameters

u : Num[Array, "Nypts Nxpts"] Physical-space field or Chebyshev coefficients. nu_hyper : float Hyperviscosity coefficient (≥ 0). dt : float Time step for the damping (≥ 0). power : int Diffusion order (> 0). Default 4 (biharmonic). spectral : bool If True, u is treated as spectral coefficients.

Returns

Num[Array, "Nypts Nxpts"]

Source code in spectraldiffx/_src/chebyshev/filters.py 
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
def hyperviscosity(
    self,
    u: Num[Array, "Nypts Nxpts"],
    nu_hyper: float,
    dt: float,
    power: int = 4,
    spectral: bool = False,
) -> Num[Array, "Nypts Nxpts"]:
    """Separable 2D hyperviscosity filter.

        F(kₓ, kᵧ) = exp(−ν_h kₓᵖ Δt) · exp(−ν_h kᵧᵖ Δt)

    Parameters
    ----------
    u : Num[Array, "Nypts Nxpts"]
        Physical-space field or Chebyshev coefficients.
    nu_hyper : float
        Hyperviscosity coefficient (≥ 0).
    dt : float
        Time step for the damping (≥ 0).
    power : int
        Diffusion order (> 0).  Default 4 (biharmonic).
    spectral : bool
        If ``True``, ``u`` is treated as spectral coefficients.

    Returns
    -------
    Num[Array, "Nypts Nxpts"]
    """
    if nu_hyper < 0:
        raise ValueError(f"nu_hyper must be >= 0, got {nu_hyper}")
    if power <= 0:
        raise ValueError(f"power must be > 0, got {power}")
    a = u if spectral else self.grid.transform(u)
    Nx, Ny = self.grid.Nx, self.grid.Ny
    nx_modes = Nx + 1 if self.grid.node_type == "gauss-lobatto" else Nx
    ny_modes = Ny + 1 if self.grid.node_type == "gauss-lobatto" else Ny
    real_dtype = _coeff_dtype(a)
    kx = jnp.arange(nx_modes, dtype=real_dtype)
    ky = jnp.arange(ny_modes, dtype=real_dtype)
    # Separable: exp(−ν kₓᵖ Δt) · exp(−ν kᵧᵖ Δt)
    Fx = jnp.exp(-nu_hyper * kx**power * dt)
    Fy = jnp.exp(-nu_hyper * ky**power * dt)
    a_f = a * Fy[:, None] * Fx[None, :]
    return a_f if spectral else self.grid.transform(a_f, inverse=True)