Smaller L in twodturb_stochasticforcing.jl does not work

[natschil]

Setting L = 1\pi in twodturb_stochasticforcing.jl gives a NaN solution, is this a bug somewhere?

[glwagner]

It's not a bug, but a limitation in the example. The value of dt is not independent from L. I was able to get the example to run without blowing up with L=π by setting dt=0.001 on line 12.

Let me know if that works for you.

[glwagner]

As a side note, NaNs typically result from either the time-step being too large, or the viscous/diffusive terms being too small.

[natschil]

ah I see, thanks. This is a bit suprising to me as I would expect a smaller domain to not need a smaller timestep, but perhaps it is related to the forcing somehow?

Btw, is there some kind of scaling of the velocity fields vs.u and vs.v ? I'm trying to advect some particles by them but I'm getting very small velocity values which don't seem to be consistent with the movement of vortices in the vorticity field.

[navidcy]

@natschil, with a quick glance at your question let me say the following:

Indeed @glwagner suggestion would "heal" your problem. But, as you mention, thinking of the CFL criterion it's very counterintuitive that you need to reduce dt for smaller domain.

I believe that the problem you are pointing out stems from a different numerical instability that has nothing to do with CFL. Instead, it seems you've hit a numerical instability that comes from the hyper viscous term -\nu*𝛁^4.

There is a different numerical stability criterion (other than CLF) that has to be satisfied for the hyper viscous term. For a hyper-viscosity of nnu=2 as the one used in the example, this reads:

nu < nucritical = 2.79 / (k_max^4*dt)

For the default values of the example: L = 2π, dt=0.005, n=256, we have that k_max=maximum(sqrt(g.Krsq))=181 and therefore:

nucritical = 5.2e-7

When you change domain size to L=π and keep everything the same you have k_max=362 which implies

nucritical = 3.2e-8

so the value nu=1e-7>nucritical and leads to instability.

You can "heal" this by either reducing nu. Otherwise you can reducing dt to, e.g., 0.001 (as @glwagner suggests) or reducing n to, e.g., 128 so that the nucritical becomes larger than 1e-7.

(Chris Bretherton has a brief discussion on the subject here: https://atmos.uw.edu/academics/classes/2012Q2/581/lect/lect21.pdf but if you come across a more detailed reference please let me know.)

[navidcy]

So it turns out in this case that it's not that "the viscous/diffusive terms being too small" (as @glwagner points out) but in fact that they were very big!

Although my first reaction would be exactly what @glwagner suggested, I think there is a bit more to it. I agree that this is very counter-intuitive.

[navidcy]

@natschil: Regarding you question about scaling. No there is not any scaling imposed on velocities. Can you perhaps show me what you are trying to do and point out where the problem is through an example script? I'll try to dig into it.

[navidcy]

@natschil: A way around these numerical instabilities related to hyper viscosity is to use semi-implicit time stepper ETDRK4 which solves for the linear terms (that include the hyper viscosity) exactly.

See:

Cox, S. M., & Matthews, P. C. (2002). Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2), 430-455.

Kassam, A. K., & Trefethen, L. N. (2005). Fourth-order time-stepping for stiff PDEs. SIAM Journal on Scientific Computing, 26(4), 1214-1233.

or some notes I've written up (which explain how we implement this time stepper in FourierFlows.jl and can be found here: https://github.com/navidcy/ETDRK4_notes

[natschil]

@navidcy thanks for the detailed response! I will try using ETDRK4, I didn't try it before as most of the examples seemed to use some form of RK4. Actually I'm trying to solve Navier-Stokes for the case without hyperviscosity, am I correct that I get this by specifying nnu=1 and not specifying mu and nmu?

For some context as to what I am trying to do: I would like to have a turbulent 2d flow in which has different vortices at various times in order to test some methods for computing Lagrangian coherent structures (cf https://coherentstructures.github.io/CoherentStructures.jl/latest/). For this I need a velocity field with which I advect particles, what I do is that I save the velocity field at various points and then interpolate and solve an ode to get trajectories (there are of course better ways to do this).

If I do this I get something like:

https://gist.github.com/natschil/2b4713b4b0d821c5bc84d07f20aeb07d

Here I get maximum(abs.(all_us)) = 0.0134 (and similar for v), but a the vortices seem to move quite a bit than that. So I was thinking that maybe u and v are somehow scaled, but if this is not the case maybe I'm doing something else wrong.

PS: What kind of dealiasing do you use ?

[navidcy]

1. Without hyper viscosity but with normal viscosity? Then yes, nnu=1 would give you normal Laplacian viscosity with coefficient nu. Also you should set mu=0 to remove linear drag.
   Mind you that with just viscosity + stochastic forcing it will take a really long time for the flow to equilibrate. Characteristically, if you expect to have vortices with scales 1/k_vort you will need to run the simulation for ~5*1/(\nu*k_vort^2). (See the discussion surrounding fig. 7 here.)

2. Regarding you saving the velocity fields: Remember that flow fields are solved in Fourier space and the only dynamic variable that is evolved is, in this case, the fourier transform of the relative vorticity sol=qh. Thus, you need to call TwoDTurb.updatevars!(prob) before saving the flow fields (i.e. just before line 77 in your .jl script) instead of at line 80. This way v.u and v.v are updated before they are saved.

I believe this is where your problems stem from. Let me know if that works.

3. There is no dealiasing by default. You can specifically call dealias!(sol, grid) function after each time-step. (We are planning to add an easier way to do that giving dealias! as an argument to the time-stepper constructor but we haven't got around doing that yet.)
   The dealias! function will zero out the top 1/3-rd higher wavenumbers: see https://github.com/FourierFlows/FourierFlows.jl/blob/master/src/domains.jl for details.

Let me know if these clear out your issues and if there is anything else you want.

[glwagner]

@navidcy @natschil it would also be possible to "use dealiasing" by adding dealias!(sol, g) at the end of the forcing function (between lines 36 and 37 here):

https://gist.github.com/natschil/2b4713b4b0d821c5bc84d07f20aeb07d

This is more correct because you dealias between substeps (which could matter for multistep schemes like RK4).

Note that the dealiasing functions are part of FourierFlows:

https://github.com/FourierFlows/FourierFlows.jl/blob/10281209059e8378a9c0838652312eacea903cc4/src/domains.jl#L171

It would be pretty interesting to write a new module for GeophysicalFlows that advects passive particles while solving a 2D turbulence problem, hooking into the FourierFlows timesteppers. I don't think that would be very difficult (famous last words). Though if I understood correctly you are actually specifically interested in interpolating the velocity fields (in a fancy way?) to find tracer trajectories, because that is actually the topic you are researching?

[natschil]

@navidcy :

Without hyper viscosity but with normal viscosity?

Yes, with normal viscosity. Thanks for the comment about removing drag, I was doing it wrong until now. Maybe the additional drag is why my particles move too slowly -- I'm guessing that in Navier-Stokes+drag the vorticity no longer obeys an advection-diffusion equation like it does in the ordinary Navier-Stokes equation, so I shouldn't expect the vorticity to be advected-diffused by the flow. Since I want the vortices to be as Lagrangian as possible, it probably makes sense to just bite the bullet and run things for a long time until they are more or less statistically stable (unless the stochastic forcing part is what makes it slow, and there is an implementation of deterministic forcing I can use somewhere and this is faster). Note that if you have any better suggestions as to how to get Lagrangian vortices of reasonable size with e.g. hyper/o-viscosity I'm all ears :)

The dealias! function will zero out the top 1/3-rd higher wavenumbers: see
Have you thought about using 3/2-dealiasing? here you don't lose the highest wavenumbers, but you pay for it with a larger Fourier transform.

you need to call TwoDTurb.updatevars!(prob) before saving the flow fields (i.e. just before line 77 in your .jl script) instead of at line 80. This way v.u and v.v are updated before they are saved.

Thanks! I will try this tomorrow.

@glwagner: Though if I understood correctly you are actually specifically interested in interpolating the velocity fields (in a fancy way?) to find tracer trajectories, because that is actually the topic you are researching?

No, but I need to seed particles (or more accurately a finite-difference stencil of particles) at different points of space and time, and the easiest way to do this is to save the velocity fields and interpolate.

[navidcy]

@natschil, I have thoughts on how you could get a mature equilibrated 2D turbulent flow. I don't really understand why you object to linear drag and want to have an advection-diffusion equation. But if that's what you need then yes, plain viscosity is the way to go.

I think though that this discussion starts drifting away from being an issue of the package. I am more than happy to resume the discussion over email (I'm sure @glwagner would be happy to join). What do you think?

Before doing that though let's sort out the issue. Let me know if what I suggested regarding (i) hyper viscosity coefficient value or using ETDRK4 and (ii) calling TwoDTurb.updatevars!(prob) before saving u and v have resolved your issues. If so, we'll close this issue and continue the discussion over email.

[natschil]

calling TwoDTurb.updatevars!(prob) before saving u and v have resolved your issues

Yes, that seems to have fixed it! Thanks for the help!