Forcing implementation in physical space

[dlemasqu]

Dear all,

My name is Daphné, I'm a Lecturer in St Andrews and I'd like to use GeophysicalFlows.jl for idealized GFD problems (thank you for developing this great tool!).
I'm interested in simulating Single layer QG dynamics on a beta plane, with a narrow-band, isotropic, stochastic forcing centred at a given wavenumber, but localized in latitude, e.g., multiplied by a gaussian in y in physical space.

I managed to install and run the examples you provide on Github, but I did not manage to "localize" the forcing in the y direction.
I tried to modify the example script for "Forced-dissipative barotropic QG beta-plane turbulence". I modified the function calcF! and tried to transform to physical space, multiply by a gaussian, and transform back to spectral space (attached)

, but this fails due to DimensionMismatch. Is this because of dealiasing?

I'm probably doing something very stupid as I'm new to both Julia and GeophysicalFlows, so I was wondering if you would be able to provide some help to implement this forcing?

Thank you very much for your help!

Daphné

[navidcy]

Hi @dlemasqu!

The problem is not with the dealiasing but I do thing that it stems from that line.
But the issue is that you have the @. at the front which implies that all methods and operations are to be applied element-wise. And thus you are trying to apply the rfft element-wise? I think that's the issue...

This MWE:

using GeophysicalFlows, Random, GLMakie

n = 128
L = 2π

forcing_wavenumber = 14.0 * 2π/L  # the forcing wavenumber, `k_f`, for a spectrum that is a ring in wavenumber space
forcing_bandwidth  = 1.5  * 2π/L  # the width of the forcing spectrum, `δ_f`

grid = TwoDGrid(nx=n, Lx=L)

K = @. sqrt(grid.Krsq)            # a 2D array with the total wavenumber

forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
forcing_spectrum[grid.Krsq .== 0] .= 0 # ensure forcing has zero domain-average

X, Y = gridpoints(grid)

σy = 0.5
y0 = 0.1

gaussian_mask = @. exp( - (Y - y0)^2 / 2σy^2)

function calcF!(Fh, grid)
  randn!(Fh)
  @. Fh *= sqrt(forcing_spectrum)
  F = irfft(Fh, grid.nx)

  @. F *= gaussian_mask

  Fh .= rfft(F)

  return nothing
end

Fh = zeros(Complex{eltype(grid)}, (grid.nkr, grid.nl))

calcF!(Fh, grid)

heatmap(irfft(Fh, grid.nx))

produces

If you wanted to apply the timestepper filter you'd have to do something like:

  Fh .= prob.timestepper.filter .* rfft(F)

without the @. in front!

[navidcy]

On another note, calling irfft and rfft, e.g. on line

  F = irfft(Fh, grid.nx)

within the calcF! function implies that you are creating a new temporary arrays every time you want a new forcing realisation (which will be every tilmestep). This is a lot of memory allocations...

julia> @time calcF!(Fh, grid)
  0.001167 seconds (879 allocations: 580.609 KiB)

A way to make this code more efficient would be to create the array F once and then just modify it in place. Something like:

using LinearAlgebra: mul!, ldiv!

F = zeros(grid.nx, grid.ny)

function calcF!(Fh, grid)
  randn!(Fh)
  @. Fh *= sqrt(forcing_spectrum)
  ldiv!(F, grid.rfftplan, Fh)

  @. F *= gaussian_mask

  mul!(Fh, grid.rfftplan, F)

  return nothing
end

(See https://fourierflows.github.io/FourierFlowsDocumentation/dev/grids/ for how Fourier transforms with mul! and ldiv! come about.)

With this modification, the calcF! method is 2x faster with 10x less memory allocations

julia> @time calcF!(Fh, grid)
  0.000501 seconds (867 allocations: 61.641 KiB)

[dlemasqu]

Thanks so much Navid for looking into this and for the pedagogy, it makes sense to me now and works like a charm!