EVPs: set_state() in multi-d problems runs into Nyquist mode zeroing

[bpbrown]

For multi-dimensional Eigenvalue Problems solved using our current best-practices approach (see examples/evp_1d_rayleigh_benard/rayleigh_benard_evp.py) we run into an edge case with set_modes() when trying to plot eigenmodes.

Best practices are to set up a small dimension in the wave-perturbation direction of size Nx=2, and length Lx = 2 * np.pi / kx, which puts the mode to be studied in the second index. When we do this, and then use set_state, the associated fields have values in coefficient space, but the first transform to grid space zeros out all of the data.

Here is a minimal example that reproduces this, using Rayleigh-Benard convection with stress-free boundaries (where the modes are analytic).:

import numpy as np
import dedalus.public as d3

# Parameters
Nz = 64
Rayleigh = 27/4*np.pi**4*(1+1e-3)
Prandtl = 1
kx = kx_crit = np.pi/np.sqrt(2)

# Parameters
Lz = 1
# Build Fourier basis for x with prescribed kx as the fundamental mode
Nx = 2 # change to Nx = 4 patches around Nyquist
Lx = 2 * np.pi / kx

# Bases
coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.complex128)
xbasis = d3.ComplexFourier(coords['x'], size=Nx, bounds=(0, Lx))
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz))

# Fields
omega = dist.Field(name='omega')
p = dist.Field(name='p', bases=(xbasis,zbasis))
b = dist.Field(name='b', bases=(xbasis,zbasis))
u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
tau_p = dist.Field(name='tau_p')
tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)

# Substitutions
kappa = (Rayleigh * Prandtl)**(-1/2)
nu = (Rayleigh / Prandtl)**(-1/2)

ex, ez = coords.unit_vector_fields(dist)
lift_basis = zbasis.derivative_basis(1)
lift = lambda A: d3.Lift(A, lift_basis, -1)
grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction
grad_b = d3.grad(b) + ez*lift(tau_b1) # First-order reduction
dt = lambda A: -1j*omega*A
e = grad_u + d3.trans(grad_u)

# Problem
problem = d3.EVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals(), eigenvalue=omega)
problem.add_equation("trace(grad_u) + tau_p = 0")
problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) - ez@u = 0")
problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = 0")
problem.add_equation("b(z=0) = 0")
problem.add_equation('ez@u(z=0) = 0')
problem.add_equation('ez@(ex@e(z=0)) = 0')
problem.add_equation("b(z=Lz) = 0")
problem.add_equation('ez@u(z=Lz) = 0')
problem.add_equation('ez@(ex@e(z=Lz)) = 0')
problem.add_equation("integ(p) = 0") # Pressure gauge

# Solver
solver = problem.build_solver(entry_cutoff=0)
solver.solve_sparse(solver.subproblems[1], 10, target=0)
i = np.argsort(solver.eigenvalues.imag)[-1]
print('omega = {:}'.format(solver.eigenvalues[i]))
solver.set_state(i, 0)
print('b[c] max amp = {:}'.format(np.max(np.abs(b['c']))))
print('b[g] max amp = {:}'.format(np.max(np.abs(b['g']))))
print('b[c] max amp = {:}'.format(np.max(np.abs(b['c']))))

The result of this script is:

omega = (3.924211175128554e-20+0.0002884588085073608j)
b[c] max amp = 0.48998225974907084
b[g] max amp = 0.0
b[c] max amp = 0.0

The correct mode is found (an unstable growing mode), but the eigenfunction is zero when evaluated in grid space and remains that way on return to coeff space.

This is related to zeroing of Nyquist modes.

If Nx = 4 is used instead of Nx = 2, then correct behaviour is instead found:

omega = (3.924211175128554e-20+0.0002884588085073608j)
b[c] max amp = 0.48998225974907084
b[g] max amp = 0.5852474523582276
b[c] max amp = 0.4899822597490708

Now the eigenmodes are not zeroed on transform, retain the same coeff amplitudes on return to coeff space, and match both eigenvalue and coeff amplitude with the original implementation. The individual eigenmode structures now also match expectations.

[kburns]

It's unclear what we want to do to change things here. We could stop zeroing the Nyquist mode for ComplexFourier, which is generally what is used with EVPs. But I think we do want to continue dropping the Nyquist mode for RealFourier since otherwise there will be weird matrix shape errors since only the cosine mode is supported there. That would mean real and complex Fourier have different behavior... but maybe that isn't so bad?