FourierFlows / FourierFlows.jl

Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains
https://bit.ly/FourierFlows
MIT License
207 stars 29 forks source link

Bug in `FourierFlows.getaliasedwavenumbers`? #284

Closed navidcy closed 3 years ago

navidcy commented 3 years ago

I feel that https://github.com/FourierFlows/FourierFlows.jl/blob/a4bcb4644169a125156f2a8426348a5cb0cb7da8/src/domains.jl#L312 contains a typo and it should instead be Int(nk/2+1).

glwagner commented 3 years ago

Doesn't nx=nk?

navidcy commented 3 years ago

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 commented 3 years ago

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

glwagner commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

I don't see nx in scope:

https://github.com/FourierFlows/FourierFlows.jl/blob/a4bcb4644169a125156f2a8426348a5cb0cb7da8/src/domains.jl#L302-L316

What am I missing?

glwagner commented 3 years ago

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 commented 3 years ago

I don't see nx in scope:

https://github.com/FourierFlows/FourierFlows.jl/blob/a4bcb4644169a125156f2a8426348a5cb0cb7da8/src/domains.jl#L302-L316

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 commented 3 years ago

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 commented 3 years ago

What about the case aliasfraction=0?

navidcy commented 3 years ago

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 commented 3 years ago

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.