From d09f649731f39a281279620c2cbc1972c079d758 Mon Sep 17 00:00:00 2001
From: Ali Ramadhan <ali.hh.ramadhan@gmail.com>
Date: Tue, 29 Jan 2019 17:42:51 -0500
Subject: [PATCH] Rayleigh Taylor instability (2D example).

---
 examples/rayleigh_taylor_2d.jl | 65 ++++++++++++++++++++++++++++++++++
 1 file changed, 65 insertions(+)
 create mode 100644 examples/rayleigh_taylor_2d.jl

diff --git a/examples/rayleigh_taylor_2d.jl b/examples/rayleigh_taylor_2d.jl
new file mode 100644
index 0000000000..3b2d25825c
--- /dev/null
+++ b/examples/rayleigh_taylor_2d.jl
@@ -0,0 +1,65 @@
+using Plots
+using Oceananigans
+
+function make_movies(problem::Problem, R::SavedFields, Nt, Δt)
+    g = problem.g
+
+    print("Creating tracer movie... ($(Int(Nt/R.ΔR)) frames)\n")
+
+    Plots.gr()
+
+    # animU = @animate for tidx in 1:Int(Nt/ΔR)
+    #     print("\rframe = $tidx / $(Int(Nt/ΔR))   ")
+    #     Plots.heatmap(g.xC ./ 1000, g.zC ./ 1000, rotl90(R.u[tidx, :, 1, :]), color=:balance,
+    #                   clims=(-0.01, 0.01),
+    #                   title="u-velocity @ t=$(tidx*ΔR*Δt)")
+    # end
+    # mp4(animU, "uvel_$(round(Int, time())).mp4", fps = 30)
+
+    # animT = @animate for tidx in 1:Int(Nt/R.ΔR)
+    #     print("\rframe = $tidx / $(Int(Nt/R.ΔR))   ")
+    #     Plots.heatmap(g.xC ./ 1000, g.zC ./ 1000, rotl90(R.T[tidx, :, 1, :]) .- 283, color=:balance,
+    #                   clims=(-0.02, 0.02),
+    #                   title="T change @ t=$(tidx*R.ΔR*Δt)")
+    # end
+    # mp4(animT, "rayleigh_taylor_2d_$(round(Int, time())).mp4", fps = 30)
+
+    default(dpi=300)
+
+    animT = @animate for tidx in 1:Int(Nt/R.ΔR)
+        print("\rframe = $tidx / $(Int(Nt/R.ΔR))   ")
+        Plots.contourf(1:g.Nx, 1:g.Nz, rotl90(R.T[tidx, :, 1, :]) .- 283, color=:balance,
+                      clims=(-0.02, 0.02), title="T change @ t=$(tidx*R.ΔR*Δt)")
+    end
+    mp4(animT, "rayleigh_taylor_2d_$(round(Int, time())).mp4", fps = 30)
+
+    # animρ = @animate for tidx in 1:Int(Nt/ΔR)
+    #     print("\rframe = $tidx / $(Int(Nt/ΔR))   ")
+    #     Plots.heatmap(g.xC ./ 1000, g.zC ./ 1000, rotl90(R.ρ[tidx, :, 1, :]) .- eos.ρ₀, color=:balance,
+    #                   clims=(-0.001, 0.001),
+    #                   # clims=(-maximum(R.ρ[tidx, :, 1, :] .- eos.ρ₀), maximum(R.ρ[tidx, :, 1, :] .- eos.ρ₀)),
+    #                   title="delta rho @ t=$(tidx*ΔR*Δt)")
+    # end
+    # mp4(animρ, "tracer_δρ_$(round(Int, time())).mp4", fps = 30)
+end
+
+function rayleigh_taylor_2d()
+    Nx, Ny, Nz = 512, 1, 256
+    Lx, Ly, Lz = 10000, 1, 5000
+    Nt, Δt = 10000, 10
+    ΔR = 10
+
+    problem = Problem((Nx, Ny, Nz), (Lx, Ly, Lz))
+
+    ΔT = 0.02  # Temperature difference.
+    ε = 0.0001  # Small temperature perturbation.
+
+    @. problem.tr.T.data[:, :, 1:Int(problem.g.Nz/2)]   = 283 - (ΔT/2) + ε*rand()
+    @. problem.tr.T.data[:, :, Int(problem.g.Nz/2):end] = 283 + (ΔT/2) + ε*rand()
+
+    R  = SavedFields(problem.g, Nt, ΔR)
+
+    time_stepping!(problem; Nt=Nt, Δt=Δt, R=R)
+
+    make_movies(problem, R, Nt, Δt)
+end