Skip to content

Navier-Stokes (in streamfunction-vorticity formulation)¤

exponax.stepper.NavierStokesVorticity ¤

Bases: BaseStepper

Source code in exponax/stepper/_navier_stokes.py 
  8
  9
 10
 11
 12
 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
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 scale 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 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}.")

        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 scale 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 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 
 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
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 scale 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 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}.")

    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"]

Performs a check

Source code in exponax/_base_stepper.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Performs a check
    """
    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 
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
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
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.
        """
        if num_spatial_dims != 2:
            raise ValueError(f"Expected num_spatial_dims = 2, got {num_spatial_dims}.")
        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.
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
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.
    """
    if num_spatial_dims != 2:
        raise ValueError(f"Expected num_spatial_dims = 2, got {num_spatial_dims}.")
    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"]

Performs a check

Source code in exponax/_base_stepper.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def __call__(
    self,
    u: Float[Array, "C ... N"],
) -> Float[Array, "C ... N"]:
    """
    Performs a check
    """
    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)