Merge branch 'main' of https://github.com/johnryantaylor/FDSE

Author: johnryantaylor

Project3/README.md

@@ -13,6 +13,7 @@ Here, we will examine the basic stability properties of a stratified shear flow,
 Start by considering a stratified shear flow of the form:
 
 $$\mathbf{U}=S_0h \text{tanh}\left(\frac{z-LZ/2}{h}\right)\hat{\imath},$$
+
 and
 
 $$B=N_0^2h \text{tanh}\left(\frac{z-LZ/2}{h}\right),$$
@@ -46,7 +47,9 @@ Now, let's simulate K-H instability in Oceananigans. The script `KH.jl` simulate
 
 See if you can verify the growth rate and the most unstable mode that you found from `linstab.jl` in the Oceananigans simulations. Note that the size of the computational domain in the x-direction quantizes the number of wavelengths associated with the perturbations. You should adjust the size of the domain in `KH.jl` to be an integer multiple of the wavelength of the most unstable mode that you found from `linstabl.jl` (be sure to match the other parameters if you changed them). Also, the perturbation that is applied to the velocity field (the parameter `kick` in `KH.jl`) is very large, and you should make it smaller to ensure that the linear approximation used in the linear stability analysis is valid.
 
-Once you are happy that `KH.jl` is capturing the most unstable mode, you can verify that the perturbation growth rate is comparable to the value that you got in `linstab.jl`. An easy way to do this is to plot the logarithm of the vertical velocity variance, $w^2$, as a function of time. Since the basic state has $w=0$, the vertical velocity does not have a contribution from the basic state and the amplitude of $w^2$ should grow like $e^{2\sigma t}$. You will need to add a line of code to calculate and then save the vertical velocity variance. You can do this either in `KH.jl` and save the vertical velocity variance to a file, or do the calculation in `plot_KH.jl` (the latter approach is probably easier). You might also like to plot the evolution of the perturbation kinetic energy (including both $u$ and $w$), where the perturbation is the departure from an x-average. Since $u$ has a non-zero $x$-average, we will need to remove the x-average of $u$ to get the perturbation, $u'$. To calculate an x-average of an array (say `u_xz`), install and use the Statistics package (`] add Statistics` followed by `using Statistics`), and then use `mean(u_xz, dims = 1)` which averages the array `u_xz` over its first dimension (in this case x). Plot a timeseries of the vertical velocity variance or the perturbation kinetic energy and overlay a line or curve indicating the growth rate from `linstab.jl`.
+Once you are happy that `KH.jl` is capturing the most unstable mode, you can verify that the perturbation growth rate is comparable to the value that you got in `linstab.jl`. An easy way to do this is to plot the logarithm of the vertical velocity variance, $w^2$, as a function of time. Since the basic state has $w=0$, the vertical velocity does not have a contribution from the basic state and the amplitude of $w^2$ should grow like $e^{2\sigma t}$. You will need to add a line of code to calculate and then save the vertical velocity variance. You can do this either in `KH.jl` and save the vertical velocity variance to a file, or do the calculation in `plot_KH.jl` (the latter approach is probably easier). 
+
+Plot a timeseries of the vertical velocity variance and overlay a line or curve indicating the growth rate from `linstab.jl`. You will probably find that the growth rate measured in `KH.jl` is smaller than the growth rate predicted in `linstab.jl`. This is because the shear layer diffuses in `KH.jl`. You can get a better match by increasing the Reynolds number and changing the advection scheme from `UpwindBiased()` to `Centered(order=4)`. The resulting 4th order centered finite difference scheme has less numerical diffusion than the upwind scheme that we have been using. Note, however, that the simulation will be under-resolved when the flow becomes nonlinear.
 
 Finally, try making a plot of the gradient Richardson number as a function of depth and time where the shear and buoyancy frequency are calculated from x-averages of the `u` and `b` fields. Note that Oceananigans includes gradient operators, e.g. $\partial z (u)$ returns the derivative of `u` in the vertical direction. How does the value of the gradient Richardson number compare with what you might have expected based on the Miles-Howard theorem.