Skip to content

Navier-Stokes¤

2D: Streamfunction-Vorticity Formulation¤

In 2D, the vorticity is a scalar, so the streamfunction-vorticity formulation reduces the system to a single-channel PDE — making it the natural choice for 2D spectral methods.

exponax.stepper.NavierStokesVorticity ¤

Bases: BaseStepper

Source code in exponax/stepper/_navier_stokes.py 
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 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
152
153
class NavierStokesVorticity(BaseStepper):
    diffusivity: float
    vorticity_convection_scale: float
    drag: float
    dealiasing_fraction: float

    def __init__(
        self,
        num_spatial_dims: int,
        domain_extent: float,
        num_points: int,
        dt: float,
        *,
        diffusivity: float = 0.01,
        vorticity_convection_scale: float = 1.0,
        drag: float = 0.0,
        order: int = 2,
        dealiasing_fraction: float = 2 / 3,
        num_circle_points: int = 16,
        circle_radius: float = 1.0,
    ):
        """
        Timestepper for the 2d Navier-Stokes equation on periodic boundary
        conditions in streamfunction-vorticity formulation. The equation reads

        ```
            uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu
        ```

        with `u` the vorticity. On the right-hand side the first term is a drag
        with coefficient `λ` and the second term is a diffusion with coefficient
        `ν`. The operation on the left-hand-side `([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u` is
        the "vorticity" convection which is scaled by `b`. It consists of the
        solution to the Poisson problem via the inverse Laplacian `Δ⁻¹` and the
        gradient `∇` of the streamfunction. The term `[1, -1]ᵀ ⊙` negates the
        second component of the gradient.

        We can map the vorticity to a (two-channel) velocity field by `∇
        (Δ⁻¹u)`.

        The expected temporal behavior is that the initial vorticity field
        continues to swirl but decays over time.

        Let `U = ‖∇ (Δ⁻¹u)‖` denote the magnitude of the velocity field, then
        the Reynolds number of the problem is `Re = U L / ν` with `L` the
        `domain_extent`.

        **Arguments:**

        - `num_spatial_dims`: The number of spatial dimensions `d`.
        - `domain_extent`: The size of the domain `L`; in higher dimensions
            the domain is assumed to be a scaled hypercube `Ω = (0, L)ᵈ`.
        - `num_points`: The number of points `N` used to discretize the
            domain. This **includes** the left boundary point and **excludes**
            the right boundary point. In higher dimensions; the number of points
            in each dimension is the same. Hence, the total number of degrees of
            freedom is `Nᵈ`.
        - `dt`: The timestep size `Δt` between two consecutive states.
        - `diffusivity`: The diffusivity coefficient `ν`. This affects the
            Reynolds number. The lower the diffusivity, the "more turbulent".
            Default is `0.01`.
        - `vorticity_convection_scale`: The scaling factor for the vorticity
            convection term. Default is `1.0`.
        - `drag`: The drag coefficient `λ`. Default is `0.0`.
        - `order`: The order of the Exponential Time Differencing Runge
            Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
            solves the linear part of the equation. Use higher values for higher
            accuracy and stability. The default choice of `2` is a good
            compromise for single precision floats.
        - `dealiasing_fraction`: The fraction of the wavenumbers to keep
            before evaluating the nonlinearity. The default 2/3 corresponds to
            Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
            2/3.
        - `num_circle_points`: How many points to use in the complex contour
            integral method to compute the coefficients of the exponential time
            differencing Runge Kutta method. Default: 16.
        - `circle_radius`: The radius of the contour used to compute the
            coefficients of the exponential time differencing Runge Kutta
            method. Default: 1.0.

        **Notes:**

        - The Reynolds number is a measure of whether the problem is dominated
            by diffusive or convective effects. The higher the Reynolds number,
            the stronger the effect of the convective. Since this term is the
            nonlinear one, the higher the Reynolds number, the worse the ETDRK
            methods become in comparison to other approaches. That is because
            those methods are better for semi-linear PDEs in which the difficult
            part is the linear one.
        - The higher the Reynolds number, the smaller the timestep size must
            be to ensure stability.

        **Good Values:**

        - `domain_extent = 1`, `num_points=50`, `dt=0.01`,
            `diffusivity=0.0003`, together with an initial condition in which
            only the first few wavenumbers are excited gives a nice decaying
            turbulence demo.
        - Use the repeated stepper to perform 10 substeps to have faster
            dynamics.
        """
        if num_spatial_dims != 2:
            raise ValueError(
                f"Expected num_spatial_dims = 2, got {num_spatial_dims}. "
                "For 3D, use NavierStokesVelocity instead."
            )

        self.diffusivity = diffusivity
        self.vorticity_convection_scale = vorticity_convection_scale
        self.drag = drag
        self.dealiasing_fraction = dealiasing_fraction
        super().__init__(
            num_spatial_dims=num_spatial_dims,
            domain_extent=domain_extent,
            num_points=num_points,
            dt=dt,
            num_channels=1,
            order=order,
            num_circle_points=num_circle_points,
            circle_radius=circle_radius,
        )

    def _build_linear_operator(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> Complex[Array, "1 ... (N//2)+1"]:
        return self.diffusivity * build_laplace_operator(
            derivative_operator, order=2
        ) + self.drag * build_laplace_operator(derivative_operator, order=0)

    def _build_nonlinear_fun(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> VorticityConvection2d:
        return VorticityConvection2d(
            self.num_spatial_dims,
            self.num_points,
            convection_scale=self.vorticity_convection_scale,
            derivative_operator=derivative_operator,
            dealiasing_fraction=self.dealiasing_fraction,
        )
__init__ ¤
__init__(
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    vorticity_convection_scale: float = 1.0,
    drag: float = 0.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0
)

Timestepper for the 2d Navier-Stokes equation on periodic boundary conditions in streamfunction-vorticity formulation. The equation reads

    uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu

with u the vorticity. On the right-hand side the first term is a drag with coefficient λ and the second term is a diffusion with coefficient ν. The operation on the left-hand-side ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u is the "vorticity" convection which is scaled by b. It consists of the solution to the Poisson problem via the inverse Laplacian Δ⁻¹ and the gradient ∇ of the streamfunction. The term [1, -1]ᵀ ⊙ negates the second component of the gradient.

We can map the vorticity to a (two-channel) velocity field by ∇ (Δ⁻¹u).

The expected temporal behavior is that the initial vorticity field continues to swirl but decays over time.

Let U = ‖∇ (Δ⁻¹u)‖ denote the magnitude of the velocity field, then the Reynolds number of the problem is Re = U L / ν with L the domain_extent.

Arguments:

  • num_spatial_dims: The number of spatial dimensions d.
  • domain_extent: The size of the domain L; in higher dimensions the domain is assumed to be a scaled hypercube Ω = (0, L)ᵈ.
  • num_points: The number of points N used to discretize the domain. This includes the left boundary point and excludes the right boundary point. In higher dimensions; the number of points in each dimension is the same. Hence, the total number of degrees of freedom is Nᵈ.
  • dt: The timestep size Δt between two consecutive states.
  • diffusivity: The diffusivity coefficient ν. This affects the Reynolds number. The lower the diffusivity, the "more turbulent". Default is 0.01.
  • vorticity_convection_scale: The scaling factor for the vorticity convection term. Default is 1.0.
  • drag: The drag coefficient λ. Default is 0.0.
  • order: The order of the Exponential Time Differencing Runge Kutta method. Must be one of {0, 1, 2, 3, 4}. The option 0 only solves the linear part of the equation. Use higher values for higher accuracy and stability. The default choice of 2 is a good compromise for single precision floats.
  • dealiasing_fraction: The fraction of the wavenumbers to keep before evaluating the nonlinearity. The default 2/3 corresponds to Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default: 2/3.
  • num_circle_points: How many points to use in the complex contour integral method to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 16.
  • circle_radius: The radius of the contour used to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 1.0.

Notes:

  • The Reynolds number is a measure of whether the problem is dominated by diffusive or convective effects. The higher the Reynolds number, the stronger the effect of the convective. Since this term is the nonlinear one, the higher the Reynolds number, the worse the ETDRK methods become in comparison to other approaches. That is because those methods are better for semi-linear PDEs in which the difficult part is the linear one.
  • The higher the Reynolds number, the smaller the timestep size must be to ensure stability.

Good Values:

  • domain_extent = 1, num_points=50, dt=0.01, diffusivity=0.0003, together with an initial condition in which only the first few wavenumbers are excited gives a nice decaying turbulence demo.
  • Use the repeated stepper to perform 10 substeps to have faster dynamics.
Source code in exponax/stepper/_navier_stokes.py 
 19
 20
 21
 22
 23
 24
 25
 26
 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
def __init__(
    self,
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    vorticity_convection_scale: float = 1.0,
    drag: float = 0.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0,
):
    """
    Timestepper for the 2d Navier-Stokes equation on periodic boundary
    conditions in streamfunction-vorticity formulation. The equation reads

    ```
        uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu
    ```

    with `u` the vorticity. On the right-hand side the first term is a drag
    with coefficient `λ` and the second term is a diffusion with coefficient
    `ν`. The operation on the left-hand-side `([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u` is
    the "vorticity" convection which is scaled by `b`. It consists of the
    solution to the Poisson problem via the inverse Laplacian `Δ⁻¹` and the
    gradient `∇` of the streamfunction. The term `[1, -1]ᵀ ⊙` negates the
    second component of the gradient.

    We can map the vorticity to a (two-channel) velocity field by `∇
    (Δ⁻¹u)`.

    The expected temporal behavior is that the initial vorticity field
    continues to swirl but decays over time.

    Let `U = ‖∇ (Δ⁻¹u)‖` denote the magnitude of the velocity field, then
    the Reynolds number of the problem is `Re = U L / ν` with `L` the
    `domain_extent`.

    **Arguments:**

    - `num_spatial_dims`: The number of spatial dimensions `d`.
    - `domain_extent`: The size of the domain `L`; in higher dimensions
        the domain is assumed to be a scaled hypercube `Ω = (0, L)ᵈ`.
    - `num_points`: The number of points `N` used to discretize the
        domain. This **includes** the left boundary point and **excludes**
        the right boundary point. In higher dimensions; the number of points
        in each dimension is the same. Hence, the total number of degrees of
        freedom is `Nᵈ`.
    - `dt`: The timestep size `Δt` between two consecutive states.
    - `diffusivity`: The diffusivity coefficient `ν`. This affects the
        Reynolds number. The lower the diffusivity, the "more turbulent".
        Default is `0.01`.
    - `vorticity_convection_scale`: The scaling factor for the vorticity
        convection term. Default is `1.0`.
    - `drag`: The drag coefficient `λ`. Default is `0.0`.
    - `order`: The order of the Exponential Time Differencing Runge
        Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
        solves the linear part of the equation. Use higher values for higher
        accuracy and stability. The default choice of `2` is a good
        compromise for single precision floats.
    - `dealiasing_fraction`: The fraction of the wavenumbers to keep
        before evaluating the nonlinearity. The default 2/3 corresponds to
        Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
        2/3.
    - `num_circle_points`: How many points to use in the complex contour
        integral method to compute the coefficients of the exponential time
        differencing Runge Kutta method. Default: 16.
    - `circle_radius`: The radius of the contour used to compute the
        coefficients of the exponential time differencing Runge Kutta
        method. Default: 1.0.

    **Notes:**

    - The Reynolds number is a measure of whether the problem is dominated
        by diffusive or convective effects. The higher the Reynolds number,
        the stronger the effect of the convective. Since this term is the
        nonlinear one, the higher the Reynolds number, the worse the ETDRK
        methods become in comparison to other approaches. That is because
        those methods are better for semi-linear PDEs in which the difficult
        part is the linear one.
    - The higher the Reynolds number, the smaller the timestep size must
        be to ensure stability.

    **Good Values:**

    - `domain_extent = 1`, `num_points=50`, `dt=0.01`,
        `diffusivity=0.0003`, together with an initial condition in which
        only the first few wavenumbers are excited gives a nice decaying
        turbulence demo.
    - Use the repeated stepper to perform 10 substeps to have faster
        dynamics.
    """
    if num_spatial_dims != 2:
        raise ValueError(
            f"Expected num_spatial_dims = 2, got {num_spatial_dims}. "
            "For 3D, use NavierStokesVelocity instead."
        )

    self.diffusivity = diffusivity
    self.vorticity_convection_scale = vorticity_convection_scale
    self.drag = drag
    self.dealiasing_fraction = dealiasing_fraction
    super().__init__(
        num_spatial_dims=num_spatial_dims,
        domain_extent=domain_extent,
        num_points=num_points,
        dt=dt,
        num_channels=1,
        order=order,
        num_circle_points=num_circle_points,
        circle_radius=circle_radius,
    )
__call__ ¤
__call__(
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]

Perform one step of the time integration for a single state.

Arguments:

  • u: The state vector, shape (C, ..., N,).

Returns:

  • u_next: The state vector after one step, shape (C, ..., N,).

Tip

Use this call method together with exponax.rollout to efficiently produce temporal trajectories.

Info

For batched operation, use jax.vmap on this function.

Source code in exponax/_base_stepper.py 
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
269
270
271
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Perform one step of the time integration for a single state.

    **Arguments:**

    - `u`: The state vector, shape `(C, ..., N,)`.

    **Returns:**

    - `u_next`: The state vector after one step, shape `(C, ..., N,)`.

    !!! tip
        Use this call method together with `exponax.rollout` to efficiently
        produce temporal trajectories.

    !!! info
        For batched operation, use `jax.vmap` on this function.
    """
    expected_shape = (self.num_channels,) + spatial_shape(
        self.num_spatial_dims, self.num_points
    )
    if u.shape != expected_shape:
        raise ValueError(
            f"""Expected shape {expected_shape}, got {u.shape}. For batched
             operation use `jax.vmap` on this function."""
        )
    return self.step(u)

exponax.stepper.KolmogorovFlowVorticity ¤

Bases: BaseStepper

Source code in exponax/stepper/_navier_stokes.py 
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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
class KolmogorovFlowVorticity(BaseStepper):
    diffusivity: float
    convection_scale: float
    drag: float
    injection_mode: int
    injection_scale: float
    dealiasing_fraction: float

    def __init__(
        self,
        num_spatial_dims: int,
        domain_extent: float,
        num_points: int,
        dt: float,
        *,
        diffusivity: float = 0.001,
        convection_scale: float = 1.0,
        drag: float = -0.1,
        injection_mode: int = 4,
        injection_scale: float = 1.0,
        order: int = 2,
        dealiasing_fraction: float = 2 / 3,
        num_circle_points: int = 16,
        circle_radius: float = 1.0,
    ):
        """
        Timestepper for the 2d Kolmogorov flow equation on periodic boundary
        conditions in streamfunction-vorticity formulation. The equation reads

        ```
            uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu + f
        ```

        For a detailed description of the terms, see the documentation of the
        `NavierStokesVorticity` stepper. The only difference is the additional
        forcing term `f` which injects new energy into the system. For a
        Kolmogorov flow **in primary variables**, it has the form

        ```
            f₀ = γ sin(k (2π/L) x₁)

            f₁ = 0
        ```

        In words, only the first channel is forced at a specific wavenumber over
        the second axis. Since this stepper considers the
        streamfunction-vorticity formulation, we take its curl to get

        ```
            f = -k (2π/L) γ cos(k (2π/L) x₁)
        ```

        The expected temporal behavior is that the initial vorticity field first
        is excited into a noisy striped pattern. This pattern breaks up and a
        turbulent spatio-temporal chaos emerges.

        A negative drag coefficient `λ` is needed to remove some of the energy
        piling up in low modes.

        According to

            Chandler, G.J. and Kerswell, R.R. (2013) ‘Invariant recurrent
            solutions embedded in a turbulent two-dimensional Kolmogorov flow’,
            Journal of Fluid Mechanics, 722, pp. 554–595.
            doi:10.1017/jfm.2013.122.

        equation (2.5), the Reynolds number of the Kolmogorov flow is given by

            Re = √ζ / ν √(L / (2π))³

        with `ζ` being the scaling of the Kolmogorov forcing, i.e., the
        `injection_scale`. Hence, in the case of `L = 2π`, `ζ = 1`, the Reynolds
        number is `Re = 1 / ν`. If one uses the default value of `ν = 0.001`,
        the Reynolds number is `Re = 1000` which also corresponds to the main
        experiments in

            Kochkov, D., Smith, J.A., Alieva, A., Wang, Q., Brenner, M.P. and
            Hoyer, S., 2021. Machine learning–accelerated computational fluid
            dynamics. Proceedings of the National Academy of Sciences, 118(21),
            p.e2101784118.

        together with `injection_mode = 4`. Note that they required a resolution
        of `num_points = 2048` (=> 2048^2 = 4.2M degrees of freedom in 2d) to
        fully resolve all scales at that Reynolds number. Using `Re = 0.01`
        which corresponds to `ν = 0.01` can be a good starting for
        `num_points=128`.


        **Arguments:**

        - `num_spatial_dims`: The number of spatial dimensions `d`.
        - `domain_extent`: The size of the domain `L`; in higher dimensions
            the domain is assumed to be a scaled hypercube `Ω = (0, L)ᵈ`.
        - `num_points`: The number of points `N` used to discretize the
            domain. This **includes** the left boundary point and **excludes**
            the right boundary point. In higher dimensions; the number of points
            in each dimension is the same. Hence, the total number of degrees of
            freedom is `Nᵈ`.
        - `dt`: The timestep size `Δt` between two consecutive states.
        - `diffusivity`: The diffusivity coefficient `ν`. This affects the
            Reynolds number. The lower the diffusivity, the "more turbulent".
            Default is `0.001`.
        - `convection_scale`: The scaling factor for the vorticity
            convection term. Default is `1.0`.
        - `drag`: The drag coefficient `λ`. Default is `-0.1`.
        - `injection_mode`: The mode of the injection. Default is `4`.
        - `injection_scale`: The scaling factor for the injection. Default
            is `1.0`.
        - `order`: The order of the Exponential Time Differencing Runge
            Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
            solves the linear part of the equation. Use higher values for higher
            accuracy and stability. The default choice of `2` is a good
            compromise for single precision floats.
        - `dealiasing_fraction`: The fraction of the wavenumbers to keep
            before evaluating the nonlinearity. The default 2/3 corresponds to
            Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
            2/3.
        - `num_circle_points`: How many points to use in the complex contour
            integral method to compute the coefficients of the exponential time
            differencing Runge Kutta method. Default: 16.
        - `circle_radius`: The radius of the contour used to compute the
            coefficients of the exponential time differencing Runge Kutta
            method. Default: 1.0.

        **Notes:**

        - Interesting pointer to explore:
          http://trieste-conf.itp.ac.ru/Boffetta.pdf
        """
        if num_spatial_dims != 2:
            raise ValueError(
                f"Expected num_spatial_dims = 2, got {num_spatial_dims}. "
                "For 3D, use KolmogorovFlowVelocity instead."
            )
        self.diffusivity = diffusivity
        self.convection_scale = convection_scale
        self.drag = drag
        self.injection_mode = injection_mode
        self.injection_scale = injection_scale
        self.dealiasing_fraction = dealiasing_fraction
        super().__init__(
            num_spatial_dims=num_spatial_dims,
            domain_extent=domain_extent,
            num_points=num_points,
            dt=dt,
            num_channels=1,
            order=order,
            num_circle_points=num_circle_points,
            circle_radius=circle_radius,
        )

    def _build_linear_operator(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> Complex[Array, "1 ... (N//2)+1"]:
        return self.diffusivity * build_laplace_operator(
            derivative_operator, order=2
        ) + self.drag * build_laplace_operator(derivative_operator, order=0)

    def _build_nonlinear_fun(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> VorticityConvection2dKolmogorov:
        return VorticityConvection2dKolmogorov(
            self.num_spatial_dims,
            self.num_points,
            convection_scale=self.convection_scale,
            injection_mode=self.injection_mode,
            injection_scale=self.injection_scale,
            derivative_operator=derivative_operator,
            dealiasing_fraction=self.dealiasing_fraction,
        )
__init__ ¤
__init__(
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.001,
    convection_scale: float = 1.0,
    drag: float = -0.1,
    injection_mode: int = 4,
    injection_scale: float = 1.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0
)

Timestepper for the 2d Kolmogorov flow equation on periodic boundary conditions in streamfunction-vorticity formulation. The equation reads

    uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu + f

For a detailed description of the terms, see the documentation of the NavierStokesVorticity stepper. The only difference is the additional forcing term f which injects new energy into the system. For a Kolmogorov flow in primary variables, it has the form

    f₀ = γ sin(k (2π/L) x₁)

    f₁ = 0

In words, only the first channel is forced at a specific wavenumber over the second axis. Since this stepper considers the streamfunction-vorticity formulation, we take its curl to get

    f = -k (2π/L) γ cos(k (2π/L) x₁)

The expected temporal behavior is that the initial vorticity field first is excited into a noisy striped pattern. This pattern breaks up and a turbulent spatio-temporal chaos emerges.

A negative drag coefficient λ is needed to remove some of the energy piling up in low modes.

According to

Chandler, G.J. and Kerswell, R.R. (2013) ‘Invariant recurrent
solutions embedded in a turbulent two-dimensional Kolmogorov flow’,
Journal of Fluid Mechanics, 722, pp. 554–595.
doi:10.1017/jfm.2013.122.

equation (2.5), the Reynolds number of the Kolmogorov flow is given by

Re = √ζ / ν √(L / (2π))³

with ζ being the scaling of the Kolmogorov forcing, i.e., the injection_scale. Hence, in the case of L = 2π, ζ = 1, the Reynolds number is Re = 1 / ν. If one uses the default value of ν = 0.001, the Reynolds number is Re = 1000 which also corresponds to the main experiments in

Kochkov, D., Smith, J.A., Alieva, A., Wang, Q., Brenner, M.P. and
Hoyer, S., 2021. Machine learning–accelerated computational fluid
dynamics. Proceedings of the National Academy of Sciences, 118(21),
p.e2101784118.

together with injection_mode = 4. Note that they required a resolution of num_points = 2048 (=> 2048^2 = 4.2M degrees of freedom in 2d) to fully resolve all scales at that Reynolds number. Using Re = 0.01 which corresponds to ν = 0.01 can be a good starting for num_points=128.

Arguments:

  • num_spatial_dims: The number of spatial dimensions d.
  • domain_extent: The size of the domain L; in higher dimensions the domain is assumed to be a scaled hypercube Ω = (0, L)ᵈ.
  • num_points: The number of points N used to discretize the domain. This includes the left boundary point and excludes the right boundary point. In higher dimensions; the number of points in each dimension is the same. Hence, the total number of degrees of freedom is Nᵈ.
  • dt: The timestep size Δt between two consecutive states.
  • diffusivity: The diffusivity coefficient ν. This affects the Reynolds number. The lower the diffusivity, the "more turbulent". Default is 0.001.
  • convection_scale: The scaling factor for the vorticity convection term. Default is 1.0.
  • drag: The drag coefficient λ. Default is -0.1.
  • injection_mode: The mode of the injection. Default is 4.
  • injection_scale: The scaling factor for the injection. Default is 1.0.
  • order: The order of the Exponential Time Differencing Runge Kutta method. Must be one of {0, 1, 2, 3, 4}. The option 0 only solves the linear part of the equation. Use higher values for higher accuracy and stability. The default choice of 2 is a good compromise for single precision floats.
  • dealiasing_fraction: The fraction of the wavenumbers to keep before evaluating the nonlinearity. The default 2/3 corresponds to Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default: 2/3.
  • num_circle_points: How many points to use in the complex contour integral method to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 16.
  • circle_radius: The radius of the contour used to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 1.0.

Notes:

  • Interesting pointer to explore: http://trieste-conf.itp.ac.ru/Boffetta.pdf
Source code in exponax/stepper/_navier_stokes.py 
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
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
def __init__(
    self,
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.001,
    convection_scale: float = 1.0,
    drag: float = -0.1,
    injection_mode: int = 4,
    injection_scale: float = 1.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0,
):
    """
    Timestepper for the 2d Kolmogorov flow equation on periodic boundary
    conditions in streamfunction-vorticity formulation. The equation reads

    ```
        uₜ + b ([1, -1]ᵀ ⊙ ∇(Δ⁻¹u)) ⋅ ∇u = λu + ν Δu + f
    ```

    For a detailed description of the terms, see the documentation of the
    `NavierStokesVorticity` stepper. The only difference is the additional
    forcing term `f` which injects new energy into the system. For a
    Kolmogorov flow **in primary variables**, it has the form

    ```
        f₀ = γ sin(k (2π/L) x₁)

        f₁ = 0
    ```

    In words, only the first channel is forced at a specific wavenumber over
    the second axis. Since this stepper considers the
    streamfunction-vorticity formulation, we take its curl to get

    ```
        f = -k (2π/L) γ cos(k (2π/L) x₁)
    ```

    The expected temporal behavior is that the initial vorticity field first
    is excited into a noisy striped pattern. This pattern breaks up and a
    turbulent spatio-temporal chaos emerges.

    A negative drag coefficient `λ` is needed to remove some of the energy
    piling up in low modes.

    According to

        Chandler, G.J. and Kerswell, R.R. (2013) ‘Invariant recurrent
        solutions embedded in a turbulent two-dimensional Kolmogorov flow’,
        Journal of Fluid Mechanics, 722, pp. 554–595.
        doi:10.1017/jfm.2013.122.

    equation (2.5), the Reynolds number of the Kolmogorov flow is given by

        Re = √ζ / ν √(L / (2π))³

    with `ζ` being the scaling of the Kolmogorov forcing, i.e., the
    `injection_scale`. Hence, in the case of `L = 2π`, `ζ = 1`, the Reynolds
    number is `Re = 1 / ν`. If one uses the default value of `ν = 0.001`,
    the Reynolds number is `Re = 1000` which also corresponds to the main
    experiments in

        Kochkov, D., Smith, J.A., Alieva, A., Wang, Q., Brenner, M.P. and
        Hoyer, S., 2021. Machine learning–accelerated computational fluid
        dynamics. Proceedings of the National Academy of Sciences, 118(21),
        p.e2101784118.

    together with `injection_mode = 4`. Note that they required a resolution
    of `num_points = 2048` (=> 2048^2 = 4.2M degrees of freedom in 2d) to
    fully resolve all scales at that Reynolds number. Using `Re = 0.01`
    which corresponds to `ν = 0.01` can be a good starting for
    `num_points=128`.


    **Arguments:**

    - `num_spatial_dims`: The number of spatial dimensions `d`.
    - `domain_extent`: The size of the domain `L`; in higher dimensions
        the domain is assumed to be a scaled hypercube `Ω = (0, L)ᵈ`.
    - `num_points`: The number of points `N` used to discretize the
        domain. This **includes** the left boundary point and **excludes**
        the right boundary point. In higher dimensions; the number of points
        in each dimension is the same. Hence, the total number of degrees of
        freedom is `Nᵈ`.
    - `dt`: The timestep size `Δt` between two consecutive states.
    - `diffusivity`: The diffusivity coefficient `ν`. This affects the
        Reynolds number. The lower the diffusivity, the "more turbulent".
        Default is `0.001`.
    - `convection_scale`: The scaling factor for the vorticity
        convection term. Default is `1.0`.
    - `drag`: The drag coefficient `λ`. Default is `-0.1`.
    - `injection_mode`: The mode of the injection. Default is `4`.
    - `injection_scale`: The scaling factor for the injection. Default
        is `1.0`.
    - `order`: The order of the Exponential Time Differencing Runge
        Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
        solves the linear part of the equation. Use higher values for higher
        accuracy and stability. The default choice of `2` is a good
        compromise for single precision floats.
    - `dealiasing_fraction`: The fraction of the wavenumbers to keep
        before evaluating the nonlinearity. The default 2/3 corresponds to
        Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
        2/3.
    - `num_circle_points`: How many points to use in the complex contour
        integral method to compute the coefficients of the exponential time
        differencing Runge Kutta method. Default: 16.
    - `circle_radius`: The radius of the contour used to compute the
        coefficients of the exponential time differencing Runge Kutta
        method. Default: 1.0.

    **Notes:**

    - Interesting pointer to explore:
      http://trieste-conf.itp.ac.ru/Boffetta.pdf
    """
    if num_spatial_dims != 2:
        raise ValueError(
            f"Expected num_spatial_dims = 2, got {num_spatial_dims}. "
            "For 3D, use KolmogorovFlowVelocity instead."
        )
    self.diffusivity = diffusivity
    self.convection_scale = convection_scale
    self.drag = drag
    self.injection_mode = injection_mode
    self.injection_scale = injection_scale
    self.dealiasing_fraction = dealiasing_fraction
    super().__init__(
        num_spatial_dims=num_spatial_dims,
        domain_extent=domain_extent,
        num_points=num_points,
        dt=dt,
        num_channels=1,
        order=order,
        num_circle_points=num_circle_points,
        circle_radius=circle_radius,
    )
__call__ ¤
__call__(
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]

Perform one step of the time integration for a single state.

Arguments:

  • u: The state vector, shape (C, ..., N,).

Returns:

  • u_next: The state vector after one step, shape (C, ..., N,).

Tip

Use this call method together with exponax.rollout to efficiently produce temporal trajectories.

Info

For batched operation, use jax.vmap on this function.

Source code in exponax/_base_stepper.py 
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
269
270
271
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Perform one step of the time integration for a single state.

    **Arguments:**

    - `u`: The state vector, shape `(C, ..., N,)`.

    **Returns:**

    - `u_next`: The state vector after one step, shape `(C, ..., N,)`.

    !!! tip
        Use this call method together with `exponax.rollout` to efficiently
        produce temporal trajectories.

    !!! info
        For batched operation, use `jax.vmap` on this function.
    """
    expected_shape = (self.num_channels,) + spatial_shape(
        self.num_spatial_dims, self.num_points
    )
    if u.shape != expected_shape:
        raise ValueError(
            f"""Expected shape {expected_shape}, got {u.shape}. For batched
             operation use `jax.vmap` on this function."""
        )
    return self.step(u)

3D: Velocity Formulation with Leray Projection¤

In 3D, the vorticity is also a three-component vector, offering no reduction in degrees of freedom. The velocity formulation with Leray projection is preferred instead, using the rotational form of the convection term.

exponax.stepper.NavierStokesVelocity ¤

Bases: BaseStepper

Source code in exponax/stepper/_navier_stokes.py 
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
class NavierStokesVelocity(BaseStepper):
    diffusivity: float
    drag: float
    dealiasing_fraction: float

    def __init__(
        self,
        num_spatial_dims: int,
        domain_extent: float,
        num_points: int,
        dt: float,
        *,
        diffusivity: float = 0.01,
        drag: float = 0.0,
        order: int = 2,
        dealiasing_fraction: float = 2 / 3,
        num_circle_points: int = 16,
        circle_radius: float = 1.0,
    ):
        """
        Timestepper for the 3d incompressible Navier-Stokes equations on
        periodic boundary conditions in velocity formulation. The equation
        reads

        ```
            uₜ = ν Δu + λu + 𝒫(u × ω)
        ```

        with `u` the three-channel velocity field and `ω = ∇ × u` the
        vorticity. The term `𝒫` denotes the Leray projection which enforces
        the incompressibility constraint `∇ ⋅ u = 0` by removing gradient
        components (pressure and kinetic energy gradient). The first term on
        the right-hand side is a diffusion with coefficient `ν` and the second
        term is an optional drag with coefficient `λ`.

        The nonlinear term uses the rotational form which exploits the vector
        identity `(u ⋅ ∇)u = ∇(|u|²/2) + ω × u` and the fact that the Leray
        projection annihilates all gradient fields.

        The Reynolds number of the problem is `Re = U L / ν` with `U` a
        characteristic velocity scale and `L` the `domain_extent`.

        **Arguments:**

        - `num_spatial_dims`: The number of spatial dimensions `d`. Must be
            `3`.
        - `domain_extent`: The size of the domain `L`; the domain is assumed
            to be a scaled hypercube `Ω = (0, L)³`.
        - `num_points`: The number of points `N` used to discretize the
            domain. This **includes** the left boundary point and **excludes**
            the right boundary point. The number of points in each dimension
            is the same. Hence, the total number of degrees of freedom is
            `N³`.
        - `dt`: The timestep size `Δt` between two consecutive states.
        - `diffusivity`: The diffusivity (viscosity) coefficient `ν`. This
            affects the Reynolds number. The lower the diffusivity, the "more
            turbulent". Default is `0.01`.
        - `drag`: The drag coefficient `λ`. Default is `0.0`.
        - `order`: The order of the Exponential Time Differencing Runge
            Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
            solves the linear part of the equation. Use higher values for
            higher accuracy and stability. The default choice of `2` is a good
            compromise for single precision floats.
        - `dealiasing_fraction`: The fraction of the wavenumbers to keep
            before evaluating the nonlinearity. The default 2/3 corresponds to
            Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
            2/3.
        - `num_circle_points`: How many points to use in the complex contour
            integral method to compute the coefficients of the exponential time
            differencing Runge Kutta method. Default: 16.
        - `circle_radius`: The radius of the contour used to compute the
            coefficients of the exponential time differencing Runge Kutta
            method. Default: 1.0.

        **Notes:**

        - In 3d, the velocity formulation is preferred over the vorticity
            formulation because the vorticity is also a three-component vector,
            offering no reduction in degrees of freedom. For 2d, use
            `NavierStokesVorticity` which solves for the scalar vorticity
            instead.
        - The nonlinear term uses the rotational form `𝒫(u × ω)` rather
            than the convective form `𝒫(-(u ⋅ ∇)u)`. Both are equivalent at
            the continuous level, but the rotational form is preferred for
            pseudo-spectral methods because it requires fewer FFTs. In the
            convective form, computing `(u ⋅ ∇)u` requires for each of the 3
            velocity components a dot product of `u` with its gradient, i.e.,
            `Σⱼ uⱼ ∂uᵢ/∂xⱼ`. That amounts to 3 × 3 = 9 physical-space
            multiplications, each requiring an inverse FFT for the operand and
            a forward FFT for the result. In the rotational form, the curl
            `ω = ∇ × u` is free in Fourier space (a cross product with `ik`),
            and the subsequent cross product `u × ω` in physical space
            involves only 6 multiplications (two per output component).
            Additionally, the gradient terms `∇(|u|²/2 + p)` are implicitly
            eliminated by the Leray projection without ever being computed.
            Since the FFTs are the most computationally demanding operations
            in higher dimensions (scaling as `O(N³ log N)` in 3d), reducing
            their count directly improves performance.
        - The higher the Reynolds number, the smaller the timestep size must
            be to ensure stability.
        """
        if num_spatial_dims != 3:
            raise ValueError(
                f"Expected num_spatial_dims = 3, got {num_spatial_dims}. "
                "For 2D, use NavierStokesVorticity instead."
            )

        self.diffusivity = diffusivity
        self.drag = drag
        self.dealiasing_fraction = dealiasing_fraction

        super().__init__(
            num_spatial_dims=num_spatial_dims,
            domain_extent=domain_extent,
            num_points=num_points,
            dt=dt,
            num_channels=3,
            order=order,
            num_circle_points=num_circle_points,
            circle_radius=circle_radius,
        )

    def _build_linear_operator(self, derivative_operator):
        return self.diffusivity * build_laplace_operator(
            derivative_operator, order=2
        ) + self.drag * build_laplace_operator(derivative_operator, order=0)

    def _build_nonlinear_fun(self, derivative_operator):
        return ProjectedConvection3d(
            num_spatial_dims=self.num_spatial_dims,
            num_points=self.num_points,
            derivative_operator=derivative_operator,
            dealiasing_fraction=self.dealiasing_fraction,
        )
__init__ ¤
__init__(
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    drag: float = 0.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0
)

Timestepper for the 3d incompressible Navier-Stokes equations on periodic boundary conditions in velocity formulation. The equation reads

    uₜ = ν Δu + λu + 𝒫(u × ω)

with u the three-channel velocity field and ω = ∇ × u the vorticity. The term 𝒫 denotes the Leray projection which enforces the incompressibility constraint ∇ ⋅ u = 0 by removing gradient components (pressure and kinetic energy gradient). The first term on the right-hand side is a diffusion with coefficient ν and the second term is an optional drag with coefficient λ.

The nonlinear term uses the rotational form which exploits the vector identity (u ⋅ ∇)u = ∇(|u|²/2) + ω × u and the fact that the Leray projection annihilates all gradient fields.

The Reynolds number of the problem is Re = U L / ν with U a characteristic velocity scale and L the domain_extent.

Arguments:

  • num_spatial_dims: The number of spatial dimensions d. Must be 3.
  • domain_extent: The size of the domain L; the domain is assumed to be a scaled hypercube Ω = (0, L)³.
  • num_points: The number of points N used to discretize the domain. This includes the left boundary point and excludes the right boundary point. The number of points in each dimension is the same. Hence, the total number of degrees of freedom is N³.
  • dt: The timestep size Δt between two consecutive states.
  • diffusivity: The diffusivity (viscosity) coefficient ν. This affects the Reynolds number. The lower the diffusivity, the "more turbulent". Default is 0.01.
  • drag: The drag coefficient λ. Default is 0.0.
  • order: The order of the Exponential Time Differencing Runge Kutta method. Must be one of {0, 1, 2, 3, 4}. The option 0 only solves the linear part of the equation. Use higher values for higher accuracy and stability. The default choice of 2 is a good compromise for single precision floats.
  • dealiasing_fraction: The fraction of the wavenumbers to keep before evaluating the nonlinearity. The default 2/3 corresponds to Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default: 2/3.
  • num_circle_points: How many points to use in the complex contour integral method to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 16.
  • circle_radius: The radius of the contour used to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 1.0.

Notes:

  • In 3d, the velocity formulation is preferred over the vorticity formulation because the vorticity is also a three-component vector, offering no reduction in degrees of freedom. For 2d, use NavierStokesVorticity which solves for the scalar vorticity instead.
  • The nonlinear term uses the rotational form 𝒫(u × ω) rather than the convective form 𝒫(-(u ⋅ ∇)u). Both are equivalent at the continuous level, but the rotational form is preferred for pseudo-spectral methods because it requires fewer FFTs. In the convective form, computing (u ⋅ ∇)u requires for each of the 3 velocity components a dot product of u with its gradient, i.e., Σⱼ uⱼ ∂uᵢ/∂xⱼ. That amounts to 3 × 3 = 9 physical-space multiplications, each requiring an inverse FFT for the operand and a forward FFT for the result. In the rotational form, the curl ω = ∇ × u is free in Fourier space (a cross product with ik), and the subsequent cross product u × ω in physical space involves only 6 multiplications (two per output component). Additionally, the gradient terms ∇(|u|²/2 + p) are implicitly eliminated by the Leray projection without ever being computed. Since the FFTs are the most computationally demanding operations in higher dimensions (scaling as O(N³ log N) in 3d), reducing their count directly improves performance.
  • The higher the Reynolds number, the smaller the timestep size must be to ensure stability.
Source code in exponax/stepper/_navier_stokes.py 
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
def __init__(
    self,
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    drag: float = 0.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0,
):
    """
    Timestepper for the 3d incompressible Navier-Stokes equations on
    periodic boundary conditions in velocity formulation. The equation
    reads

    ```
        uₜ = ν Δu + λu + 𝒫(u × ω)
    ```

    with `u` the three-channel velocity field and `ω = ∇ × u` the
    vorticity. The term `𝒫` denotes the Leray projection which enforces
    the incompressibility constraint `∇ ⋅ u = 0` by removing gradient
    components (pressure and kinetic energy gradient). The first term on
    the right-hand side is a diffusion with coefficient `ν` and the second
    term is an optional drag with coefficient `λ`.

    The nonlinear term uses the rotational form which exploits the vector
    identity `(u ⋅ ∇)u = ∇(|u|²/2) + ω × u` and the fact that the Leray
    projection annihilates all gradient fields.

    The Reynolds number of the problem is `Re = U L / ν` with `U` a
    characteristic velocity scale and `L` the `domain_extent`.

    **Arguments:**

    - `num_spatial_dims`: The number of spatial dimensions `d`. Must be
        `3`.
    - `domain_extent`: The size of the domain `L`; the domain is assumed
        to be a scaled hypercube `Ω = (0, L)³`.
    - `num_points`: The number of points `N` used to discretize the
        domain. This **includes** the left boundary point and **excludes**
        the right boundary point. The number of points in each dimension
        is the same. Hence, the total number of degrees of freedom is
        `N³`.
    - `dt`: The timestep size `Δt` between two consecutive states.
    - `diffusivity`: The diffusivity (viscosity) coefficient `ν`. This
        affects the Reynolds number. The lower the diffusivity, the "more
        turbulent". Default is `0.01`.
    - `drag`: The drag coefficient `λ`. Default is `0.0`.
    - `order`: The order of the Exponential Time Differencing Runge
        Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
        solves the linear part of the equation. Use higher values for
        higher accuracy and stability. The default choice of `2` is a good
        compromise for single precision floats.
    - `dealiasing_fraction`: The fraction of the wavenumbers to keep
        before evaluating the nonlinearity. The default 2/3 corresponds to
        Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
        2/3.
    - `num_circle_points`: How many points to use in the complex contour
        integral method to compute the coefficients of the exponential time
        differencing Runge Kutta method. Default: 16.
    - `circle_radius`: The radius of the contour used to compute the
        coefficients of the exponential time differencing Runge Kutta
        method. Default: 1.0.

    **Notes:**

    - In 3d, the velocity formulation is preferred over the vorticity
        formulation because the vorticity is also a three-component vector,
        offering no reduction in degrees of freedom. For 2d, use
        `NavierStokesVorticity` which solves for the scalar vorticity
        instead.
    - The nonlinear term uses the rotational form `𝒫(u × ω)` rather
        than the convective form `𝒫(-(u ⋅ ∇)u)`. Both are equivalent at
        the continuous level, but the rotational form is preferred for
        pseudo-spectral methods because it requires fewer FFTs. In the
        convective form, computing `(u ⋅ ∇)u` requires for each of the 3
        velocity components a dot product of `u` with its gradient, i.e.,
        `Σⱼ uⱼ ∂uᵢ/∂xⱼ`. That amounts to 3 × 3 = 9 physical-space
        multiplications, each requiring an inverse FFT for the operand and
        a forward FFT for the result. In the rotational form, the curl
        `ω = ∇ × u` is free in Fourier space (a cross product with `ik`),
        and the subsequent cross product `u × ω` in physical space
        involves only 6 multiplications (two per output component).
        Additionally, the gradient terms `∇(|u|²/2 + p)` are implicitly
        eliminated by the Leray projection without ever being computed.
        Since the FFTs are the most computationally demanding operations
        in higher dimensions (scaling as `O(N³ log N)` in 3d), reducing
        their count directly improves performance.
    - The higher the Reynolds number, the smaller the timestep size must
        be to ensure stability.
    """
    if num_spatial_dims != 3:
        raise ValueError(
            f"Expected num_spatial_dims = 3, got {num_spatial_dims}. "
            "For 2D, use NavierStokesVorticity instead."
        )

    self.diffusivity = diffusivity
    self.drag = drag
    self.dealiasing_fraction = dealiasing_fraction

    super().__init__(
        num_spatial_dims=num_spatial_dims,
        domain_extent=domain_extent,
        num_points=num_points,
        dt=dt,
        num_channels=3,
        order=order,
        num_circle_points=num_circle_points,
        circle_radius=circle_radius,
    )
__call__ ¤
__call__(
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]

Perform one step of the time integration for a single state.

Arguments:

  • u: The state vector, shape (C, ..., N,).

Returns:

  • u_next: The state vector after one step, shape (C, ..., N,).

Tip

Use this call method together with exponax.rollout to efficiently produce temporal trajectories.

Info

For batched operation, use jax.vmap on this function.

Source code in exponax/_base_stepper.py 
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
269
270
271
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Perform one step of the time integration for a single state.

    **Arguments:**

    - `u`: The state vector, shape `(C, ..., N,)`.

    **Returns:**

    - `u_next`: The state vector after one step, shape `(C, ..., N,)`.

    !!! tip
        Use this call method together with `exponax.rollout` to efficiently
        produce temporal trajectories.

    !!! info
        For batched operation, use `jax.vmap` on this function.
    """
    expected_shape = (self.num_channels,) + spatial_shape(
        self.num_spatial_dims, self.num_points
    )
    if u.shape != expected_shape:
        raise ValueError(
            f"""Expected shape {expected_shape}, got {u.shape}. For batched
             operation use `jax.vmap` on this function."""
        )
    return self.step(u)

exponax.stepper.KolmogorovFlowVelocity ¤

Bases: BaseStepper

Source code in exponax/stepper/_navier_stokes.py 
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
class KolmogorovFlowVelocity(BaseStepper):
    diffusivity: float
    drag: float
    injection_mode: int
    injection_scale: float
    dealiasing_fraction: float

    def __init__(
        self,
        num_spatial_dims: int,
        domain_extent: float,
        num_points: int,
        dt: float,
        *,
        diffusivity: float = 0.01,
        drag: float = 0.0,
        injection_mode: int = 4,
        injection_scale: float = 1.0,
        order: int = 2,
        dealiasing_fraction: float = 2 / 3,
        num_circle_points: int = 16,
        circle_radius: float = 1.0,
    ):
        """
        Timestepper for the 3d Kolmogorov flow equation on periodic boundary
        conditions in velocity formulation. The equation reads

        ```
            uₜ = ν Δu + λu + 𝒫(u × ω) + f
        ```

        For a detailed description of the terms, see the documentation of the
        `NavierStokesVelocity` stepper. The only difference is the additional
        forcing term `f` which injects new energy into the system. It has the
        form

        ```
            f₀ = γ sin(k (2π/L) x₁)

            f₁ = 0

            f₂ = 0
        ```

        In words, only the first velocity channel is forced at a specific
        wavenumber over the second spatial axis. This forcing is divergence-
        free because the forced component does not vary along its own
        direction.

        The expected temporal behavior is that the velocity field develops
        shear layers which become unstable and break up into turbulent
        spatio-temporal chaos.

        A negative drag coefficient `λ` is needed to remove some of the energy
        piling up in low modes.

        **Arguments:**

        - `num_spatial_dims`: The number of spatial dimensions `d`. Must be
            `3`.
        - `domain_extent`: The size of the domain `L`; the domain is assumed
            to be a scaled hypercube `Ω = (0, L)³`.
        - `num_points`: The number of points `N` used to discretize the
            domain. This **includes** the left boundary point and **excludes**
            the right boundary point. The number of points in each dimension
            is the same. Hence, the total number of degrees of freedom is
            `N³`.
        - `dt`: The timestep size `Δt` between two consecutive states.
        - `diffusivity`: The diffusivity (viscosity) coefficient `ν`. This
            affects the Reynolds number. The lower the diffusivity, the "more
            turbulent". Default is `0.01`.
        - `drag`: The drag coefficient `λ`. Default is `0.0`.
        - `injection_mode`: The wavenumber `k` at which energy is injected.
            Default is `4`.
        - `injection_scale`: The intensity `γ` of the injection term. Default
            is `1.0`.
        - `order`: The order of the Exponential Time Differencing Runge
            Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
            solves the linear part of the equation. Use higher values for
            higher accuracy and stability. The default choice of `2` is a good
            compromise for single precision floats.
        - `dealiasing_fraction`: The fraction of the wavenumbers to keep
            before evaluating the nonlinearity. The default 2/3 corresponds to
            Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
            2/3.
        - `num_circle_points`: How many points to use in the complex contour
            integral method to compute the coefficients of the exponential time
            differencing Runge Kutta method. Default: 16.
        - `circle_radius`: The radius of the contour used to compute the
            coefficients of the exponential time differencing Runge Kutta
            method. Default: 1.0.
        """
        if num_spatial_dims != 3:
            raise ValueError(
                f"Expected num_spatial_dims = 3, got {num_spatial_dims}. "
                "For 2D, use KolmogorovFlowVorticity instead."
            )
        self.diffusivity = diffusivity
        self.drag = drag
        self.injection_mode = injection_mode
        self.injection_scale = injection_scale
        self.dealiasing_fraction = dealiasing_fraction
        super().__init__(
            num_spatial_dims=num_spatial_dims,
            domain_extent=domain_extent,
            num_points=num_points,
            dt=dt,
            num_channels=3,
            order=order,
            num_circle_points=num_circle_points,
            circle_radius=circle_radius,
        )

    def _build_linear_operator(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> Complex[Array, "1 ... (N//2)+1"]:
        return self.diffusivity * build_laplace_operator(
            derivative_operator, order=2
        ) + self.drag * build_laplace_operator(derivative_operator, order=0)

    def _build_nonlinear_fun(
        self,
        derivative_operator: Complex[Array, "D ... (N//2)+1"],
    ) -> ProjectedConvection3dKolmogorov:
        return ProjectedConvection3dKolmogorov(
            self.num_spatial_dims,
            self.num_points,
            injection_mode=self.injection_mode,
            injection_scale=self.injection_scale,
            derivative_operator=derivative_operator,
            dealiasing_fraction=self.dealiasing_fraction,
        )
__init__ ¤
__init__(
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    drag: float = 0.0,
    injection_mode: int = 4,
    injection_scale: float = 1.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0
)

Timestepper for the 3d Kolmogorov flow equation on periodic boundary conditions in velocity formulation. The equation reads

    uₜ = ν Δu + λu + 𝒫(u × ω) + f

For a detailed description of the terms, see the documentation of the NavierStokesVelocity stepper. The only difference is the additional forcing term f which injects new energy into the system. It has the form

    f₀ = γ sin(k (2π/L) x₁)

    f₁ = 0

    f₂ = 0

In words, only the first velocity channel is forced at a specific wavenumber over the second spatial axis. This forcing is divergence- free because the forced component does not vary along its own direction.

The expected temporal behavior is that the velocity field develops shear layers which become unstable and break up into turbulent spatio-temporal chaos.

A negative drag coefficient λ is needed to remove some of the energy piling up in low modes.

Arguments:

  • num_spatial_dims: The number of spatial dimensions d. Must be 3.
  • domain_extent: The size of the domain L; the domain is assumed to be a scaled hypercube Ω = (0, L)³.
  • num_points: The number of points N used to discretize the domain. This includes the left boundary point and excludes the right boundary point. The number of points in each dimension is the same. Hence, the total number of degrees of freedom is N³.
  • dt: The timestep size Δt between two consecutive states.
  • diffusivity: The diffusivity (viscosity) coefficient ν. This affects the Reynolds number. The lower the diffusivity, the "more turbulent". Default is 0.01.
  • drag: The drag coefficient λ. Default is 0.0.
  • injection_mode: The wavenumber k at which energy is injected. Default is 4.
  • injection_scale: The intensity γ of the injection term. Default is 1.0.
  • order: The order of the Exponential Time Differencing Runge Kutta method. Must be one of {0, 1, 2, 3, 4}. The option 0 only solves the linear part of the equation. Use higher values for higher accuracy and stability. The default choice of 2 is a good compromise for single precision floats.
  • dealiasing_fraction: The fraction of the wavenumbers to keep before evaluating the nonlinearity. The default 2/3 corresponds to Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default: 2/3.
  • num_circle_points: How many points to use in the complex contour integral method to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 16.
  • circle_radius: The radius of the contour used to compute the coefficients of the exponential time differencing Runge Kutta method. Default: 1.0.
Source code in exponax/stepper/_navier_stokes.py 
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
def __init__(
    self,
    num_spatial_dims: int,
    domain_extent: float,
    num_points: int,
    dt: float,
    *,
    diffusivity: float = 0.01,
    drag: float = 0.0,
    injection_mode: int = 4,
    injection_scale: float = 1.0,
    order: int = 2,
    dealiasing_fraction: float = 2 / 3,
    num_circle_points: int = 16,
    circle_radius: float = 1.0,
):
    """
    Timestepper for the 3d Kolmogorov flow equation on periodic boundary
    conditions in velocity formulation. The equation reads

    ```
        uₜ = ν Δu + λu + 𝒫(u × ω) + f
    ```

    For a detailed description of the terms, see the documentation of the
    `NavierStokesVelocity` stepper. The only difference is the additional
    forcing term `f` which injects new energy into the system. It has the
    form

    ```
        f₀ = γ sin(k (2π/L) x₁)

        f₁ = 0

        f₂ = 0
    ```

    In words, only the first velocity channel is forced at a specific
    wavenumber over the second spatial axis. This forcing is divergence-
    free because the forced component does not vary along its own
    direction.

    The expected temporal behavior is that the velocity field develops
    shear layers which become unstable and break up into turbulent
    spatio-temporal chaos.

    A negative drag coefficient `λ` is needed to remove some of the energy
    piling up in low modes.

    **Arguments:**

    - `num_spatial_dims`: The number of spatial dimensions `d`. Must be
        `3`.
    - `domain_extent`: The size of the domain `L`; the domain is assumed
        to be a scaled hypercube `Ω = (0, L)³`.
    - `num_points`: The number of points `N` used to discretize the
        domain. This **includes** the left boundary point and **excludes**
        the right boundary point. The number of points in each dimension
        is the same. Hence, the total number of degrees of freedom is
        `N³`.
    - `dt`: The timestep size `Δt` between two consecutive states.
    - `diffusivity`: The diffusivity (viscosity) coefficient `ν`. This
        affects the Reynolds number. The lower the diffusivity, the "more
        turbulent". Default is `0.01`.
    - `drag`: The drag coefficient `λ`. Default is `0.0`.
    - `injection_mode`: The wavenumber `k` at which energy is injected.
        Default is `4`.
    - `injection_scale`: The intensity `γ` of the injection term. Default
        is `1.0`.
    - `order`: The order of the Exponential Time Differencing Runge
        Kutta method. Must be one of {0, 1, 2, 3, 4}. The option `0` only
        solves the linear part of the equation. Use higher values for
        higher accuracy and stability. The default choice of `2` is a good
        compromise for single precision floats.
    - `dealiasing_fraction`: The fraction of the wavenumbers to keep
        before evaluating the nonlinearity. The default 2/3 corresponds to
        Orszag's 2/3 rule. To fully eliminate aliasing, use 1/2. Default:
        2/3.
    - `num_circle_points`: How many points to use in the complex contour
        integral method to compute the coefficients of the exponential time
        differencing Runge Kutta method. Default: 16.
    - `circle_radius`: The radius of the contour used to compute the
        coefficients of the exponential time differencing Runge Kutta
        method. Default: 1.0.
    """
    if num_spatial_dims != 3:
        raise ValueError(
            f"Expected num_spatial_dims = 3, got {num_spatial_dims}. "
            "For 2D, use KolmogorovFlowVorticity instead."
        )
    self.diffusivity = diffusivity
    self.drag = drag
    self.injection_mode = injection_mode
    self.injection_scale = injection_scale
    self.dealiasing_fraction = dealiasing_fraction
    super().__init__(
        num_spatial_dims=num_spatial_dims,
        domain_extent=domain_extent,
        num_points=num_points,
        dt=dt,
        num_channels=3,
        order=order,
        num_circle_points=num_circle_points,
        circle_radius=circle_radius,
    )
__call__ ¤
__call__(
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]

Perform one step of the time integration for a single state.

Arguments:

  • u: The state vector, shape (C, ..., N,).

Returns:

  • u_next: The state vector after one step, shape (C, ..., N,).

Tip

Use this call method together with exponax.rollout to efficiently produce temporal trajectories.

Info

For batched operation, use jax.vmap on this function.

Source code in exponax/_base_stepper.py 
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
269
270
271
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Perform one step of the time integration for a single state.

    **Arguments:**

    - `u`: The state vector, shape `(C, ..., N,)`.

    **Returns:**

    - `u_next`: The state vector after one step, shape `(C, ..., N,)`.

    !!! tip
        Use this call method together with `exponax.rollout` to efficiently
        produce temporal trajectories.

    !!! info
        For batched operation, use `jax.vmap` on this function.
    """
    expected_shape = (self.num_channels,) + spatial_shape(
        self.num_spatial_dims, self.num_points
    )
    if u.shape != expected_shape:
        raise ValueError(
            f"""Expected shape {expected_shape}, got {u.shape}. For batched
             operation use `jax.vmap` on this function."""
        )
    return self.step(u)