-
Notifications
You must be signed in to change notification settings - Fork 207
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Dry rising thermal bubble verification experiment (#31)
Dry rising thermal bubble verification experiment
- Loading branch information
Showing
8 changed files
with
261 additions
and
29 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -25,4 +25,6 @@ notes/*.log | |
| notes/*.gz | ||
|
|
||
| # Various output | ||
| *.png | ||
| *.gif | ||
| *.mp4 | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,99 @@ | ||
| using JULES, Oceananigans | ||
| using Plots | ||
| using Printf | ||
|
|
||
| Nx = Ny = 1 | ||
| Nz = 32 | ||
| L = 10e3 | ||
|
|
||
| grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(2, 2, 2), | ||
| x=(0, L), y=(0, L), z=(0, L)) | ||
| ##### | ||
| ##### Dry adiabatic atmosphere | ||
| ##### | ||
|
|
||
| gas = IdealGas() | ||
| Rᵈ, cₚ, cᵥ = gas.Rᵈ, gas.cₚ, gas.cᵥ | ||
|
|
||
| g = 9.80665 | ||
| pₛ = 100000 | ||
| θₛ = 300 | ||
| πₛ = 1 | ||
|
|
||
| θ(x, y, z) = θₛ | ||
| π(x, y, z) = πₛ - g*z / (cₚ*θₛ) | ||
| p(x, y, z) = pₛ * π(x, y, z)^(cₚ/Rᵈ) | ||
| ρ(x, y, z) = pₛ / (Rᵈ * θ(x, y, z)) * π(x, y, z)^(cᵥ/Rᵈ) | ||
| Θ(x, y, z) = ρ(x, y, z) * θ(x, y, z) | ||
|
|
||
| ##### | ||
| ##### Create model | ||
| ##### | ||
|
|
||
| model = CompressibleModel(grid=grid, buoyancy=gas, reference_pressure=pₛ, | ||
| prognostic_temperature=ModifiedPotentialTemperature(), | ||
| tracers=(:Θᵐ,)) | ||
|
|
||
| ##### | ||
| ##### Set initial conditions | ||
| ##### | ||
|
|
||
| set!(model.density, ρ) | ||
| set!(model.tracers.Θᵐ, θ) | ||
|
|
||
| # Initial profiles | ||
| Θ₀_prof = model.tracers.Θᵐ[1, 1, 1:Nz] | ||
| ρ₀_prof = model.density[1, 1, 1:Nz] | ||
|
|
||
| # Arrays we will use to store a time series of ρw(z = L/2). | ||
| times = [model.clock.time] | ||
| ρw_ts = [model.momenta.ρw[1, 1, Int(Nz/2)]] | ||
|
|
||
| ##### | ||
| ##### Time step and keep plotting vertical profiles of ρθ′, ρw, and ρ′. | ||
| ##### | ||
|
|
||
| # @animate for i=1:100 | ||
| while model.clock.time < 500 | ||
| time_step!(model; Δt=0.5, Nt=10) | ||
|
|
||
| @show model.clock.time | ||
| Θ_prof = model.tracers.Θᵐ[1, 1, 1:Nz] .- Θ₀_prof | ||
| W_prof = model.momenta.ρw[1, 1, 1:Nz+1] | ||
| ρ_prof = model.density[1, 1, 1:Nz] .- ρ₀_prof | ||
|
|
||
| Θ_plot = plot(Θ_prof, grid.zC, xlim=(-2.5, 2.5), xlabel="Theta_prime", ylabel="z (m)", label="") | ||
| W_plot = plot(W_prof, grid.zF, xlim=(-1.2, 1.2), xlabel="rho*w", label="") | ||
| ρ_plot = plot(ρ_prof, grid.zC, xlim=(-0.01, 0.01), xlabel="rho", label="") | ||
|
|
||
| t_str = @sprintf("t = %d s", model.clock.time) | ||
| display(plot(Θ_plot, W_plot, ρ_plot, title=["" "$t_str" ""], layout=(1, 3), show=true)) | ||
|
|
||
| push!(times, model.clock.time) | ||
| push!(ρw_ts, model.momenta.ρw[1, 1, Int(Nz/2)]) | ||
| end | ||
|
|
||
| ##### | ||
| ##### Plot the initial profile and the hydrostatically balanced profile. | ||
| ##### | ||
|
|
||
| ρ∞_prof = model.density[1, 1, 1:Nz] | ||
|
|
||
| θ₀_prof = Θ₀_prof ./ ρ₀_prof | ||
| θ∞_prof = model.tracers.Θᵐ[1, 1, 1:Nz] ./ ρ∞_prof | ||
|
|
||
| θ_plot = plot(θ₀_prof, grid.zC, xlabel="theta (K)", ylabel="z (m)", label="initial", legend=:topleft) | ||
| plot!(θ_plot, θ∞_prof, grid.zC, label="balanced") | ||
|
|
||
| ρ_plot = plot(ρ₀_prof, grid.zC, xlabel="rho (kg/m³)", ylabel="z (m)", label="initial") | ||
| plot!(ρ_plot, ρ∞_prof, grid.zC, label="balanced") | ||
|
|
||
| display(plot(θ_plot, ρ_plot, show=true)) | ||
|
|
||
|
|
||
| ##### | ||
| ##### Plot timeseries of maximum ρw | ||
| ##### | ||
|
|
||
| plot(times, ρw_ts, linewidth=2, xlabel="time (s)", ylabel="rho*w(z=$(Int(L/2)) m)", label="", show=true) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
138 changes: 138 additions & 0 deletions
138
verification/dry_rising_thermal_bubble/dry_rising_thermal_bubble.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,138 @@ | ||
| """ | ||
| This example sets up a dry, warm thermal bubble perturbation in a uniform | ||
| lateral mean flow which buoyantly rises. | ||
| """ | ||
|
|
||
| using Printf | ||
| using Plots | ||
| using Oceananigans | ||
| using JULES | ||
| using JULES: Π | ||
|
|
||
| const km = 1000 | ||
| const hPa = 100 | ||
|
|
||
| Lx = 20km | ||
| Lz = 10km | ||
|
|
||
| Δ = 125 # grid spacing [m] | ||
|
|
||
| Nx = Int(Lx/Δ) | ||
| Ny = 1 | ||
| Nz = Int(Lz/Δ) | ||
|
|
||
| grid = RegularCartesianGrid(size=(Nx, Ny, Nz), halo=(2, 2, 2), | ||
| x=(-Lx/2, Lx/2), y=(-Lx/2, Lx/2), z=(0, Lz)) | ||
|
|
||
| ##### | ||
| ##### Dry thermal bubble perturbation | ||
| ##### | ||
|
|
||
| xᶜ, zᶜ = 0km, 2km | ||
| xʳ, zʳ = 2km, 2km | ||
|
|
||
| @inline L(x, y, z) = √(((x - xᶜ) / xʳ)^2 + ((z - zᶜ) / zʳ)^2) | ||
| @inline function θ′(x, y, z) | ||
| ℓ = L(x, y, z) | ||
| return (ℓ <= 1) * 2cos(π/2 * ℓ)^2 | ||
| end | ||
|
|
||
| ##### | ||
| ##### Set up model | ||
| ##### | ||
|
|
||
| gas = IdealGas() | ||
| Rᵈ, cₚ, cᵥ = gas.Rᵈ, gas.cₚ, gas.cᵥ | ||
|
|
||
| g = 9.80665 | ||
| pₛ = 1000hPa | ||
| θₛ = 300 | ||
| πₛ = 1 | ||
|
|
||
| θ₀(x, y, z) = θₛ | ||
| π₀(x, y, z) = πₛ - g*z / (cₚ*θₛ) | ||
| p₀(x, y, z) = pₛ * π₀(x, y, z)^(cₚ/Rᵈ) | ||
| ρ₀(x, y, z) = pₛ / (Rᵈ * θ₀(x, y, z)) * π₀(x, y, z)^(cᵥ/Rᵈ) | ||
| Θ₀(x, y, z) = ρ₀(x, y, z) * θ₀(x, y, z) | ||
|
|
||
| const τ⁻¹ = 1 # Damping/relaxation time scale [s⁻¹]. This is very strong damping. | ||
| const Δμ = 0.1Lz # Sponge layer width [m] set to 10% of the domain height. | ||
| @inline μ(z, Lz) = τ⁻¹ * exp(-(Lz-z) / Δμ) | ||
|
|
||
| @inline Fw(i, j, k, grid, t, Ũ, C̃, p) = @inbounds (t <= 500) * -μ(grid.zF[k], grid.Lz) * Ũ.ρw[i, j, k] | ||
| forcing = ModelForcing(w=Fw) | ||
|
|
||
| model = CompressibleModel(grid=grid, buoyancy=gas, reference_pressure=pₛ, | ||
| prognostic_temperature=ModifiedPotentialTemperature(), | ||
| tracers=(:Θᵐ,), forcing=forcing) | ||
|
|
||
| set!(model.density, ρ₀) | ||
| set!(model.tracers.Θᵐ, Θ₀) | ||
|
|
||
| ##### | ||
| ##### Run a dry adiabatic atmosphere to hydrostatic balance | ||
| ##### | ||
|
|
||
| while model.clock.time < 500 | ||
| @printf("t = %.2f s\n", model.clock.time) | ||
| time_step!(model, Δt=0.2, Nt=100) | ||
| end | ||
|
|
||
| ##### | ||
| ##### Now add the cold bubble perturbation. | ||
| ##### | ||
|
|
||
| ρʰᵈ = model.density.data[1:Nx, 1, 1:Nz] | ||
| Θʰᵈ = model.tracers.Θᵐ.data[1:Nx, 1, 1:Nz] | ||
|
|
||
| xC, zC = grid.xC, grid.zC | ||
| ρ, Θ = model.density, model.tracers.Θᵐ | ||
| for k in 1:Nz, i in 1:Nx | ||
| θ = Θ[i, 1, k] / ρ[i, 1, k] + θ′(xC[i], 0, zC[k]) | ||
| π = Π(i, 1, k, grid, gas, Θ) | ||
|
|
||
| ρ[i, 1, k] = pₛ / (Rᵈ*θ) * π^(cᵥ/Rᵈ) | ||
| Θ[i, 1, k] = ρ[i, 1, k] * θ | ||
| end | ||
|
|
||
| ρ_plot = contour(model.grid.xC ./ km, model.grid.zC ./ km, rotr90(ρ.data[1:Nx, 1, 1:Nz] .- ρʰᵈ), | ||
| fill=true, levels=10, xlims=(-5, 5), clims=(-0.008, 0.008), color=:balance, dpi=200) | ||
| savefig(ρ_plot, "rho_prime_initial_condition.png") | ||
|
|
||
| θ_slice = rotr90(Θ.data[1:Nx, 1, 1:Nz] ./ ρ.data[1:Nx, 1, 1:Nz]) | ||
| Θ_plot = contour(model.grid.xC ./ km, model.grid.zC ./ km, θ_slice, | ||
| fill=true, levels=10, xlims=(-5, 5), color=:thermal, dpi=200) | ||
| savefig(Θ_plot, "theta_initial_condition.png") | ||
|
|
||
| ##### | ||
| ##### Watch the thermal bubble rise! | ||
| ##### | ||
|
|
||
| for n in 1:200 | ||
| time_step!(model, Δt=0.1, Nt=50) | ||
|
|
||
| @printf("t = %.2f s\n", model.clock.time) | ||
| xC, yC, zC = model.grid.xC ./ km, model.grid.yC ./ km, model.grid.zC ./ km | ||
| xF, yF, zF = model.grid.xF ./ km, model.grid.yF ./ km, model.grid.zF ./ km | ||
|
|
||
| j = 1 | ||
| u_slice = rotr90(model.momenta.ρu.data[1:Nx, j, 1:Nz] ./ model.density.data[1:Nx, j, 1:Nz]) | ||
| w_slice = rotr90(model.momenta.ρw.data[1:Nx, j, 1:Nz] ./ model.density.data[1:Nx, j, 1:Nz]) | ||
| ρ_slice = rotr90(model.density.data[1:Nx, j, 1:Nz] .- ρʰᵈ) | ||
| θ_slice = rotr90(model.tracers.Θᵐ.data[1:Nx, j, 1:Nz] ./ model.density.data[1:Nx, j, 1:Nz]) | ||
|
|
||
| u_title = @sprintf("u, t = %d s", round(Int, model.clock.time)) | ||
| pu = contour(xC, zC, u_slice, title=u_title, fill=true, levels=10, xlims=(-5, 5), color=:balance, clims=(-10, 10)) | ||
| pw = contour(xC, zC, w_slice, title="w", fill=true, levels=10, xlims=(-5, 5), color=:balance, clims=(-10, 10)) | ||
| pρ = contour(xC, zC, ρ_slice, title="rho_prime", fill=true, levels=10, xlims=(-5, 5), color=:balance, clims=(-0.006, 0.006)) | ||
| pθ = contour(xC, zC, θ_slice, title="theta", fill=true, levels=10, xlims=(-5, 5), color=:thermal, clims=(299.9, 302)) | ||
|
|
||
| p = plot(pu, pw, pρ, pθ, layout=(2, 2), dpi=200, show=true) | ||
| savefig(p, @sprintf("thermal_bubble_%03d.png", n)) | ||
| end | ||
|
|
||
| θ_1000 = (model.tracers.Θᵐ.data[1:Nx, 1, 1:Nz] ./ model.density.data[1:Nx, 1, 1:Nz]) .- θₛ | ||
| w_1000 = (model.momenta.ρw.data[1:Nx, 1, 1:Nz] ./ model.density.data[1:Nx, 1, 1:Nz]) | ||
|
|
||
| @printf("θ′: min=%.2f, max=%.2f\n", minimum(θ_1000), maximum(θ_1000)) | ||
| @printf("w: min=%.2f, max=%.2f\n", minimum(w_1000), maximum(w_1000)) |