Bug in `FourierFlows.getaliasedwavenumbers`?

[navidcy]

I feel that

FourierFlows.jl/src/domains.jl

Line 312 in a4bcb46

kalias = (aliasfraction > 0) ? (iL:iR) : Int(nx/2 + 1)

contains a typo and it should instead be Int(nk/2+1).

[glwagner]

Doesn't nx=nk?

[navidcy]

Yes, but nx is not an argument of the getaliasedwavenumbers(). And also, we use this getaliasedwavenumbers() function to find the aliased wavenumbers in y and z...

[navidcy]

Perhaps when this function was first made we were only thinking of x?

[glwagner]

Yes, but nx is not an argument of the getaliasedwavenumbers(). And also, we use this getaliasedwavenumbers() function to find the aliased wavenumbers in y and z...

It's independent of the direction. We have nx=nk, ny=nl, and nz=nm. This function is 1D.

But it certainly would throw an error that nx is not an argument to the function.

I was just confused about what the issue was. The issue is not that nk is used rather than nx; the issue is that nx is not available in the function.

[navidcy]

Yes, but the function was not throwing any error. However, when I called it as

malias, _ = getaliasedwavenumbers(nm, nm, aliasfraction)

to get the aliased wavenumbers in z it was using nx I think and was giving me out nonsense if nz ≠ nx...

Sorry, I should have been more explicit in the bug description.

[glwagner]

I don't see nx in scope:

FourierFlows.jl/src/domains.jl

Lines 302 to 316 in a4bcb46

	function getaliasedwavenumbers(nk, nkr, aliasfraction)
	# Index endpoints for aliased i, j wavenumbers
	# 1/3 aliasfraction => upper 1/6 of +/- wavenumbers (1/3 total) are set to 0 after performing fft.
	# 1/2 aliasfraction => upper 1/4 of +/- wavenumbers (1/2 total) are set to 0 after performing fft.
	L = (1 - aliasfraction)/2 # (1 - 1/3) / 2 + 1 = 1/3.
	R = (1 + aliasfraction)/2 # (1 + 1/3) / 2 - 1 = 2/3.
	iL = floor(Int, L * nk) + 1
	iR = ceil(Int, R * nk)
	aliasfraction < 1 || error("`aliasfraction` must be less than 1") # aliasfraction=1 is not sensible.
	kalias = (aliasfraction > 0) ? (iL:iR) : Int(nx/2 + 1)
	kralias = (aliasfraction > 0) ? (iL:nkr) : nkr
	return kalias, kralias
	end

What am I missing?

[glwagner]

Perhaps we do not have a test that accesses the second branch in line 312? It would seem that we usually call line 313 in which case there is no error.

Fine to merge either way of course.

[navidcy]

I don't see nx in scope:

FourierFlows.jl/src/domains.jl

Lines 302 to 316 in a4bcb46

	function getaliasedwavenumbers(nk, nkr, aliasfraction)
	# Index endpoints for aliased i, j wavenumbers
	# 1/3 aliasfraction => upper 1/6 of +/- wavenumbers (1/3 total) are set to 0 after performing fft.
	# 1/2 aliasfraction => upper 1/4 of +/- wavenumbers (1/2 total) are set to 0 after performing fft.
	L = (1 - aliasfraction)/2 # (1 - 1/3) / 2 + 1 = 1/3.
	R = (1 + aliasfraction)/2 # (1 + 1/3) / 2 - 1 = 2/3.
	iL = floor(Int, L * nk) + 1
	iR = ceil(Int, R * nk)
	aliasfraction < 1 || error("`aliasfraction` must be less than 1") # aliasfraction=1 is not sensible.
	kalias = (aliasfraction > 0) ? (iL:iR) : Int(nx/2 + 1)
	kralias = (aliasfraction > 0) ? (iL:nkr) : nkr
	return kalias, kralias
	end

What am I missing?

You are asking why the function does not complain? I guess because we always call it from within the grid constructor and nx is in the global scope of the constructor? It takes the value from there?

[navidcy]

Actually now that I'm thinking about it, what is this second branch? I'd expect that a better structure of this function should be:

  aliasfraction < 1 || error("`aliasfraction` must be less than 1") # aliasfraction=1 is not sensible.
  aliasfraction >= 0 || error("`aliasfraction` must be non-negative")

kalias = iL:iR
kralias = iL:nkr

[glwagner]

What about the case aliasfraction=0?

[navidcy]

Shouldn't the second branch then be:

   kalias = (aliasfraction > 0) ? (iL:iR) : (1:Int(nx/2 + 1))
  kralias = (aliasfraction > 0) ? (iL:nkr) : (1:nkr)

because as it is now, it returns just the number nkr or nk.

[glwagner]

Hmm I think that's reversed logic: if aliasfraction = 0 then no wavenumbers are aliased. I guess nx/2 + 1 is a hack that just zeros out the highest wavenumber. In reality maybe it should return kalias=nothing or something.