Eady promblem(EVP): all eigenvalues are real (no growth)

[HXC303]

Hello everyone,

I'm trying to solve the classical Eady baroclinic instability problem using Dedalus v3 (EVP) as described in Vallis (2006), Atmospheric and Oceanic Fluid Dynamics, Chapter 6., but all the eigenvalues I obtain are real — no imaginary parts, hence no growth rate.

Any suggestions or code references would be greatly appreciated!
Thanks in advance.

Below is a simplified version of my script:

`# Parameters
H = 1
Nz = 100
Ld = 10000
k_list = np.arange(0, 4/Ld, 0.01/Ld)
dtype = np.complex128

Bases

zcoord = d3.Coordinate('z')
dist = d3.Distributor(zcoord, dtype=dtype)
zbasis = d3.Chebyshev(zcoord, size=Nz, bounds=(0, H), dealias=3/2)

Fields

psi = dist.Field(name='psi', bases=zbasis)
tau_1 = dist.Field(name='tau_1')
tau_2 = dist.Field(name='tau_2')
c = dist.Field(name='c')

Substitutions

z = dist.local_grid(zbasis)
z0 = dist.Field(bases=zbasis)
z0['g'] = z
miu = dist.Field()
dz = lambda A: d3.Differentiate(A, zcoord)
lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
psiz = dz(psi) + lift(tau_1) # First-order reduction
psizz = dz(psiz) + lift(tau_2)

Problem

problem = d3.EVP([psi, tau_1, tau_2], eigenvalue=c, namespace=locals())
problem.add_equation("psizz - miu**2psi = 0")
problem.add_equation("cpsiz(z=0) + psi(z=0) = 0") # z = 0
problem.add_equation("(c-1)*psiz(z=H) + psi(z=H) = 0") # z = H

Solver

solver = problem.build_solver()
evals_all = np.zeros((len(k_list), Nz+2))
for i, ki in enumerate(k_list):
miu['g'] = ki * Ld
solver.solve_dense(solver.subproblems[0], rebuild_matrices=True)
evals_i = np.sort(solver.eigenvalues)
evals_all[i, :] = evals_i

cr = evals_all.real
ci = evals_all.imag
`
This is my results：

The following is the correct result form (Vallis, 2006) fig6.10: