Skip to content

Spherical Filters

SphericalFilter1D

Bases: Module

1D spectral filter on the Gauss–Legendre latitude grid.

Filters multiply each Legendre coefficient by a kernel F(l):

c̃ₗ = F(l) · cₗ

Attributes

grid : SphericalGrid1D Underlying 1D Gauss–Legendre grid.

Examples

import jax.numpy as jnp grid = SphericalGrid1D.from_N_L(N=64, L=jnp.pi) flt = SphericalFilter1D(grid=grid) u = jnp.cos(grid.x) + 1e-3 * jnp.cos(60 * grid.x) u_smooth = flt.exponential_filter(u)

Source code in spectraldiffx/_src/spherical/filters.py 
 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
class SphericalFilter1D(eqx.Module):
    """1D spectral filter on the Gauss–Legendre latitude grid.

    Filters multiply each Legendre coefficient by a kernel F(l):

        c̃ₗ = F(l) · cₗ

    Attributes
    ----------
    grid : SphericalGrid1D
        Underlying 1D Gauss–Legendre grid.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = SphericalGrid1D.from_N_L(N=64, L=jnp.pi)
    >>> flt = SphericalFilter1D(grid=grid)
    >>> u = jnp.cos(grid.x) + 1e-3 * jnp.cos(60 * grid.x)
    >>> u_smooth = flt.exponential_filter(u)
    """

    grid: SphericalGrid1D

    def exponential_filter(
        self,
        u: Num[Array, "N"],
        alpha: float = 36.0,
        power: int = 16,
        spectral: bool = False,
    ) -> Num[Array, "N"]:
        """Exponential filter in Legendre coefficient space:

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

        Near unity for low l, falls sharply near lₘₐₓ = N−1.

        Parameters
        ----------
        u : Num[Array, "N"]
            Physical field or Legendre coefficients (if ``spectral=True``).
        alpha : float
            Damping coefficient (≥ 0).  Default 36.0 (≈ ε_64 damping at lₘₐₓ).
        power : int
            Filter sharpness (> 0, even).  Default 16.
        spectral : bool
            If ``True``, treat ``u`` as Legendre coefficients.

        Returns
        -------
        Num[Array, "N"]
        """
        if alpha < 0:
            raise ValueError(f"alpha must be >= 0, got {alpha}")
        if power <= 0:
            raise ValueError(f"power must be > 0, got {power}")
        c = u if spectral else self.grid.transform(u)
        l = self.grid.l
        l_max = float(self.grid.N - 1)
        l_max_safe = jnp.where(l_max == 0, 1.0, l_max)
        mask = jnp.exp(-alpha * (l / l_max_safe) ** power)
        c_f = c * mask
        return c_f if spectral else self.grid.transform(c_f, inverse=True)

    def hyperviscosity(
        self,
        u: Num[Array, "N"],
        nu_hyper: float,
        dt: float,
        power: int = 4,
        spectral: bool = False,
    ) -> Num[Array, "N"]:
        """Hyperviscous damping driven by the Laplace–Beltrami eigenvalue:

            F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

        where R = ``grid.L / π`` is the sphere radius.  Simulates
        high-order diffusion ``∂u/∂t = (−1)^{p+1} ν_h ∇^p u`` with the
        damping rate independent of R for fixed ν_h.

        Parameters
        ----------
        u : Num[Array, "N"]
            Physical field or Legendre coefficients.
        nu_hyper : float
            Hyperviscosity coefficient (≥ 0).
        dt : float
            Time step (≥ 0).
        power : int
            Laplacian power p (> 0, even).  Default 4 (biharmonic).
        spectral : bool
            If ``True``, treat ``u`` as Legendre coefficients.

        Returns
        -------
        Num[Array, "N"]
        """
        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}")
        c = u if spectral else self.grid.transform(u)
        R = self.grid.L / jnp.pi
        l = self.grid.l
        eigenval = l * (l + 1) / (R**2)
        mask = jnp.exp(-nu_hyper * eigenval ** (power / 2.0) * dt)
        c_f = c * mask
        return c_f if spectral else self.grid.transform(c_f, inverse=True)

Functions

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

Exponential filter in Legendre coefficient space:

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

Near unity for low l, falls sharply near lₘₐₓ = N−1.

Parameters

u : Num[Array, "N"] Physical field or Legendre coefficients (if spectral=True). alpha : float Damping coefficient (≥ 0). Default 36.0 (≈ ε_64 damping at lₘₐₓ). power : int Filter sharpness (> 0, even). Default 16. spectral : bool If True, treat u as Legendre coefficients.

Returns

Num[Array, "N"]

Source code in spectraldiffx/_src/spherical/filters.py 
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
def exponential_filter(
    self,
    u: Num[Array, "N"],
    alpha: float = 36.0,
    power: int = 16,
    spectral: bool = False,
) -> Num[Array, "N"]:
    """Exponential filter in Legendre coefficient space:

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

    Near unity for low l, falls sharply near lₘₐₓ = N−1.

    Parameters
    ----------
    u : Num[Array, "N"]
        Physical field or Legendre coefficients (if ``spectral=True``).
    alpha : float
        Damping coefficient (≥ 0).  Default 36.0 (≈ ε_64 damping at lₘₐₓ).
    power : int
        Filter sharpness (> 0, even).  Default 16.
    spectral : bool
        If ``True``, treat ``u`` as Legendre coefficients.

    Returns
    -------
    Num[Array, "N"]
    """
    if alpha < 0:
        raise ValueError(f"alpha must be >= 0, got {alpha}")
    if power <= 0:
        raise ValueError(f"power must be > 0, got {power}")
    c = u if spectral else self.grid.transform(u)
    l = self.grid.l
    l_max = float(self.grid.N - 1)
    l_max_safe = jnp.where(l_max == 0, 1.0, l_max)
    mask = jnp.exp(-alpha * (l / l_max_safe) ** power)
    c_f = c * mask
    return c_f if spectral else self.grid.transform(c_f, inverse=True)

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

Hyperviscous damping driven by the Laplace–Beltrami eigenvalue:

F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

where R = grid.L / π is the sphere radius. Simulates high-order diffusion ∂u/∂t = (−1)^{p+1} ν_h ∇^p u with the damping rate independent of R for fixed ν_h.

Parameters

u : Num[Array, "N"] Physical field or Legendre coefficients. nu_hyper : float Hyperviscosity coefficient (≥ 0). dt : float Time step (≥ 0). power : int Laplacian power p (> 0, even). Default 4 (biharmonic). spectral : bool If True, treat u as Legendre coefficients.

Returns

Num[Array, "N"]

Source code in spectraldiffx/_src/spherical/filters.py 
 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
def hyperviscosity(
    self,
    u: Num[Array, "N"],
    nu_hyper: float,
    dt: float,
    power: int = 4,
    spectral: bool = False,
) -> Num[Array, "N"]:
    """Hyperviscous damping driven by the Laplace–Beltrami eigenvalue:

        F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

    where R = ``grid.L / π`` is the sphere radius.  Simulates
    high-order diffusion ``∂u/∂t = (−1)^{p+1} ν_h ∇^p u`` with the
    damping rate independent of R for fixed ν_h.

    Parameters
    ----------
    u : Num[Array, "N"]
        Physical field or Legendre coefficients.
    nu_hyper : float
        Hyperviscosity coefficient (≥ 0).
    dt : float
        Time step (≥ 0).
    power : int
        Laplacian power p (> 0, even).  Default 4 (biharmonic).
    spectral : bool
        If ``True``, treat ``u`` as Legendre coefficients.

    Returns
    -------
    Num[Array, "N"]
    """
    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}")
    c = u if spectral else self.grid.transform(u)
    R = self.grid.L / jnp.pi
    l = self.grid.l
    eigenval = l * (l + 1) / (R**2)
    mask = jnp.exp(-nu_hyper * eigenval ** (power / 2.0) * dt)
    c_f = c * mask
    return c_f if spectral else self.grid.transform(c_f, inverse=True)

SphericalFilter2D

Bases: Module

2D spectral filter on the full sphere lat-lon grid.

Applies multiplicative kernels F(l) in spherical-harmonic-coefficient space (broadcast across all m for each l).

Attributes

grid : SphericalGrid2D Underlying 2D lat-lon grid.

Examples

import jax.numpy as jnp grid = SphericalGrid2D.from_N_L(Nx=64, Ny=32) flt = SphericalFilter2D(grid=grid) PHI, THETA = grid.X u = jnp.sin(THETA) * jnp.cos(4 * PHI) u_smooth = flt.exponential_filter(u, alpha=16.0, power=8)

Source code in spectraldiffx/_src/spherical/filters.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
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
class SphericalFilter2D(eqx.Module):
    """2D spectral filter on the full sphere lat-lon grid.

    Applies multiplicative kernels F(l) in spherical-harmonic-coefficient
    space (broadcast across all m for each l).

    Attributes
    ----------
    grid : SphericalGrid2D
        Underlying 2D lat-lon grid.

    Examples
    --------
    >>> import jax.numpy as jnp
    >>> grid = SphericalGrid2D.from_N_L(Nx=64, Ny=32)
    >>> flt = SphericalFilter2D(grid=grid)
    >>> PHI, THETA = grid.X
    >>> u = jnp.sin(THETA) * jnp.cos(4 * PHI)
    >>> u_smooth = flt.exponential_filter(u, alpha=16.0, power=8)
    """

    grid: SphericalGrid2D

    def exponential_filter(
        self,
        u: Num[Array, "Nlat Nlon"],
        alpha: float = 36.0,
        power: int = 16,
        spectral: bool = False,
    ) -> Num[Array, "Nlat Nlon"]:
        """Exponential filter in (l, m) space:

            F(l, m) = exp(−α · (l / lₘₐₓ)ᵖ)

        Parameters
        ----------
        u : Num[Array, "Nlat Nlon"]
            Physical field or SHT coefficients.
        alpha : float
            Damping coefficient (≥ 0).
        power : int
            Filter sharpness (> 0).
        spectral : bool
            If ``True``, treat ``u`` as SHT 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}")
        u_hat = u if spectral else self.grid.transform(u)
        l = self.grid.l  # (Nl,)
        l_max = float(self.grid.Ny - 1)
        l_max_safe = jnp.where(l_max == 0.0, 1.0, l_max)
        mask = jnp.exp(-alpha * (l / l_max_safe) ** power)  # (Nl,)
        u_hat_f = u_hat * mask[:, None]  # broadcast over m
        return u_hat_f if spectral else self.grid.transform(u_hat_f, inverse=True)

    def hyperviscosity(
        self,
        u: Num[Array, "Nlat Nlon"],
        nu_hyper: float,
        dt: float,
        power: int = 4,
        spectral: bool = False,
    ) -> Num[Array, "Nlat Nlon"]:
        """Hyperviscous damping using the Laplace–Beltrami eigenvalue:

            F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

        where R = ``grid.Ly / π``.  Applied identically to every m at a
        given l (isotropic on the sphere).

        Parameters
        ----------
        u : Num[Array, "Nlat Nlon"]
            Physical field or SHT coefficients.
        nu_hyper : float
            Hyperviscosity coefficient (≥ 0).
        dt : float
            Time step (≥ 0).
        power : int
            Laplacian power p (> 0).  Default 4 (biharmonic damping).
        spectral : bool
            If ``True``, treat ``u`` as SHT 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}")
        u_hat = u if spectral else self.grid.transform(u)
        R = self.grid.Ly / jnp.pi
        l = self.grid.l  # (Nl,)
        eigenval = l * (l + 1) / (R**2)
        mask = jnp.exp(-nu_hyper * eigenval ** (power / 2.0) * dt)  # (Nl,)
        u_hat_f = u_hat * mask[:, None]
        return u_hat_f if spectral else self.grid.transform(u_hat_f, inverse=True)

Functions

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

Exponential filter in (l, m) space:

F(l, m) = exp(−α · (l / lₘₐₓ)ᵖ)
Parameters

u : Num[Array, "Nlat Nlon"] Physical field or SHT coefficients. alpha : float Damping coefficient (≥ 0). power : int Filter sharpness (> 0). spectral : bool If True, treat u as SHT coefficients.

Source code in spectraldiffx/_src/spherical/filters.py 
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
def exponential_filter(
    self,
    u: Num[Array, "Nlat Nlon"],
    alpha: float = 36.0,
    power: int = 16,
    spectral: bool = False,
) -> Num[Array, "Nlat Nlon"]:
    """Exponential filter in (l, m) space:

        F(l, m) = exp(−α · (l / lₘₐₓ)ᵖ)

    Parameters
    ----------
    u : Num[Array, "Nlat Nlon"]
        Physical field or SHT coefficients.
    alpha : float
        Damping coefficient (≥ 0).
    power : int
        Filter sharpness (> 0).
    spectral : bool
        If ``True``, treat ``u`` as SHT 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}")
    u_hat = u if spectral else self.grid.transform(u)
    l = self.grid.l  # (Nl,)
    l_max = float(self.grid.Ny - 1)
    l_max_safe = jnp.where(l_max == 0.0, 1.0, l_max)
    mask = jnp.exp(-alpha * (l / l_max_safe) ** power)  # (Nl,)
    u_hat_f = u_hat * mask[:, None]  # broadcast over m
    return u_hat_f if spectral else self.grid.transform(u_hat_f, inverse=True)

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

Hyperviscous damping using the Laplace–Beltrami eigenvalue:

F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

where R = grid.Ly / π. Applied identically to every m at a given l (isotropic on the sphere).

Parameters

u : Num[Array, "Nlat Nlon"] Physical field or SHT coefficients. nu_hyper : float Hyperviscosity coefficient (≥ 0). dt : float Time step (≥ 0). power : int Laplacian power p (> 0). Default 4 (biharmonic damping). spectral : bool If True, treat u as SHT coefficients.

Source code in spectraldiffx/_src/spherical/filters.py 
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
def hyperviscosity(
    self,
    u: Num[Array, "Nlat Nlon"],
    nu_hyper: float,
    dt: float,
    power: int = 4,
    spectral: bool = False,
) -> Num[Array, "Nlat Nlon"]:
    """Hyperviscous damping using the Laplace–Beltrami eigenvalue:

        F(l) = exp(−ν_h · [l(l+1)/R²]^(p/2) · Δt)

    where R = ``grid.Ly / π``.  Applied identically to every m at a
    given l (isotropic on the sphere).

    Parameters
    ----------
    u : Num[Array, "Nlat Nlon"]
        Physical field or SHT coefficients.
    nu_hyper : float
        Hyperviscosity coefficient (≥ 0).
    dt : float
        Time step (≥ 0).
    power : int
        Laplacian power p (> 0).  Default 4 (biharmonic damping).
    spectral : bool
        If ``True``, treat ``u`` as SHT 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}")
    u_hat = u if spectral else self.grid.transform(u)
    R = self.grid.Ly / jnp.pi
    l = self.grid.l  # (Nl,)
    eigenval = l * (l + 1) / (R**2)
    mask = jnp.exp(-nu_hyper * eigenval ** (power / 2.0) * dt)  # (Nl,)
    u_hat_f = u_hat * mask[:, None]
    return u_hat_f if spectral else self.grid.transform(u_hat_f, inverse=True)