[example] Internal wave in saturated conditions #192
Conversation
…th/Breeze.jl into glw/saturated-internal-wave
|
I'm trying to fix this example because the AI generation doesn't even run ... it's not that good tbh. What do we want to demonstrate? How a small-amplitude wavepacket moves and what dispersion relationship would follow? |
…th/Breeze.jl into glw/saturated-internal-wave
| #= | ||
| @info "Simulation complete" | ||
|
|
||
| # ## Diagnostics: observed vs predicted buoyancy frequency | ||
|
|
||
| w_series = FieldTimeSeries(filename, "w") | ||
| θ_series = FieldTimeSeries(filename, "θ") | ||
| q_series = FieldTimeSeries(filename, "qᵗ") | ||
|
|
||
| times = w_series.times | ||
| Nt = length(times) | ||
| i_mid = Nx ÷ 2 | ||
| k_mid = Nz ÷ 2 | ||
|
|
||
| @allowscalar begin | ||
| w_signal = [w_series[i_mid, 1, k_mid, n] for n in 1:Nt] | ||
| end | ||
|
|
||
| function dominant_frequency(times, signal) | ||
| demeaned = signal .- mean(signal) | ||
| dt = mean(diff(times)) | ||
| spectrum = abs.(rfft(demeaned)) | ||
| freqs = (0:length(spectrum)-1) ./ (dt * length(signal)) | ||
| idx = argmax(spectrum[2:end]) + 1 # skip zero frequency | ||
| return 2π * freqs[idx] | ||
| end | ||
|
|
||
| ω_observed = dominant_frequency(times, w_signal) | ||
| period_observed = 2π / ω_observed | ||
|
|
||
| @info "Wave frequencies" ω_dry ω_moist ω_observed | ||
|
|
||
| # ## Visualization | ||
|
|
||
| x = range(grid.x[1], grid.x[2]; length = Nx) | ||
| z = range(grid.z[1], grid.z[2]; length = Nz) | ||
| time_minutes = times ./ 60 | ||
|
|
||
| w_snapshot = Array(permutedims(w_series[:, 1, :, end])) | ||
| w_amplitude = maximum(abs, w_signal) | ||
|
|
||
| dry_fit = w_amplitude * sin.(ω_dry .* (times .- times[1])) | ||
| moist_fit = w_amplitude * sin.(ω_moist .* (times .- times[1])) | ||
|
|
||
| fig = Figure(size = (900, 950), fontsize = 14) | ||
|
|
||
| ax_field = Axis(fig[1, 1]; xlabel = "x (km)", ylabel = "z (km)", | ||
| title = "Vertical velocity w at t = $(prettytime(times[end]))") | ||
| hm = heatmap!(ax_field, x ./ 1_000, z ./ 1_000, w_snapshot; | ||
| colormap = :balance, colorrange = (-w_amplitude, w_amplitude)) | ||
| Colorbar(fig[1, 2], hm, label = "w (m s⁻¹)") | ||
|
|
||
| ax_ts = Axis(fig[2, 1]; xlabel = "time (min)", ylabel = "w (m s⁻¹)", | ||
| title = "Oscillation at domain center") | ||
| lines!(ax_ts, time_minutes, w_signal; label = "simulation", color = :black) | ||
| lines!(ax_ts, time_minutes, dry_fit; label = "dry N", color = :royalblue, linestyle = :dash) | ||
| lines!(ax_ts, time_minutes, moist_fit; label = "saturated N", color = :firebrick, linestyle = :dot) | ||
| axislegend(ax_ts; position = :rt) | ||
|
|
||
| ax_freq = Axis(fig[3, 1]; xlabel = "", ylabel = "ω (s⁻¹)", | ||
| title = "Dominant frequency comparison", | ||
| xticks = (1:3, ["observed", "dry", "saturated"])) | ||
| barplot!(ax_freq, 1:3, [ω_observed, ω_dry, ω_moist]; | ||
| color = (:gray50, :royalblue, :firebrick)) | ||
|
|
||
| fig | ||
| =# |
There was a problem hiding this comment.
This is messing up with the rendering of the documentation: https://numericalearth.github.io/BreezeDocumentation/previews/PR192/literated/moist_internal_wave/###-Diagnostics:-observed-vs-predicted-buoyancy-frequency
|
@glwagner I think we are committing over each other. I can't see the things I was doing yesterday... I'll leave it to you? |
|
I don't think Greg has pushed in a while: https://github.com/NumericalEarth/Breeze.jl/activity?ref=glw/saturated-internal-wave 🤔 |
|
There is a chance I committed over myself? Perhaps I need to clear my head... |
|
Maybe? There are no force-pushes at all on this branch, and you're the only one who has pushed over the last couple of days 😅 |
|
This might be the culprit ;) |
|
I think I was working on two files! |
This PR adds a new example that builds an internal wave in saturated conditions. It also adds some features that help compute the saturation specific humidity in various configurations:
the addition of a utility for the last case should be useful for idealized problems I think.