Dealias sol before computing nonlinear terms

[navidcy]

This PR ensures that before any calculation of quadratic nonlinear terms the state vector sol is properly dealiased.

Closes #247.

[glwagner]

Are there any docstrings that are updated? This changes default behavior and makes it so that filters can't be used without setting dealias=0 when constructing the grid. I think some of the examples use filters so they might need to be updated as well?

[navidcy]

Are there any docstrings that are updated? This changes default behavior and makes it so that filters can't be used without setting dealias=0 when constructing the grid. I think some of the examples use filters so they might need to be updated as well?

• This changes the equations that the module solve for I would say.. not "changes the default behavior". It just computes the nonlinear terms without any aliasing from high, spurious wavenumbers. Of course one could argue that a user may wanna simulate a flow with aliasing. To do that they may will have to either rewrite calcN! or set dealias=0 when constructing the grid. I'm simply arguing here that we were computing nonlinear terms wrong before and now they are correct.
• filters can be used even without setting dealias=0 when constructing the grid
• I agree that examples should be enhanced to better demonstrate/explain aliasing. I was thinking that this should be more appropriate in the FourierFlows.jl/Docs; see Docs section on dealiasing/filtered time-steppers FourierFlows.jl#283. What do you think?

[glwagner]

I think the exponential filters we have defined assume that the flow is not dealiased or won't work as expected with dealiasing. We could get rid of these?

[navidcy]

Agree -- I think we should remove the FilteredTimeSteppers. Perhaps in a different PR? But I do wanna give the users an easy way to apply filters if they want to. See FourierFlows/FourierFlows.jl#100.

[glwagner]

There a number of examples that use filtering:

• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/8b082dd0834317644a9a213cc4fd050dc83611af/examples/singlelayerqg_betadecay.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/bafe2bcd6c02afd8a693abd404a4cbd7da7d0305/examples/singlelayerqg_decaying_topography.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/bafe2bcd6c02afd8a693abd404a4cbd7da7d0305/examples/surfaceqg_decaying.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/1a4dda891359092e493fc02864a245483183e80e/examples/singlelayerqg_betaforced.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/1a4dda891359092e493fc02864a245483183e80e/examples/barotropicqgql_betaforced.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/5575c266aee329c5cc5c3a072fdc568f2ba6c5b9/examples/twodnavierstokes_decaying.jl
• https://github.com/FourierFlows/GeophysicalFlows.jl/blob/bafe2bcd6c02afd8a693abd404a4cbd7da7d0305/examples/multilayerqg_2layer.jl

So this seems like a major change to the code. It might be warranted, but should be considered carefully.

[navidcy]

I still feel that there is a disconnect in communication.

I am trying to say that regardless of filtering, hyper viscosity, and what not, dealiasing should be -- in principle -- done. Not doing dealiasing brings errors from aliasing. Oftentimes we can get away with thous with filtering or viscosity. But even in cases where filtering/viscosity seems to work "OK", I can't see why doing dealiasing on top of everything else shouldn't be done. In fact, I'm trying to argue that this is the proper way.

If we forget about examples, code changes, and technicalities that should come along with such a change, do you agree with the above statement?

[glwagner]

I believe that the papers which introduce filtering do not dealias.

I agree that aliasing errors are real, however. I wouldn't try to enter an argument about whether or not a numerical computation without dealiasing is "improper" or "incorrect".

Not dealiasing is essentially a performance optimization which seems to rely on an assumption that power in high wavenumbers is negligible due to explicit dissipation. It's a common approximation in many spectral calculations --- proper or not.

I'd be fine with making dealiasing default and even getting rid of all the examples that use filtering. It's a bit of work to get there and does reduce available features, so there's some cause to pause and consider.

[glwagner]

If we forget about examples, code changes, and technicalities that should come along with such a change, do you agree with the above statement?

No, I don't agree that dealiasing "should" be done in all circumstances.

I think dealiasing is necessary to eliminate aliasing errors. Aliasing errors often compromise solutions with significant power at high wave numbers. This often occurs in marginally resolved simulations whose resolution is close to being as low as necessary to resolve the dynamics of a particular flow.

[navidcy]

I believe that the papers which introduce filtering do not dealias.

I agree that aliasing errors are real, however. I wouldn't try to enter an argument about whether or not a numerical computation without dealiasing is "improper" or "incorrect".

Not dealiasing is essentially a performance optimization which seems to rely on an assumption that power in high wavenumbers is negligible due to explicit dissipation. It's a common approximation in many spectral calculations --- proper or not.

I'd be fine with making dealiasing default and even getting rid of all the examples that use filtering. It's a bit of work to get there and does reduce available features, so there's some cause to pause and consider.

Oh, I hope the disagreement was not due the use of "proper" in the PR title... I can happily take that back! It's not about being proper or not. My main point was that there always exists some aliasing error and if users decide to ignore it because it doesn't seem to be important in the particular calculation or because they want the code to run faster then fair enough -- they can do it.

[navidcy]

If we forget about examples, code changes, and technicalities that should come along with such a change, do you agree with the above statement?

No, I don't agree that dealiasing "should" be done in all circumstances.

I think dealiasing is necessary to eliminate aliasing errors. Aliasing errors often compromise solutions with significant power at high wave numbers. This often occurs in marginally resolved simulations whose resolution is close to being as low as necessary to resolve the dynamics of a particular flow.

OK, I see this point. Now I'm skeptical about this whole PR...

[glwagner]

It's probably a positive change. Users looking for a speed up can set dealias=0 and interpret their resolution as the full resolution. Most people perform computations at marginal resolution. Filtering seems like an alright alternative to dealiasing, but hard to interpret. The conservative approach is to dealias so its a good default.

I think dealiasing by default is clearer though using the 3/2 rule, because people won't get confused about the resolution their running at.

[navidcy]

Yes, but that requires padding, etc... Also, I never understood how FFTs are optimised then because e.g. 128*3/2 = 192 and 192 is not a power of 2....

I know some argue for double padding, e.g. 128->256! But that seems a lot of waste, in particular when you think of 3D sols...

[navidcy]

An alternative to this PR would be to demonstrate dealiasing in an example. So that could, e.g.,

dealias!(sol, grid)
stepforward!(prob)

But question that arises (in my head at least) is: Is this enough for multi-step time-stepping? E.g., RK4 involves 4 calls of calcN! but dealiasing only happens once.

[glwagner]

These days FFTs are optimized for grids that are products of many possible small numbers like 2, 3, 5, 7, and possibly more.

[glwagner]

An alternative to this PR would be to demonstrate dealiasing in an example. So that could, e.g.,

dealias!(sol, grid)
stepforward!(prob)

But question that arises (in my head at least) is: Is this enough for multi-step time-stepping? E.g., RK4 involves 4 calls of calcN! but dealiasing only happens once.

I don't think it's enough. To completely eliminate aliasing errors we have to dealias after every nonlinear product I think.

[navidcy]

Then the example could demonstrate how the user can import calcN! and redefine it, e.g.

import GeophysicalFlows.TwoDNavierStokes.calcN!

function calcN!(N, sol, t, clock, vars, params, grid)
  dealias!(sol, grid)
  
  TwoDNavierStokes.calcN_advection!(N, sol, t, clock, vars, params, grid)

  TwoDNavierStokes.addforcing!(N, sol, t, clock, vars, params, grid)
      
  return nothing
end

Or this gets to much into the nitty-gritty of the code?

[glwagner]

You could set the default aliased_fraction=0 and then illustrate the effect of aliased_fraction=1/3 and explain how this implies a reduced effective resolution? If the simulation isn't well resolved then the aliasing errors might be apparent, illustrating the value of dealiasing.

[navidcy]

Let's first wait for FourierFlows/FourierFlows.jl#285 first.