Upgrading eigentools scripts

Version 2 of eigentools has made significant changes to the API and will necessitate some changes (for the better, we hope) to the user experience. The guiding principle behind the new API is that one should no longer need to touch the Dedalus EVP object that defines the eigenvalue problem at hand.

Most importantly, no changes need to be made to the underlying Dedalus EVP object.

Basic eigenproblem usage

Choosing a sparse or dense solve is no longer done when instantiating Eigenproblem objects. Instead, this is a choice at solve time:

EP = Eigenproblem(string, reject=True)
EP.solve(sparse=False)

Also, notice that rejection of spurious modes is now done automatically with EP.solve if reject=True is selected at instantiation time. Note that although in the above code, we explicitly set reject=True, this is unnecessary, as it is the default. The EP.reject_spurious() function has been removed

In addition, solving again with different parameters has been greatly simplified from the previous version. You now simply pass a dictionary with the parameters you wish to change to solve itself. Let’s revisit the simple waves-on-a-string problem from the getting started page, but add a parameter, c2, the wave speed squared.

Here, we solve twice, once with c1 = 1 and once with c2 = 2. Given the dispersion relation for this problem is \(\omega^2 = c^2 k\) and our eigenvalue omega is really \(\omega^2\), we expect the eigenvalues for the second solve to be twice those for the first.

import numpy as np
from eigentools import Eigenproblem
import dedalus.public as de

Nx = 128
x = de.Chebyshev('x',Nx, interval=(-1, 1))
d = de.Domain([x])

string = de.EVP(d, ['u','u_x'], eigenvalue='omega')
string.parameters['c2'] = 1
string.add_equation("omega*u + c2*dx(u_x) = 0")
string.add_equation("u_x - dx(u) = 0")
string.add_bc("left(u) = 0")
string.add_bc("right(u) = 0")
EP = Eigenproblem(string)
EP.solve(sparse=False)
evals_c1 = EP.evalues
EP.solve(sparse=False, parameters={'c2':2})
evals_c2 = EP.evalues

print(np.allclose(evals_c2, 2*evals_c1))

Getting eigenmodes

Getting eigenmodes has also been simplified and significantly extended. Previously, getting an eigenmode corresponding to an eigenvalue required using the set_state() method on the underlying EVP object. In keeping with the principle of not needing to manipulate the EVP, we provide a new .eigenmode(index), where index is the mode number corresponding to the eigenvalue index in EP.evalues. By default, with mode rejection on, these are the “good” eigenmodes. .eigenmode(index) returns a Dedalus FieldSystem object, with a Dedalus Field for each field in the eigenmode:

emode = EP.eigenmode(0)
print([f.name for f in emode.fields])
u = emode.fields[0]
u_x = emode.fields[1]

Finding critical parameters

This has been considerably cleaned up. The two major things to note are that

  1. one no longer needs to create a shim function to translate between an x-y grid and the parameter names within the EVP.

  2. The parameter grid is no longer defined inside CriticalFinder, but is instead created by the user and passed in

For example, here are the relevant changes necessary for the MRI test problem.

First, replace

EP = Eigenproblem(mri, sparse=True)

# create a shim function to translate (x, y) to the parameters for the eigenvalue problem:
def shim(x,y):
    iRm = 1/x
    iRe = (iRm*Pm)
    print("Rm = {}; Re = {}; Pm = {}".format(1/iRm, 1/iRe, Pm))
    gr, indx, freq = EP.growth_rate({"Q":y,"iRm":iRm,"iR":iRe})
    ret = gr+1j*freq
    return ret

cf = CriticalFinder(shim, comm)

with

EP = Eigenproblem(mri)

cf = CriticalFinder(EP, ("Q", "Rm"), comm, find_freq=False)

Important: note that find_freq is specified at instantiation rather than when calling cf.crit_finder later.

Once this is done, the grid generation changes from

mins = np.array((4.6, 0.5))
maxs = np.array((5.5, 1.5))
ns   = np.array((10,10))
logs = np.array((False, False))

cf.grid_generator(mins, maxs, ns, logs=logs)

to

nx = 20
ny = 20
xpoints = np.linspace(0.5, 1.5, nx)
ypoints = np.linspace(4.6, 5.5, ny)

cf.grid_generator((xpoints, ypoints), sparse=True)