Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
bpbrown opened this issue Aug 31, 2023 · 1 comment
Open
Assignees
Labels
bug Something isn't working

Comments

@bpbrown
Copy link
Contributor

bpbrown commented Aug 31, 2023

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.

@bpbrown bpbrown added the bug Something isn't working label Aug 31, 2023
@kburns
Copy link
Member

kburns commented Jan 8, 2024

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

3 participants