<jupyter_start><jupyter_text>Time averaged velocities<jupyter_code>for var in velocities:
    ds[var + "_mean"] = ds.sel(t=slice(15.0, None))[var].mean('t')
ds.sel(t=slice(11.0, None))<jupyter_output><empty_output><jupyter_text>Fluctuation velocities<jupyter_code>for var in velocities:
    ds["u" + var[0]] = ds[var + "_mean"] - ds[var]
ds<jupyter_output><empty_output><jupyter_text>Friction velocityComputed as a function of time![](http://glossary.ametsoc.org/w/images/7/75/Ams2001glos-Fe18.gif) or using the current variable names $u_* = [\langle{u'v'}\rangle^2 + \langle{w'v'}\rangle^2]^{1/4}$<jupyter_code>ds_w = ds.sel(y=ds.y.min())
avg_coords = ('x', 'z')

uv = (ds_w.ux * ds_w.uy).mean(avg_coords)
wv = (ds_w.uz * ds_w.uy).mean(avg_coords)

urms = urms1 = (ds.ux**2 + ds.uy**2 + ds.uz**2).mean(('x', 'y', 'z'))**(1/2)

ustar = ustar1 = (uv**2 + wv**2)**(1/4)
plt.figure(figsize=(5, 4))
(urms).plot()
plt.ylabel("$u_{rms}$")
plt.xlabel("$t$")
plt.title("Spatially averaged RMS velocity")<jupyter_output><empty_output><jupyter_text>Momentum fluxAt the end of the simulation$\langle u'v' \rangle / u_{rms}$<jupyter_code>ds_end = ds.sel(t=ds.t.max())
ds_end
avg_coords = ('x', 'z')

# mom_flux = (ds_end.ux * ds_end.uy + ds_end.uz * ds_end.uy).mean(avg_coords)
mom_flux = (ds_end.ux * ds_end.uy).mean(avg_coords)
ux2 = (ds_end.ux ** 2).mean(avg_coords)
urms = urms2 = (ds_end.ux**2 + ds_end.uy**2 + ds_end.uz**2).mean(avg_coords)**(1/2)

ubulk = ds_end.x_velocity.mean()

ds_w = ds_end.sel(y=ds.y.min())
uv = (ds_w.ux * ds_w.uy).mean(avg_coords)
wv = (ds_w.uz * ds_w.uy).mean(avg_coords)

ustar = ustar2 = (uv**2 + wv**2)**(1/4)

plt.figure(figsize=(3, 4))
# (- mom_flux / ustar**2).plot.line(y="y", yscale="log", label=r"$- \langle u'v'\rangle / u_{*}^2$")
# (ux2 / ustar**2).plot.line(y="y", yscale="log", label=r"$- \langle u'^2\rangle / u_{*}^2$")
(- mom_flux / ubulk**2).plot(y="y")

# plt.xlabel(r"$- \langle u'v' + w'v' \rangle / u_{*}^2$")
# plt.xlabel(r"$- \langle u'w'\rangle / u_{*}^2$")
plt.xlabel(r"$- \langle u'w'\rangle / |u|^2$")

plt.ylabel("$z/H$")
plt.title("Horizontally averaged momentum flux at $t=21.6$")
# plt.legend()<jupyter_output><empty_output><jupyter_text>Combined plot<jupyter_code>fig, axes = plt.subplots(1, 2, figsize=(10, 4))
ax1, ax2 = axes.ravel()

(urms1).plot(ax=ax1)
ax1.set_ylabel("$u_{rms}$")
ax1.set_xlabel("$t$")
ax1.set_title("Spatially averaged RMS velocity")

(- mom_flux / ubulk**2).plot.line(ax=ax2, y="y", label=r"$- \langle u'v'\rangle / u_{*}^2$")

# plt.xlabel(r"$- \langle u'v' + w'v' \rangle / u_{*}^2$")
plt.xlabel(r"$- \langle u'w'\rangle / \langle u \rangle^2$")
plt.ylabel("$z/H$")
plt.xlim(0, 0.002)
plt.title("Normalised mean momentum flux at $t=21.6$")
# plt.legend()

# plt.savefig("stats.eps")
plt.savefig("output/urms_t-uw_t.png")<jupyter_output><empty_output><jupyter_text>Averaged velocity profile<jupyter_code>avg_coords = ('x', 'z')

ux_mean = ds_end.x_velocity.mean(avg_coords)

(ux_mean / ustar).plot(y="y", xscale="log", yscale="log")
plt.xlabel(r"$\langle u / u_* \rangle$")
plt.ylabel("$z/H$")
nu = 1e-5
yplus = ds.y * ustar / nu
plt.loglog(yplus, (ux_mean / ustar))
plt.ylabel(r"$u^+$")
plt.xlabel("$y^+$")

kappa = 0.41
c = 60
plt.loglog(yplus, np.log(yplus) / kappa + c, 'k--')<jupyter_output><empty_output><jupyter_text>Averaged wall-normal velocity profile<jupyter_code>avg_coords = ('x', 'z', 't')

uy_mean = ds.y_velocity.mean(avg_coords)
line, = (uy_mean).plot(y="y", xscale="linear", yscale="linear")
line.set_marker('x')
plt.xlabel(r"$\langle w \rangle$")
plt.ylabel("$z/H$")
plt.savefig("output/w_z.png")
avg_coords = ('x', 'z', 't')

uy_mean = dsw.y_velocity.mean(avg_coords)
line, = (uy_mean).plot(y="y", xscale="linear", yscale="linear")
line.set_marker('x')
plt.xlabel(r"$\langle w \rangle$")
plt.ylabel("$z/H$")
plt.savefig("output/w_z.png")<jupyter_output><empty_output><jupyter_text>Averaged wall normal velocity **before** reducing pressure residual tolerance<jupyter_code>avg_coords = ('x', 'z')

uy_mean = ds.y_velocity.mean(avg_coords)
line, = (uy_mean).isel(t=23).plot(y="y", xscale="linear", yscale="linear")
line.set_marker('x')
plt.xlabel(r"$\langle w \rangle$")
plt.ylabel("$z/H$")
plt.title(r"pressure residual tolerance $10^{-5}$ , $t \approx 19$")
plt.savefig("output/w_z.png")<jupyter_output><empty_output><jupyter_text>Averaged-wall normal velocity **after** reducing pressure residual tolerance<jupyter_code>import xarray as xr
ds_div_check = xr.open_dataset("velocity_slices_for_div_check.nc", engine="h5netcdf")
avg_coords = ('x', 'z')

uy_mean = ds_div_check.isel(t=-1).y_velocity.mean(avg_coords)
line, = (uy_mean).plot(y="y", xscale="linear", yscale="linear")
line.set_marker('x')
plt.xlabel(r"$\langle w \rangle$")
plt.ylabel("$z/H$")
plt.title(r"pressure residual tolerance $10^{-10}$ , $t \approx 19$")
plt.savefig("output/w_z.png")<jupyter_output><empty_output><jupyter_text>Wall-shear stress from MOST<jupyter_code>avg_coords = ('x', 'z')

ux = dsw.x_velocity
uz = dsw.z_velocity

kappa = 0.41
z0, z1 = dsw.y[:2]
dz2 = (z0 + z1) / 2  # or dzw.y[:2].mean()

L = 1e-4 * 1  # roughness height

ux_dz2 = (ux.sel(y=z0) + ux.sel(y=z1)) / 2
uz_dz2 = (uz.sel(y=z0) + uz.sel(y=z1)) / 2

uh_dz2 = np.sqrt(ux_dz2 ** 2 + uz_dz2 ** 2)

tau_x = -kappa**2 * (ux_dz2 * uh_dz2) / np.log(dz2/L)**2
tau_z = -kappa**2 * (uz_dz2 * uh_dz2) / np.log(dz2/L)**2

tau = np.sqrt(tau_x**2 + tau_z**2)
tau.mean(avg_coords).plot()
plt.title(rf"$\tau$ (MOST) z = {float(dz2)}")
plt.savefig("output/tau_most.png")<jupyter_output><empty_output><jupyter_text>Wall-shear stress turbulence definition<jupyter_code>z = ds.y.min()
ds_w = ds.sel(y=z)
avg_coords = ('x', 'z')

uv = (ds_w.ux * ds_w.uy).mean(avg_coords)
wv = (ds_w.uz * ds_w.uy).mean(avg_coords)

urms = (ds.ux**2 + ds.uy**2 + ds.uz**2).mean(('x', 'y', 'z'))**(1/2)

tau_turb = (uv**2 + wv**2)**(1/2)

tau_turb.tail(t=9).plot()
plt.title(rf"$\tau$ (Turb); z = {float(z)}")
plt.savefig("output/tau_turb_z1.png")
z = dsw.y[2]
ds_w2 = dsw.sel(y=z)
avg_coords = ('x', 'z')

ux = ds_w2.x_velocity.mean("t") - ds_w2.x_velocity
uy = ds_w2.y_velocity.mean("t") - ds_w2.y_velocity
uz = ds_w2.z_velocity.mean("t") - ds_w2.z_velocity


uv = (ux * uy).mean(avg_coords)
wv = (uz * uy).mean(avg_coords)

tau_turb = (uv**2 + wv**2)**(1/2)
tau_turb.plot()
plt.title(rf"$\tau$ (Turb); z = {float(z)}")
plt.savefig("output/tau_turb_z2.png")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2323800593842788, 'fraction_numerical': 0.026879199874980466, 'mean_word_length': 5.5775950668037, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '112420', 'license': 'Apache-2.0', 'n_tokens_mistral': 3112, 'n_tokens_neox': 2970, 'n_words': 449, 'path': 'paraview/velocity_to_stats.ipynb', 'repo_name': 'ashwinvis/sandbox', 'size': '6399'}