# Kelvin-Helmholtz Instability (image: Chuck Doswell)

We will simulate the incompressible Kelvin-Helmholtz problem. We non-dimensionalize the problem by taking the box height to be one and the jump in velocity to be one. Then the Reynolds number is given by

$\mathrm{Re} = \frac{U H}{\nu} = \frac{1}{\nu}.$

We use no slip boundary conditions, and a box with aspect ratio $$L/H=2$$. The initial velocity profile is given by a hyperbolic tangent, and only a single mode is initially excited. We will also track a passive scalar which will help us visualize the instability.

First, we import the necessary modules.

:

%matplotlib inline

:

import numpy as np
import matplotlib.pyplot as plt
import h5py
import time
from IPython import display

from dedalus import public as de
from dedalus.extras import flow_tools

import logging
logger = logging.getLogger(__name__)


To perform an initial value problem (IVP) in Dedalus, you need three things:

1. A domain to solve the problem on

2. Equations to solve

3. A timestepping scheme

## Problem Domain

First, we will specify the domain. Domains are built by taking the direct product of bases. Here we are running a 2D simulation, so we will define $$x$$ and $$y$$ bases. From these, we build the domain.

:

# Aspect ratio 2
Lx, Ly = (2., 1.)
nx, ny = (192, 96)

# Create bases and domain
x_basis = de.Fourier('x', nx, interval=(0, Lx), dealias=3/2)
y_basis = de.Chebyshev('y',ny, interval=(-Ly/2, Ly/2), dealias=3/2)
domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)


The last basis ($$y$$ direction) is represented in Chebyshev polynomials. This will allow us to apply interesting boundary conditions in the $$y$$ direction. We call the other directions (in this case just $$x$$) the “horizontal” directions. The horizontal directions must be “easy” in the sense that taking derivatives cannot couple different horizontal modes. Right now, we have Fourier and Sin/Cos series implemented for the horizontal directions, and are working on implementing spherical harmonics.

## Equations

Next we will define the equations that will be solved on this domain. The equations are

$\partial_t u + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} u + \partial_x p = \frac{1}{{\rm Re}} \nabla^2 u$
$\partial_t v + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} v + \partial_y p = \frac{1}{{\rm Re}} \nabla^2 v$
$\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0$
$\partial_t s + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} s = \frac{1}{{\rm Re}{\rm Sc}} \nabla^2 s$

The equations are written such that the left-hand side (LHS) is treated implicitly, and the right-hand side (RHS) is treated explicitly. The LHS is limited to only linear terms, though linear terms can also be placed on the RHS. Since $$y$$ is our special direction in this example, we also restrict the LHS to be at most first order in derivatives with respect to $$y$$.

We also set the parameters, the Reynolds and Schmidt numbers.

:

Reynolds = 1e4
Schmidt = 1.

problem = de.IVP(domain, variables=['p','u','v','uy','vy','S','Sy'])
problem.parameters['Re'] = Reynolds
problem.parameters['Sc'] = Schmidt
problem.add_equation("dt(u) + dx(p) - 1/Re*(dx(dx(u)) + dy(uy)) = - u*dx(u) - v*uy")
problem.add_equation("dt(v) + dy(p) - 1/Re*(dx(dx(v)) + dy(vy)) = - u*dx(v) - v*vy")
problem.add_equation("dt(S) - 1/(Re*Sc)*(dx(dx(S)) + dy(Sy)) = - u*dx(S) - v*Sy")


Because we are using this first-order formalism, we define auxiliary variables uy, vy, and Sy to be the $$y$$-derivative of u, v, and S respectively.

Next, we set our boundary conditions. “Left” boundary conditions are applied at $$y=-Ly/2$$ and “right” boundary conditions are applied at $$y=+Ly/2$$.

:

problem.add_bc("left(u) = 0.5")
problem.add_bc("right(v) = 0", condition="(nx != 0)")
problem.add_bc("left(p) = 0", condition="(nx == 0)")


Note that we have a special boundary condition for the $$k_x=0$$ mode (singled out by condition="(nx==0)"). This is because the continuity equation implies $$\partial_y v=0$$ if $$k_x=0$$; thus, $$v=0$$ on the top and bottom are redundant boundary conditions. We replace one of these with a gauge choice for the pressure.

## Timestepping

We have implemented a variety of multi-step and Runge-Kutta implicit-explicit timesteppers in Dedalus. The available options can be seen in the timesteppers.py module. For this problem, we will use a third-order, four-stage Runge-Kutta integrator. Changing the timestepping algorithm is as easy as changing one line of code.

:

ts = de.timesteppers.RK443


## Initial Value Problem

We now have the three ingredients necessary to set up our IVP:

:

solver =  problem.build_solver(ts)

2022-08-26 11:24:29,478 pencil 0/1 INFO :: Building pencil matrix 1/96 (~1%) Elapsed: 0s, Remaining: 1s, Rate: 7.3e+01/s
2022-08-26 11:24:29,590 pencil 0/1 INFO :: Building pencil matrix 10/96 (~10%) Elapsed: 0s, Remaining: 1s, Rate: 8.0e+01/s
2022-08-26 11:24:29,711 pencil 0/1 INFO :: Building pencil matrix 20/96 (~21%) Elapsed: 0s, Remaining: 1s, Rate: 8.1e+01/s
2022-08-26 11:24:29,833 pencil 0/1 INFO :: Building pencil matrix 30/96 (~31%) Elapsed: 0s, Remaining: 1s, Rate: 8.1e+01/s
2022-08-26 11:24:29,963 pencil 0/1 INFO :: Building pencil matrix 40/96 (~42%) Elapsed: 0s, Remaining: 1s, Rate: 8.0e+01/s
2022-08-26 11:24:30,090 pencil 0/1 INFO :: Building pencil matrix 50/96 (~52%) Elapsed: 1s, Remaining: 1s, Rate: 8.0e+01/s
2022-08-26 11:24:30,211 pencil 0/1 INFO :: Building pencil matrix 60/96 (~62%) Elapsed: 1s, Remaining: 0s, Rate: 8.0e+01/s
2022-08-26 11:24:30,335 pencil 0/1 INFO :: Building pencil matrix 70/96 (~73%) Elapsed: 1s, Remaining: 0s, Rate: 8.0e+01/s
2022-08-26 11:24:30,474 pencil 0/1 INFO :: Building pencil matrix 80/96 (~83%) Elapsed: 1s, Remaining: 0s, Rate: 7.9e+01/s
2022-08-26 11:24:30,606 pencil 0/1 INFO :: Building pencil matrix 90/96 (~94%) Elapsed: 1s, Remaining: 0s, Rate: 7.9e+01/s
2022-08-26 11:24:30,686 pencil 0/1 INFO :: Building pencil matrix 96/96 (~100%) Elapsed: 1s, Remaining: 0s, Rate: 7.9e+01/s


Now we set our initial conditions. We set the horizontal velocity and scalar field to tanh profiles, and using a single-mode initial perturbation in $$v$$.

:

x = domain.grid(0)
y = domain.grid(1)
u = solver.state['u']
uy = solver.state['uy']
v = solver.state['v']
vy = solver.state['vy']
S = solver.state['S']
Sy = solver.state['Sy']

a = 0.05
sigma = 0.2
flow = -0.5
amp = -0.2
u['g'] = flow*np.tanh(y/a)
v['g'] = amp*np.sin(2.0*np.pi*x/Lx)*np.exp(-(y*y)/(sigma*sigma))
S['g'] = 0.5*(1+np.tanh(y/a))
u.differentiate('y',out=uy)
v.differentiate('y',out=vy)
S.differentiate('y',out=Sy)

:

<Field 5147257408>


Now we set integration parameters and the CFL.

:

solver.stop_sim_time = 2.01
solver.stop_wall_time = np.inf
solver.stop_iteration = np.inf

initial_dt = 0.2*Lx/nx
cfl = flow_tools.CFL(solver,initial_dt,safety=0.8)


## Analysis

We have a sophisticated analysis framework in which the user specifies analysis tasks as strings. Users can output full data cubes, slices, volume averages, and more. Here we will only output a few 2D slices, and a 1D profile of the horizontally averaged concentration field. Data is output in the hdf5 file format.

:

analysis = solver.evaluator.add_file_handler('analysis_tasks', sim_dt=0.1, max_writes=50)
solver.evaluator.vars['Lx'] = Lx


## Main Loop

We now have everything set up for our simulation. In Dedalus, the user writes their own main loop.

:

# Make plot of scalar field
x = domain.grid(0,scales=domain.dealias)
y = domain.grid(1,scales=domain.dealias)
xm, ym = np.meshgrid(x,y)
fig, axis = plt.subplots(figsize=(10,5))
p = axis.pcolormesh(xm, ym, S['g'].T, cmap='RdBu_r');
axis.set_xlim([0,2.])
axis.set_ylim([-0.5,0.5])

logger.info('Starting loop')
start_time = time.time()
while solver.ok:
dt = cfl.compute_dt()
solver.step(dt)
if solver.iteration % 10 == 0:
# Update plot of scalar field
p.set_array(np.ravel(S['g'].T))
display.clear_output()
display.display(plt.gcf())
logger.info('Iteration: %i, Time: %e, dt: %e' %(solver.iteration, solver.sim_time, dt))

end_time = time.time()

p.set_array(np.ravel(S['g'].T))
display.clear_output()
# Print statistics
logger.info('Run time: %f' %(end_time-start_time))
logger.info('Iterations: %i' %solver.iteration)

2022-08-26 11:25:50,227 __main__ 0/1 INFO :: Run time: 78.498186
2022-08-26 11:25:50,227 __main__ 0/1 INFO :: Iterations: 269 ## Analysis

As an example of doing some analysis, we will load in the horizontally averaged profiles of the scalar field $$s$$ and plot them.

:

# Read in the data
y = f['/scales/y/1.0'][:]
t = f['scales']['sim_time'][:]
f.close()

S_ave = S_ave[:,0,:] # remove length-one x dimension

:

for i in range(0,21,5):
plt.plot(S_ave[i,:],y,label='t=%4.2f' %t[i])

plt.ylim([-0.5,0.5])
plt.xlim([0,1])
plt.xlabel(r'$\frac{\int \ S dx}{L_x}$',fontsize=24)
plt.ylabel(r'$y$',fontsize=24)
plt.legend(loc='lower right').draw_frame(False) [ ]: