<jupyter_start><jupyter_text>OCEAN SCIENCES 272:  Problem Set 5<jupyter_code># imports
import numpy as np
from importlib import reload

from scipy import optimize

import unyt

from oceanpy import constants
from oceanpy.forces import coriolis

from os_classes.ocea257 import chapter1, chapter2
from os_classes import plotting as os_plotting

from bokeh import plotting
from bokeh import models
from bokeh import layouts

from IPython.display import Image
from IPython.core.display import HTML<jupyter_output><empty_output><jupyter_text>PS5-1-- 3 Waves  We are to consider a set of 3 waves: $\psi_j(x,y,y) = \phi_j(y) \sin[k_j x - \omega t]$ each growing exponentially with  $\phi_j(y,t) = A_j \exp(k_j c_{ij} t)$ and $A_j = 0.001$ at $t=0$.  The wavenumbers are given by: $\lambda_j = 2\pi/k_j = j * 1000$m and the imaginary part is given by: $c_{ij} = 100, 300, 350$m/day<jupyter_code># Setup
k_1 = 2 * np.pi / (1000*unyt.m)
k_2 = 2 * np.pi / (2000*unyt.m)
k_3 = 2 * np.pi / (3000*unyt.m)
ks = [k_1, k_2, k_3]

c_1 = 100 * unyt.m/unyt.day
c_2 = 300 * unyt.m/unyt.day
c_3 = 350 * unyt.m/unyt.day
cs = [c_1, c_2, c_3]<jupyter_output><empty_output><jupyter_text>We are asked which wave reaches amplitude 1 first.  This will be the one that maximizes $k_j c_{ij}$<jupyter_code>for j, c,k in zip(range(3), cs, ks):
    print("j={},  ck={}".format(j+1, c*k))<jupyter_output>j=1,  ck=0.6283185307179586 1/day
j=2,  ck=0.942477796076938 1/day
j=3,  ck=0.7330382858376183 1/day<jupyter_text>It is clear that the $j=2$ wave grows the fastest. It will reach $A=1$ at $t = \ln(10^3) / (k_jc_{ij})$<jupyter_code>t_A2_1 = np.log(1e3) / (k_2 * c_2)
t_A2_1<jupyter_output><empty_output><jupyter_text>At $t = 7.33$ days, the other two waves have amplitudes: $A_j = A(t=0) \exp(k_j c_{ij} t)$<jupyter_code>for j, c,k in zip(range(3), cs, ks):
    A_j = 1e-3 * np.exp(c*k*t_A2_1)
    print("j={},  A_j = {}".format(j+1, A_j))<jupyter_output>j=1,  A_j = 0.09999999999999996
j=2,  A_j = 0.9999999999999998
j=3,  A_j = 0.2154434690031881<jupyter_text>So $A_1 = 0.1$ and $A_3 = 0.215$  $A_2 = 10$ occurs $\ln(10) / (k_2 c_{i2})$ days later<jupyter_code>t_A2_10 = np.log(1e4) / (k_2 * c_2)
t_A2_10
for j, c,k in zip(range(3), cs, ks):
    A_j = 1e-3 * np.exp(c*k*t_A2_10)
    print("j={},  A_j = {}".format(j+1, A_j))<jupyter_output>j=1,  A_j = 0.46415888336127775
j=2,  A_j = 10.000000000000009
j=3,  A_j = 1.2915496650148832<jupyter_text>Because the waves are all in phase for $x < 1000$m, the superposition is largest in that region. Yet, the exponential growth of the $j=2$ wave dominates and is effectively the solution at late times. Here is an interactive graph:<jupyter_code>x = np.linspace(0, 3000., 3000) * unyt.m 

psi_1 = 1e-3 * np.sin(k_1 * x)
psi_2 = 1e-3 * np.sin(k_2 * x)
psi_3 = 1e-3 * np.sin(k_3 * x)
source = plotting.ColumnDataSource(data=dict(x=x, y=psi_1+psi_2+psi_3))
# Plot
plotting.output_notebook()

p = plotting.figure(title='Three waves', x_axis_label='x',
                    #x_axis_type="log", #y_axis_type="log",
                   y_axis_label='A', plot_width=600)#, x_range=(10**-4, 1.),
                   #y_range=(1e-12, 1)))
    
p.line('x', 'y', source=source, line_width=2, color='blue')

time_slider = models.Slider(start=0, end=20, value=0, step=1, title="time (days)")

callback = models.CustomJS(args=dict(source=source, time=time_slider),
                    code="""
    const data = source.data;
    const t = time.value;
    const x = data['x']
    const y = data['y']
    const ck1 = 0.62831853
    const ck2 = 0.9424777
    const ck3 = 0.7330382
    const k1 = 0.00628319
    const k2 = 0.00314159
    const k3 = 0.0020944
    const A = 0.001
    for (var i = 0; i < x.length; i++) {
        y[i] = A*Math.sin(k1*x[i])*Math.exp(ck1*t) + A*Math.sin(k2*x[i])*Math.exp(ck2*t) + A*Math.sin(k3*x[i])*Math.exp(ck3*t);
    }
    source.change.emit();
""")

time_slider.js_on_change('value', callback)

layout = layouts.row(p, time_slider)

plotting.show(layout)<jupyter_output><empty_output><jupyter_text>----  PS5-2 (C-R 10-2)   Recall that a *necessary* condition for instability is that $\beta_0 - d^2 \bar u / dy^2$ vanish somewhere within the domain. And, likewise, a *sufficient* condition for stability is that it does not.  Another *necessary* condition for stability is that $(\bar u - \bar u_0) (\beta_0 - d^2 \bar u / dy^2)$ be positive in the domain, where $\bar u_0$ is the value of $\bar u(y)$ where the expression vanishes.  Let us now go through the various functions:  $\bar u(y) = U \left ( 1 - \frac{y^2}{L^2} \right ) \;\; [-L \le y \le +L]$ $\frac{d^2 \bar u(y)}{d y^2} = \frac{-2U}{L}$ As $\beta_0$ is positive definite, $\beta_0 - d^2 \bar u / dy^2$ will also be positive definite.  Therefore, this flow field is *stable*.  $\bar u(y) = U \sin \frac{\pi y}{L} \;\; [0 \le y \le L]$ $\frac{d^2 \bar u(y)}{d y^2} = \frac{-U \pi^2}{L^2} \sin \frac{\pi y}{L}$ As $\sin \frac{\pi y}{L}$ is positive definite for this domain, then $\beta_0 - d^2 \bar u / dy^2$ will also be positive definite.  Therefore, this flow field is *stable*.  $\bar u(y) = U \cos \frac{\pi y}{L} \;\; [0 \le y \le L]$ $\frac{d^2 \bar u(y)}{d y^2} = \frac{-U \pi^2}{L^2} \cos \frac{\pi y}{L}$  $\cos \frac{\pi y}{L}$ flips sign in this domain, so $\beta_0 - d^2 \bar u / dy^2$ can also flip sign and vanish.   This would occur when $\beta_0 = \frac{U \pi^2}{L^2} \cos \frac{\pi y}{L}$ implying $\bar u_0 = \frac{\beta_0 L^2}{\pi^2}$ This sets a condition on $U$ for instability.  $\bar u(y) = U \tanh \frac{y}{L} \;\; [-\infty \le y \le \infty]$ $\frac{d^2 \bar u(y)}{d y^2} = \frac{-2U}{L^2} \frac{\sinh \frac{y}{L}}{\cosh^3 \frac{y}{L}}$ The 2nd derivative flips sign at $y=0$.  Therefore,  $\beta_0 - d^2 \bar u / dy^2$ can also flip sign and vanish.   The largest *negative* value of $d^2 \bar u/dy^2$ is $\approx 1.4 U / L^2$ so this sets an upper limit on $\beta_0$ for instability, i.e.  $\beta_0 < 1.4 U/L^2$ ----  PS5-3 (C-R 10-4) -- Jet Stream  We will parameterize the jet stream with: $\bar u(y) = U \, \exp \left ( -\frac{y^2}{2L^2} \right )$ where $y=0$ at $45^\circ$N and $\beta_0 = 1.61 \times 10^{-11} \, \rm m^{-1} s^{-1}$, $U = 40$m/s, and $L = 570$km.  This implies $\frac{d^2 \bar u(y)}{d y^2} = \frac{U}{L^2} \, \exp \left ( -\frac{y^2}{2L^2} \right ) \left [ \frac{y^2}{L^2} -1 \right ]$  Let's plot $\beta_0 - d^2 \bar u / dy^2$ in the domain $0 \le y \le 3000$km<jupyter_code># constants
U = 40 * unyt.m / unyt.s
L = 570 * unyt.km
beta_0 = 1.61e-11 / unyt.m / unyt.s
yval = np.linspace(0, 3000., 3000) * unyt.km
d2u_dy2 = (U/L**2) * np.exp(- yval**2 / 2 / L**2) * (yval**2 / L**2 - 1)
d2u_dy2.to('1/m/s')
# Plot
plotting.output_notebook()

p = plotting.figure(title='Jet Stream', x_axis_label='y',
                   # x_axis_type="log", #y_axis_type="log",
                    plot_width=600,
                   y_axis_label='beta - d^2 u/dy^2')#, x_range=(10**-4, 1.),
                   #y_range=(1e-12, 1))

# Plot em
#p.line(y, beta_0-d2u_dy2, legend='difference', line_width=2)
p.line(yval, beta_0.value - d2u_dy2.to('1/m/s').value, legend='difference', line_width=2)



p.legend.location = "top_right"

os_plotting.bokeh_set_fontsize(p, 16)

plotting.show(p)<jupyter_output><empty_output><jupyter_text>We can see that $\beta_0 - d^2 \bar u / dy^2$ crosses 0 at $y \approx$ 630 and 1600km.  Let's solve for this numerically<jupyter_code>U0 = float(U.value)
L0 = float(L.to('m').value)
beta_00 = float(beta_0.value)
def diff_bd2(y):
    d2u_dy2 = (U0/L0**2) * np.exp(- y**2 / 2 / L0**2) * (y**2 / L0**2 - 1)
    return np.abs(beta_00 - d2u_dy2)*1e11  # To deal with small values
#
#print(U0, L0, beta_00)<jupyter_output><empty_output><jupyter_text>First 630km<jupyter_code>optimize.minimize_scalar(diff_bd2, bounds=[600.*1e3, 650.*1e3], method='bounded', options=dict(xatol=1e-10))
y0 = 635626.69 / 1e3 * unyt.km 
y0<jupyter_output><empty_output><jupyter_text>Now near 1600km<jupyter_code>optimize.minimize_scalar(diff_bd2, bounds=[1580.*1e3, 1650.*1e3], method='bounded', options=dict(xatol=1e-10))
y1 =  1606600.8401 / 1e3 * unyt.km 
y1<jupyter_output><empty_output><jupyter_text>Let's evaluate $\bar u$ at these positions:<jupyter_code>u0 = U * np.exp(-y0**2/2/L**2)
u0
u1 = U * np.exp(-y1**2/2/L**2)
u1<jupyter_output><empty_output><jupyter_text>Seems we should take the latter (u1).  Let's see if we are positive<jupyter_code>u = U * np.exp(-yval**2/2/L**2)
# Plot
plotting.output_notebook()

p = plotting.figure(title='Jet Stream', x_axis_label='y',
                   # x_axis_type="log", #y_axis_type="log",
                    plot_width=600,
                   y_axis_label='(u-u0)(beta - d^2 u/dy^2)')#, x_range=(10**-4, 1.),
                   #y_range=(1e-12, 1))

# Plot em
#p.line(y, beta_0-d2u_dy2, legend='difference', line_width=2)
p.line(yval, (u-u1)*(beta_0.value - d2u_dy2.to('1/m/s').value), legend='difference', line_width=2)



p.legend.location = "top_right"

os_plotting.bokeh_set_fontsize(p, 16)

plotting.show(p)<jupyter_output><empty_output><jupyter_text>Looks like we are positive until the first crossing.   Instability reigns! ----  PS5-4 -- More on the jet stream!  (a) How can modifications to $U,L$ or $\beta_0$ affect the instability for the first criterion?  First $\beta$.  Our expression $\beta_0 - d^2 \bar u / dy^2$ bottoms out at $\approx - 4 \times 10^{-11}$/m/s Therefore, increasing $\beta_0$ by $\approx 2.5$ would lead to stability.  Alternatively, $d^2 \bar u / dy^2 \propto U/L^2$.  Reducing this ratio so that the quantity never reaches 0 would also render the jet stream to be stable.  (b) Consider marginally (i) stable and (ii) unstable jets and then run them westward.  i.e. what is the effect of swapping the sign of $U$?  Recall: $\frac{d^2 \bar u(y)}{d y^2} = \frac{U}{L^2} \, \exp \left ( -\frac{y^2}{2L^2} \right ) \left [ \frac{y^2}{L^2} -1 \right ]$  so, for $U<0$, we have: $\frac{d^2 \bar u(y)}{d y^2} = \frac{|U|}{L^2} \, \exp \left ( -\frac{y^2}{2L^2} \right ) \left [ 1 - \frac{y^2}{L^2} \right ]$  Let's plot $d^2 \bar u / dy^2$ of the East reduced to have a peak just below $\beta_0$<jupyter_code>d2u_dy2 = (U/L**2/5) * np.exp(- yval**2 / 2 / L**2) * (yval**2 / L**2 - 1)
# Plot
plotting.output_notebook()

p = plotting.figure(title='Jet Stream West', x_axis_label='y',
                   # x_axis_type="log", #y_axis_type="log",
                    plot_width=600,
                   y_axis_label='d^2 u/dy^2')#, x_range=(10**-4, 1.),
                   #y_range=(1e-12, 1))

# Plot em
#p.line(y, beta_0-d2u_dy2, legend='difference', line_width=2)
p.line(yval, d2u_dy2.to('1/m/s').value, legend='difference', line_width=2)



p.legend.location = "top_right"

os_plotting.bokeh_set_fontsize(p, 16)

plotting.show(p)<jupyter_output><empty_output><jupyter_text>Now send it westward<jupyter_code>d2u_dy2_west = (U/L**2/5) * np.exp(- yval**2 / 2 / L**2) * (1 - yval**2 / L**2)
# Plot
plotting.output_notebook()

p = plotting.figure(title='Jet Stream West', x_axis_label='y',
                   # x_axis_type="log", #y_axis_type="log",
                    plot_width=600,
                   y_axis_label='d^2 u/dy^2')#, x_range=(10**-4, 1.),
                   #y_range=(1e-12, 1))

# Plot em
#p.line(y, beta_0-d2u_dy2, legend='difference', line_width=2)
p.line(yval, d2u_dy2_west.to('1/m/s').value, legend='difference', line_width=2)



p.legend.location = "top_right"

os_plotting.bokeh_set_fontsize(p, 16)

plotting.show(p)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18219895287958116, 'fraction_numerical': 0.07190226876090751, 'mean_word_length': 3.517540402049665, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '268114', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5321, 'n_tokens_neox': 4842, 'n_words': 1187, 'path': 'OCEA-272/hw/OCEA272_PS5.ipynb', 'repo_name': 'profxj/os_classes', 'size': '11460'}