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
one no longer needs to create a shim function to translate between an x-y grid and the parameter names within the
EVP
.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)