From 4401edbbdff07df8ded78c721f509becf26719ab Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 17 Jun 2023 20:43:25 +1000 Subject: [PATCH 1/6] linear -> exponential plus add spaces after commas --- .../Orr Somerfeld pseudospectra.ipynb | 60 +++++++++---------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/docs/notebooks/Orr Somerfeld pseudospectra.ipynb b/docs/notebooks/Orr Somerfeld pseudospectra.ipynb index b8013f4..2a80a5e 100644 --- a/docs/notebooks/Orr Somerfeld pseudospectra.ipynb +++ b/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", From d16c78a02237c3dc88576b7ef1bc67aa434016c2 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 17 Jun 2023 20:50:29 +1000 Subject: [PATCH 2/6] spaces after commas --- docs/notebooks/Mixed Layer Instability.ipynb | 32 ++++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/docs/notebooks/Mixed Layer Instability.ipynb b/docs/notebooks/Mixed Layer Instability.ipynb index 627b848..cb36e4e 100644 --- a/docs/notebooks/Mixed Layer Instability.ipynb +++ b/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)" ] }, { From a94134d3028e5da227a6fa48c0134278c986e5d7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 17 Jun 2023 20:59:43 +1000 Subject: [PATCH 3/6] code formatting --- docs/pages/getting_started.rst | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/pages/getting_started.rst b/docs/pages/getting_started.rst index c165e04..082bc92 100644 --- a/docs/pages/getting_started.rst +++ b/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 From b91e1d6c1f8ae460720bcf34b0ce9aa911dcb921 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 17 Jun 2023 21:03:15 +1000 Subject: [PATCH 4/6] fix code rendering --- docs/pages/upgrading.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/pages/upgrading.rst b/docs/pages/upgrading.rst index 60acc15..46053ef 100644 --- a/docs/pages/upgrading.rst +++ b/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) From 2185586d218758b4b533f7cfa7ce2ab254d5ab01 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 18 Jun 2023 07:01:42 +1000 Subject: [PATCH 5/6] add space after comma --- examples/mri.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/mri.py b/examples/mri.py index 9706394..261e4aa 100644 --- a/examples/mri.py +++ b/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) From e8b3224c6feaf804137c7073f6f09e01470a160a Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 18 Jun 2023 07:03:09 +1000 Subject: [PATCH 6/6] add spaces --- examples/orr_sommerfeld.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/examples/orr_sommerfeld.py b/examples/orr_sommerfeld.py index c131b5f..a22dcff 100644 --- a/examples/orr_sommerfeld.py +++ b/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() -