Fix typo in OS example

docs/notebooks/Mixed Layer Instability.ipynb

@@ -64,10 +64,10 @@
    ],
    "source": [
     "# Only need the z direction to compute eigenvalues\n",
-    "z = de.Chebyshev('z',nz, interval=[-1,0])\n",
+    "z = de.Chebyshev('z', nz, interval=[-1, 0])\n",
     "d = de.Domain([z], grid_dtype=np.complex128)\n",
     "\n",
-    "evp = de.EVP(d, ['u','v','w','p','b'],'sigma')\n",
+    "evp = de.EVP(d, ['u', 'v', 'w', 'p', 'b'], 'sigma')\n",
     "\n",
     "evp.parameters['Ri'] = Ri\n",
     "evp.parameters['δ'] = δ\n",
@@ -101,7 +101,7 @@
     "\n",
     "In this problem the time dependence is given by\n",
     "\n",
-    "$$e^{i \\sigma t},$$\n",
+    "$$e^{i \\sigma t}$$,\n",
     "\n",
     "so the growth rate is $-Im(\\sigma)$ and the frequency is $Re(\\sigma)$."
    ]
@@ -146,7 +146,7 @@
     }
    ],
    "source": [
-    "rate, indx, freq = EP.growth_rate(parameters={'Hy':Hy, 'kx': kx}, sparse=False)\n",
+    "rate, indx, freq = EP.growth_rate(parameters={'Hy' :Hy, 'kx': kx}, sparse=False)\n",
     "print(\"fastest growing mode: {} @ freq {}\".format(rate, freq))"
    ]
   },
@@ -209,8 +209,8 @@
    "outputs": [],
    "source": [
     "# define a 2D domain to plot the eigenmode\n",
-    "x = de.Fourier('x',nx)\n",
-    "d_2d = de.Domain([x,z], grid_dtype=np.float64)"
+    "x = de.Fourier('x', nx)\n",
+    "d_2d = de.Domain([x, z], grid_dtype=np.float64)"
    ]
   },
   {
@@ -272,30 +272,30 @@
     "mfig = plot_tools.MultiFigure(nrows, ncols, image, pad, margin, scale)\n",
     "fig = mfig.figure\n",
     "\n",
-    "for var in ['b','u','v','w','p']:\n",
+    "for var in ['b', 'u', 'v', 'w', 'p']:\n",
     "    fs[var]['g']\n",
     "\n",
     "# b\n",
-    "axes = mfig.add_axes(0,0, [0,0,1,1])\n",
+    "axes = mfig.add_axes(0, 0, [0,0,1,1])\n",
     "plot_tools.plot_bot_2d(fs['b'], axes=axes)\n",
     "\n",
     "# u,w\n",
-    "axes = mfig.add_axes(1,0, [0,0,1,1])\n",
+    "axes = mfig.add_axes(1, 0, [0, 0, 1, 1])\n",
     "data_slices = (slice(None), slice(None))\n",
     "xx,zz = d_2d.grids()\n",
-    "xx,zz = np.meshgrid(xx,zz)\n",
-    "axes.quiver(xx[::2,::2],zz[::2,::2], fs['u']['g'][::2,::2].T, fs['w']['g'][::2,::2].T,zorder=10)\n",
+    "xx,zz = np.meshgrid(xx, zz)\n",
+    "axes.quiver(xx[::2, ::2], zz[::2, ::2], fs['u']['g'][::2, ::2].T, fs['w']['g'][::2, ::2].T, zorder=10)\n",
     "axes.set_xlabel('x')\n",
     "axes.set_ylabel('z')\n",
     "axes.set_title('u-w vectors')\n",
     "\n",
     "# v\n",
-    "axes = mfig.add_axes(2,0, [0,0,1,1])\n",
+    "axes = mfig.add_axes(2, 0, [0, 0, 1, 1])\n",
     "plot_tools.plot_bot_2d(fs['v'], axes=axes)\n",
     "\n",
     "# eta\n",
-    "axes = mfig.add_axes(3,0, [0,0,1,1])\n",
-    "axes.plot(xx[0,:],-fs['p']['g'][:,0])\n",
+    "axes = mfig.add_axes(3, 0, [0, 0, 1, 1])\n",
+    "axes.plot(xx[0, :], -fs['p']['g'][:, 0])\n",
     "axes.set_xlabel('x')\n",
     "axes.set_ylabel(r'$\\eta$')\n",
     "\n",
@@ -507,9 +507,9 @@
     "pax.collections[0].set_clim(0,0.2)\n",
     "cax.xaxis.set_ticks_position('top')\n",
     "cax.xaxis.set_label_position('top')\n",
-    "contours = np.linspace(0,0.2,20,endpoint=False)\n",
+    "contours = np.linspace(0, 0.2, 20, endpoint=False)\n",
     "pax.contour(cf.parameter_grids[0], cf.parameter_grids[1],cf.evalue_grid.real, levels=contours, colors='k')\n",
-    "pax.figure.savefig('mixed_layer_growth_rates.png',dpi=300)"
+    "pax.figure.savefig('mixed_layer_growth_rates.png', dpi=300)"
    ]
   },
   {

docs/notebooks/Orr Somerfeld pseudospectra.ipynb

@@ -26,11 +26,11 @@
     "\n",
     "This is a very interesting problem because the linear operator for the Orr-Sommerfeld problem $\\mathcal{L}_{OS}$ is non-normal, that is $\\mathcal{L}_{OS}^\\dagger \\mathcal{L}_{OS} \\neq \\mathcal{L}_{OS}\\mathcal{L}_{OS}^\\dagger$. Non-normal operators can drive transient growth: even systems in which *all* eigenmodes are linearly stable and will decay to zero, arbitrary initial conditions can be amplified by large factors.\n",
     "\n",
-    "Here, we assume perturbations take the form $u(x,z,t) = \\hat{u} e^{i\\alpha (x - ct)}$, where $c = \\omega/\\alpha$ is the wave speed. One must take care to note that $c = c_r + i c_i$. If a perturbation has $c_i > 0$, it will linearly grow. This occurs at a Reynolds number $\\mathrm{Re} \\simeq 5772$, a result you can also find using `CriticalFinder` (it's actually a test problem for `eigentools`).\n",
+    "Here, we assume perturbations take the form $u(x,z,t) = \\hat{u} e^{i\\alpha (x - ct)}$, where $c = \\omega/\\alpha$ is the wave speed. One must take care to note that $c = c_r + i c_i$. If a perturbation has $c_i > 0$, it will grow exponentially. This occurs at Reynolds number $\\mathrm{Re} \\simeq 5772$, a result you can also find using `CriticalFinder` (it's actually a test problem for `eigentools`).\n",
     "\n",
     "However, as we will see below, even at $\\mathrm{Re} = 10000$, the growth rate is tiny: $c_i \\approx 3.7\\times10^{-3}$! Yet at those Reynolds numbers, channel flow is clearly unstable. How can this be?\n",
     "\n",
-    "The answer is that the *spectrum* of the Orr-Sommerfeld operator is quite misleading. Because $\\mathcal{L}_{OS}$ is non-normal, its eigenfuntions are not orthogonal. As they decay at different rates, the total amplitude of some arbitrary initial condition can **increase**!\n",
+    "The answer is that the *spectrum* of the Orr-Sommerfeld operator is quite misleading. Because the linear operator $\\mathcal{L}_{OS}$ is non-normal, its eigenfuntions are not orthogonal. As they decay at different rates, the total amplitude of some arbitrary initial condition can **increase**!\n",
     "\n",
     "This observation was first noted in the pioneering paper by [Reddy, Schmid, and Henningson (1993)](https://www.jstor.org/stable/2102271?seq=1). For a detailed discussion of the Orr-Sommerfeld problem including a thorough discussion of transient growth, see [Schmid and Henningson (2001)](https://link.springer.com/book/10.1007/978-1-4613-0185-1)."
    ]
@@ -80,16 +80,16 @@
     "alpha = 1.\n",
     "Re = 10000\n",
     "\n",
-    "orr_somerfeld = de.EVP(d,['w','wz','wzz','wzzz'],'c')\n",
+    "orr_somerfeld = de.EVP(d,['w', 'wz', 'wzz', 'wzzz'], 'c')\n",
     "orr_somerfeld.parameters['alpha'] = alpha\n",
     "orr_somerfeld.parameters['Re'] = Re\n",
     "orr_somerfeld.substitutions['umean']= '1 - z**2'\n",
     "orr_somerfeld.substitutions['umeanzz']= '-2'\n",
     "\n",
     "orr_somerfeld.add_equation('dz(wzzz) - 2*alpha**2*wzz + alpha**4*w -1j*alpha*Re*((umean-c)*(wzz - alpha**2*w) - umeanzz*w) = 0')\n",
-    "orr_somerfeld.add_equation('dz(w)-wz = 0')\n",
-    "orr_somerfeld.add_equation('dz(wz)-wzz = 0')\n",
-    "orr_somerfeld.add_equation('dz(wzz)-wzzz = 0')\n",
+    "orr_somerfeld.add_equation('dz(w) - wz = 0')\n",
+    "orr_somerfeld.add_equation('dz(wz) - wzz = 0')\n",
+    "orr_somerfeld.add_equation('dz(wzz) - wzzz = 0')\n",
     "\n",
     "orr_somerfeld.add_bc('left(w) = 0')\n",
     "orr_somerfeld.add_bc('right(w) = 0')\n",
@@ -199,11 +199,11 @@
     }
    ],
    "source": [
-    "plt.scatter(EP.evalues.real, EP.evalues.imag,color='red')\n",
+    "plt.scatter(EP.evalues.real, EP.evalues.imag, color='red')\n",
     "plt.axis('square')\n",
-    "plt.axhline(0,color='k',alpha=0.2)\n",
-    "plt.xlim(0,1)\n",
-    "plt.ylim(-1,0.1)\n",
+    "plt.axhline(0, color='k', alpha=0.2)\n",
+    "plt.xlim(0, 1)\n",
+    "plt.ylim(-1, 0.1)\n",
     "plt.xlabel(r\"$c_r$\")\n",
     "plt.ylabel(r\"$c_i$\")"
    ]
@@ -250,8 +250,8 @@
     "k = 100\n",
     "\n",
     "psize = 100\n",
-    "real_points = np.linspace(0,1, psize)\n",
-    "imag_points = np.linspace(-1,0.1, psize)"
+    "real_points = np.linspace(0, 1, psize)\n",
+    "imag_points = np.linspace(-1, 0.1, psize)"
    ]
   },
   {
@@ -265,7 +265,7 @@
     "\n",
     "Here, we will use the energy inner product for the Orr-Sommerfeld form of the equations,\n",
     "\n",
-    "$$<\\mathbf{X}_1|\\mathbf{X}_2>_E = \\int_{-1}^{1} w_1^\\dagger w_2 + \\frac{(\\partial_z w_1)^\\dagger \\partial_z w_2}{\\alpha^2} dz,$$\n",
+    "$$\\langle\\mathbf{X}_1|\\mathbf{X}_2\\rangle_E = \\int_{-1}^{1} w_1^\\dagger w_2 + \\frac{(\\partial_z w_1)^\\dagger \\partial_z w_2}{\\alpha^2} dz,$$\n",
     "\n",
     "with $\\mathbf{X}_{1,2}$ state vectors containing $w$ and its first three derivatives. This inner product induces a norm $||\\mathbf{X}||_E$ that is equal to the energy of the perturbation.\n",
     "\n",
@@ -280,9 +280,9 @@
    "outputs": [],
    "source": [
     "def energy_norm(X1, X2):\n",
-    "    u1 = X1['wz']/alpha\n",
+    "    u1 = X1['wz'] / alpha\n",
     "    w1 = X1['w']\n",
-    "    u2 = X2['wz']/alpha\n",
+    "    u2 = X2['wz'] / alpha\n",
     "    w2 = X2['w']\n",
     "    \n",
     "    field = (np.conj(u1)*u2 + np.conj(w1)*w2).evaluate().integrate()\n",
@@ -347,13 +347,13 @@
     }
    ],
    "source": [
-    "plt.scatter(EP.evalues.real, EP.evalues.imag,color='red')\n",
-    "plt.contour(EP.ps_real,EP.ps_imag, np.log10(EP.pseudospectrum),levels=np.arange(-8,0),colors='k')\n",
+    "plt.scatter(EP.evalues.real, EP.evalues.imag, color='red')\n",
+    "plt.contour(EP.ps_real,EP.ps_imag, np.log10(EP.pseudospectrum), levels=np.arange(-8, 0), colors='k')\n",
     "plt.axis('square')\n",
-    "plt.axhline(0,color='k',alpha=0.2)\n",
+    "plt.axhline(0, color='k', alpha=0.2)\n",
     "\n",
-    "plt.xlim(0,1)\n",
-    "plt.ylim(-1,0.1)\n",
+    "plt.xlim(0, 1)\n",
+    "plt.ylim(-1, 0.1)\n",
     "plt.xlabel(r\"$c_r$\")\n",
     "plt.ylabel(r\"$c_i$\")"
    ]
@@ -432,7 +432,7 @@
     "\n",
     "However, the energy norm in primitive variables is slightly different:\n",
     "\n",
-    "$$<\\mathbf{X}_1|\\mathbf{X}_2>_E = \\int_{-1}^{1} u_1^\\dagger u_2 + w_1^\\dagger w_2 dz,$$\n",
+    "$$\\langle\\mathbf{X}_1|\\mathbf{X}_2\\rangle_E = \\int_{-1}^{1} u_1^\\dagger u_2 + w_1^\\dagger w_2 dz,$$\n",
     "\n",
     "so we write another `energy_norm` function and then call `prim_EP.calc_ps` with the same arguments as before."
    ]
@@ -492,20 +492,20 @@
     }
    ],
    "source": [
-    "clevels = np.arange(-8,0)\n",
-    "OS_CS = plt.contour(EP.ps_real,EP.ps_imag, np.log10(EP.pseudospectrum),levels=clevels,colors='blue',linestyles='solid', alpha=0.5)\n",
-    "P_CS = plt.contour(prim_EP.ps_real,prim_EP.ps_imag, np.log10(prim_EP.pseudospectrum),levels=clevels,colors='red',linestyles='solid', alpha=0.5)\n",
+    "clevels = np.arange(-8, 0)\n",
+    "OS_CS = plt.contour(EP.ps_real, EP.ps_imag, np.log10(EP.pseudospectrum), levels=clevels, colors='blue', linestyles='solid', alpha=0.5)\n",
+    "P_CS = plt.contour(prim_EP.ps_real, prim_EP.ps_imag, np.log10(prim_EP.pseudospectrum), levels=clevels, colors='red', linestyles='solid', alpha=0.5)\n",
     "\n",
-    "plt.scatter(prim_EP.evalues.real, prim_EP.evalues.imag,color='red',zorder=3)\n",
-    "plt.scatter(EP.evalues.real, EP.evalues.imag,color='blue',marker='x',zorder=4)\n",
+    "plt.scatter(prim_EP.evalues.real, prim_EP.evalues.imag, color='red', zorder=3)\n",
+    "plt.scatter(EP.evalues.real, EP.evalues.imag, color='blue', marker='x',zorder=4)\n",
     "\n",
-    "lines = [OS_CS.collections[0],P_CS.collections[0]]\n",
+    "lines = [OS_CS.collections[0], P_CS.collections[0]]\n",
     "labels = ['O-S form', 'primitive form']\n",
     "\n",
     "plt.axis('square')\n",
-    "plt.xlim(0,1)\n",
-    "plt.ylim(-1,0.1)\n",
-    "plt.axhline(0,color='k',alpha=0.2)\n",
+    "plt.xlim(0, 1)\n",
+    "plt.ylim(-1, 0.1)\n",
+    "plt.axhline(0, color='k', alpha=0.2)\n",
     "plt.xlabel(r\"$c_r$\")\n",
     "plt.ylabel(r\"$c_i$\")\n",
     "plt.legend(lines, labels)\n",

docs/pages/getting_started.rst

@@ -11,10 +11,10 @@ This is not quite as trivial a problem as it might seem, because we are expandin
     import dedalus.public as de
 
     Nx = 128
-    x = de.Chebyshev('x',Nx, interval=(-1, 1))
+    x = de.Chebyshev('x' ,Nx, interval=(-1, 1))
     d = de.Domain([x])
 
-    string = de.EVP(d, ['u','u_x'], eigenvalue='omega')
+    string = de.EVP(d, ['u', 'u_x'], eigenvalue='omega')
     string.add_equation("omega*u + dx(u_x) = 0")
     string.add_equation("u_x - dx(u) = 0")
     string.add_bc("left(u) = 0")
@@ -24,7 +24,7 @@ This is not quite as trivial a problem as it might seem, because we are expandin
     EP.solve(sparse=False)
     ax = EP.plot_spectrum()
     print("there are {} good eigenvalues.".format(len(EP.evalues)))
-    ax.set_ylim(-1,1)
+    ax.set_ylim(-1, 1)
     ax.figure.savefig('waves_spectrum.png')
 
     ax = EP.plot_drift_ratios()
@@ -38,15 +38,15 @@ That code takes about 10 seconds to run on a 2020 Core-i7 laptop, produces about
 
 eigentools has taken a Dedalus eigenvalue problem, automatically run it at 1.5 times the specified resolution, rejected any eigenvalues that do not agree to a default precision of one part in :math:`10^{-6}` and plotted a spectrum in six extra lines of code!           
 
-Most of the plotting functions in eigentools return a `matplotlib` `axes` object, making it easy to modify the plot defaults.
-Here, we set the y-limits manually, because the eigenvalues of a string are real.
-Try removing the `ax.set_ylim(-1,1)` line and see what happens.
+Most of the plotting functions in eigentools return a :code:`matplotlib` :code:`axes` object, making it easy to modify the plot defaults.
+Here, we set the :math:y`-limits manually, because the eigenvalues of a string are real.
+Try removing the :code:`ax.set_ylim(-1, 1)` line and see what happens.
 
 Mode Rejection
 --------------
-One of the most important tasks eigentools performs is spurious mode rejection. It does so by computing the "drift ratio" [Boyd2000]_ between the eigenvalues at the given resolution and a higher resolution problem that eigentools automatically assembles. By default, the "high" resolution case is 1.5 times the given resolution, though this is user configurable via the `factor` keyword option to `Eigenproblem()`.
+One of the most important tasks eigentools performs is spurious mode rejection. It does so by computing the "drift ratio" [Boyd2000]_ between the eigenvalues at the given resolution and a higher resolution problem that eigentools automatically assembles. By default, the "high" resolution case is 1.5 times the given resolution, though this is user configurable via the :code:`factor` keyword option to `Eigenproblem()`.
 
-The drift ratio :math:`\delta` is calculated using either the **ordinal** (e.g. first mode of low resolution to first mode of high resolution) or **nearest** (mode with smallest difference between a given high mode and all low modes). In order to visualize this, `EP.plot_drift_ratios()` in the above code returns an `axes` object making a plot of the *inverse drift ratio* (:math:`1/\delta`),
+The drift ratio :math:`\delta` is calculated using either the **ordinal** (e.g. first mode of low resolution to first mode of high resolution) or **nearest** (mode with smallest difference between a given high mode and all low modes). In order to visualize this, :code:`EP.plot_drift_ratios()` in the above code returns an :code:`axes` object making a plot of the *inverse drift ratio* (:math:`1/\delta`),
 
 .. image:: ../images/waves_drift_ratios.png
            :width: 400

docs/pages/upgrading.rst

@@ -31,7 +31,7 @@ Here, we solve twice, once with :code:`c1 = 1` and once with :code:`c2 = 2`. Giv
     x = de.Chebyshev('x',Nx, interval=(-1, 1))
     d = de.Domain([x])
 
-    string = de.EVP(d, ['u','u_x'], eigenvalue='omega')
+    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")
@@ -40,7 +40,7 @@ Here, we solve twice, once with :code:`c1 = 1` and once with :code:`c2 = 2`. Giv
     EP = Eigenproblem(string)
     EP.solve(sparse=False)
     evals_c1 = EP.evalues
-    EP.solve(sparse=False, parameters={'c2':2})
+    EP.solve(sparse=False, parameters={'c2': 2})
     evals_c2 = EP.evalues
 
     print(np.allclose(evals_c2, 2*evals_c1))
@@ -49,7 +49,7 @@ Getting eigenmodes
 ==================
 
 Getting eigenmodes has also been simplified and significantly extended. Previously, getting an eigenmode corresponding to an eigenvalue required using the :code:`set_state()` method on the underlying :code:`EVP` object. In keeping with the principle of not needing to manipulate the :code:`EVP`, we provide a new :code:`.eigenmode(index)`, where :code:`index` is the mode number corresponding to the eigenvalue index in :code:`EP.evalues`. By default, with mode rejection on, these are the "good" eigenmodes.
-`.eigenmode(index)` returns a Dedalus :code:`FieldSystem` object, with a Dedalus :code:`Field` for each field in the eigenmode:
+:code:`.eigenmode(index)` returns a Dedalus :code:`FieldSystem` object, with a Dedalus :code:`Field` for each field in the eigenmode:
 
 .. code-block:: python
 
@@ -80,8 +80,8 @@ First, replace
         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
+        gr, indx, freq = EP.growth_rate({"Q":y, "iRm": iRm, "iR": iRe})
+        ret = gr + 1j*freq
         return ret
 
     cf = CriticalFinder(shim, comm)
@@ -102,7 +102,7 @@ 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))
+    ns   = np.array((10, 10))
     logs = np.array((False, False))
 
     cf.grid_generator(mins, maxs, ns, logs=logs)

examples/mri.py

@@ -94,7 +94,7 @@
     pax,cax = cf.plot_crit()
     fig = pax.figure
     # add an interpolated critical line
-    x_lim = cf.parameter_grids[0][0,np.isfinite(cf.roots)]
+    x_lim = cf.parameter_grids[0][0, np.isfinite(cf.roots)]
     x_hires = np.linspace(x_lim[0], x_lim[-1], 100)
     pax.plot(x_hires, cf.root_fn(x_hires), color='k')
     fig.savefig('{}.png'.format(file_name), dpi=300)

examples/orr_sommerfeld.py

@@ -6,7 +6,6 @@
 sigma = -1j*alpha*Re*lambda,
 
 where sigma is our eigenvalue and Lambda is Orszag's.
-
 """
 import matplotlib
 matplotlib.use('Agg')
@@ -26,17 +25,17 @@
 
 # Define the Orr-Somerfeld problem in Dedalus: 
 
-z = de.Chebyshev('z',50)
+z = de.Chebyshev('z', 50)
 d = de.Domain([z],comm=MPI.COMM_SELF)
 
-orr_somerfeld = de.EVP(d,['w','wz','wzz','wzzz'],'sigma')
+orr_somerfeld = de.EVP(d,['w', 'wz', 'wzz', 'wzzz'], 'sigma')
 orr_somerfeld.parameters['alpha'] = 1.
 orr_somerfeld.parameters['Re'] = 10000.
 
 orr_somerfeld.add_equation('dz(wzzz) - 2*alpha**2*wzz + alpha**4*w - sigma*(wzz-alpha**2*w)-1j*alpha*(Re*(1-z**2)*(wzz-alpha**2*w) + 2*Re*w) = 0 ')
-orr_somerfeld.add_equation('dz(w)-wz = 0')
-orr_somerfeld.add_equation('dz(wz)-wzz = 0')
-orr_somerfeld.add_equation('dz(wzz)-wzzz = 0')
+orr_somerfeld.add_equation('dz(w) - wz = 0')
+orr_somerfeld.add_equation('dz(wz) - wzz = 0')
+orr_somerfeld.add_equation('dz(wzz) - wzzz = 0')
 
 orr_somerfeld.add_bc('left(w) = 0')
 orr_somerfeld.add_bc('right(w) = 0')
@@ -89,4 +88,3 @@
 
     cf.save_grid('orr_sommerfeld_growth_rates')
     cf.plot_crit()
-