NSE Equations Vorticity evolution integer overflow

[Ish Kaul]

Hi All!

Apologies if this is a basic question, since I am new to dedalus. I was implementing the vorticity formulation of the Navier-Stokes equation, but I had two problems: I wanted to have periodic boundary conditions but at the same time handle the zero of the frequencies when inverting the poisson equation - this gives me an error saying the LHS operator interp_y is coupled along direction y. Second: I have a bunch of nan's when I actually evolve the vorticity in time in the main loop and so it doesn't return any result. I have copied my code here for reference:

Lx, Ly = (2., 2.)

nx, ny = (152, 152)

# Create bases and domain

x_basis = de.Fourier('x', nx, interval=(-Lx/2, Lx/2), dealias=3/2)

y_basis = de.Fourier('y',ny, interval=(-Ly/2, Ly/2), dealias=3/2)

domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)

viscosity = 0

problem = de.IVP(domain, variables = ['omega', 'psi', 'omegay'])

problem.parameters['nu'] = viscosity

problem.add_equation("dt(omega) = omegay*dx(psi) - dx(omega)*dy(psi) + nu*dx(dx(omega)) + nu*dy(dy(omega))")

problem.add_equation("-omega + 0.01*psi= dx(dx(psi)) + dy(dy(psi))")

problem.add_equation("omegay - dy(omega) = 0")

# This next line when uncommented gives the first error

# problem.add_bc("left(psi) = 0", condition = "(ny==0)")

ts = de.timesteppers.RK443

solver = problem.build_solver(ts)

x = domain.grid(0)

y = domain.grid(1)

omega = solver.state['omega']

psi = solver.state['psi']

omegay = solver.state['omegay']

# Intializing sigma and mu

sigma = 1

mu1 = 0.25

mu2=-0.25

# Initializing Gaussian monopoles on omega

dst1 = (x-mu1)**2 + (y)**2

dst2 = (x-mu2)**2 + (y)**2

omega['g'] = np.exp(-(dst1/(sigma**2 /20))) + np.exp(-(dst2 / (sigma**2 /20)))

omega.differentiate('y', out = omegay)

solver.stop_sim_time = 5.01

solver.stop_wall_time = np.inf

solver.stop_iteration = np.inf

dt = 0.01

start_time = time.time()

while solver.ok:

   print(omega['g'])

   solver.step(dt)

   if solver.iteration % 10 == 0:

      print(omega['g'])

end_time = time.time()

Here omega becomes an array of nan's with the following message: 

/usr/local/lib/python3.7/dist-packages/dedalus/core/operators.py:791: RuntimeWarning: overflow encountered in multiply np.multiply(arg0.data, arg1.data, out.data) 

/usr/local/lib/python3.7/dist-packages/dedalus/core/operators.py:510: RuntimeWarning: invalid value encountered in add np.add(arg0.data, arg1.data, out.data)

I am not sure what is going wrong and I would appreciate any help!

Thank you,

Ish

[Fryderyk Wilczynski]

Hi Ish, 

I think there are 2 things that are causing problems here.

First and foremost, unfortunately, spectral method require finite diffusion - you need to set your viscosity parameter > 0 (as small as you like, provided your resolution is fine enough).
Second, I think you need to impose conditions for the equations — you want equations to apply to every Fourier mode except the mean mode (where those equations reduce to 0 = 0), so the condition should be “(nx != 0) or (ny != 0)”,

and you then need to add an additional equation for the mean mode “(nx == 0) and (ny == 0)”. 

I have modified your code with these changes and the simulation runs without hiccups (I am attaching the script and example output).

Best wishes, 

Fryderyk

[Keaton Burns]

This looks great -- I’ll just add that the first-order formulation is only necessary for problems with a Chebyshev basis.  For a full Fourier problem like this, you can instead use “psi" as your only variable, and define “u”, “w”, and “omega” via substitutions instead of additional equations.  This should speed things up a lot!  

I’ll also emphasize that the other important change Fryderyk made, versus the original code, was that all of the linear terms were moved to the LHS of the equations.  This is necessary for them to be implicitly integrated/coupled.

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/1db507e4-a027-4971-8866-05ecb013785an%40googlegroups.com.

[Jack]

Hi Fryderyk, 

Thanks for sharing the scripts. I'm quite interested in your mp4 file. What tool do you use to make the movie? 

Jack