From eb57ff151ada34557e33487a12c33799c110388c Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 14 Oct 2019 19:15:43 -0400 Subject: [PATCH 01/31] Drafting a lid-driven cavity verification Former-commit-id: 9b22bc94e11955802293ac1f1c2c71972d7f1415 --- examples/lid_driven_cavity.jl | 55 +++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 examples/lid_driven_cavity.jl diff --git a/examples/lid_driven_cavity.jl b/examples/lid_driven_cavity.jl new file mode 100644 index 0000000000..704addb3b4 --- /dev/null +++ b/examples/lid_driven_cavity.jl @@ -0,0 +1,55 @@ +using Oceananigans, PyPlot, Printf +using Oceananigans: NoPenetrationBC, NonDimensionalModel + +Nx, Ny, Nz = 1, 64, 64 +Lx, Ly, Lz = 1, 1, 1 + +vbcs = ChannelBCs(top = BoundaryCondition(Value, 1), + bottom = BoundaryCondition(Value, 0), + north = NoPenetrationBC(), + south = NoPenetrationBC()) + +wbcs = ChannelBCs(top = NoPenetrationBC(), + bottom = NoPenetrationBC(), + north = BoundaryCondition(Value, 0), + south = BoundaryCondition(Value, 0)) + +bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) + +Re = 400 +model = NonDimensionalModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), Re=Re, Ri=0, Pr=Inf, Ro=Inf, boundary_conditions=bcs) + +fig, ax = subplots(nrows=1, ncols=1, figsize=(10, 10)) + +Δt = 1e-3 +Δ = max(model.grid.Δy, model.grid.Δz) +CFL = Δt / Δ +dCFL = (1/Re) * Δt / Δ^2 +@show CFL +@show dCFL + +while model.clock.time < 5 + time_step!(model; Δt=Δt, Nt=100, init_with_euler=model.clock.time == 0 ? true : false) + @printf("Time: %.4f\n", model.clock.time) + + y = collect(model.grid.yC) + z = collect(model.grid.zC) + v = model.velocities.v.data[1, :, :] + w = model.velocities.w.data[1, :, :] + + Δy, Δz = model.grid.Δy, model.grid.Δz + dvdz = (v[1:Ny, 2:Nz+1] - v[1:Ny, 1:Nz])/ Δz + dwdy = (w[1:Ny, 1:Nz] - w[2:Ny+1, 1:Nz])/ Δy + ζ = dwdy - dvdz + ζ = reverse(log10.(abs.(ζ)), dims=1) + + # ax.streamplot(y, z, v, w) + pcolormesh(y, z, ζ, vmin=1e-2, cmap="inferno") + # fig.colorbar(im, ax=ax) + + ax.set_title(@sprintf("Time: %.4f", model.clock.time)) + ax.set_xlabel("\$y\$"); ax.set_ylabel("\$z\$"); + ax.set_xlim([0, 1]); ax.set_ylim([-1, 0]); + ax.set_aspect(1) + gcf(); +end From 7948a803e3dfc91f38acb3a6b65b94ea0c495b05 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Tue, 15 Oct 2019 06:33:15 -0400 Subject: [PATCH 02/31] Move to verification folder and bump up to 128^2 Former-commit-id: 6ae31ced0957a2825d253c33361e61b6d17f322c --- .../lid_driven_cavity}/lid_driven_cavity.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) rename {examples => verification/lid_driven_cavity}/lid_driven_cavity.jl (85%) diff --git a/examples/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl similarity index 85% rename from examples/lid_driven_cavity.jl rename to verification/lid_driven_cavity/lid_driven_cavity.jl index 704addb3b4..f212ca2aaf 100644 --- a/examples/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,7 +1,7 @@ using Oceananigans, PyPlot, Printf using Oceananigans: NoPenetrationBC, NonDimensionalModel -Nx, Ny, Nz = 1, 64, 64 +Nx, Ny, Nz = 1, 128, 128 Lx, Ly, Lz = 1, 1, 1 vbcs = ChannelBCs(top = BoundaryCondition(Value, 1), @@ -21,14 +21,14 @@ model = NonDimensionalModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), Re=Re, Ri=0, Pr=Inf, fig, ax = subplots(nrows=1, ncols=1, figsize=(10, 10)) -Δt = 1e-3 +Δt = 2e-4 Δ = max(model.grid.Δy, model.grid.Δz) CFL = Δt / Δ dCFL = (1/Re) * Δt / Δ^2 @show CFL @show dCFL -while model.clock.time < 5 +# while model.clock.time < 1 time_step!(model; Δt=Δt, Nt=100, init_with_euler=model.clock.time == 0 ? true : false) @printf("Time: %.4f\n", model.clock.time) @@ -44,12 +44,12 @@ while model.clock.time < 5 ζ = reverse(log10.(abs.(ζ)), dims=1) # ax.streamplot(y, z, v, w) - pcolormesh(y, z, ζ, vmin=1e-2, cmap="inferno") + pcolormesh(y, z, ζ, vmin=1e-3, cmap="viridis"); title(); xlabel("y"); ylabel("z"); # fig.colorbar(im, ax=ax) - ax.set_title(@sprintf("Time: %.4f", model.clock.time)) + ax.set_title(@sprintf("Lid-driven cavity vorticity log₁₀(ζ): Re=%d, t=%.4f", Re, model.clock.time)) ax.set_xlabel("\$y\$"); ax.set_ylabel("\$z\$"); ax.set_xlim([0, 1]); ax.set_ylim([-1, 0]); ax.set_aspect(1) gcf(); -end +# end From 2e05c278709daddd61a259ab0c8dcdc27398b5a0 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Tue, 15 Oct 2019 07:55:32 -0400 Subject: [PATCH 03/31] Start comparing with Ghia et al. (1982) Former-commit-id: d98280ecbef4e960709c0918b6a58bce2c522f3f --- .../lid_driven_cavity/lid_driven_cavity.jl | 28 ++++++++++++++++--- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index f212ca2aaf..34ce646da3 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,6 +1,22 @@ using Oceananigans, PyPlot, Printf using Oceananigans: NoPenetrationBC, NonDimensionalModel +#### +#### Data from tables 1 and 2 of Ghia et al. (1982). +#### + +j̃ = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 156, 126, 129] +ỹ = [0.0, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0] + +ũ = Dict( + 100 => [0.0, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0], + 400 => [0.0, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0] +) + +#### +#### Model setup +#### + Nx, Ny, Nz = 1, 128, 128 Lx, Ly, Lz = 1, 1, 1 @@ -19,16 +35,20 @@ bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) Re = 400 model = NonDimensionalModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), Re=Re, Ri=0, Pr=Inf, Ro=Inf, boundary_conditions=bcs) +nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) +push!(model.diagnostics, nan_checker) + fig, ax = subplots(nrows=1, ncols=1, figsize=(10, 10)) -Δt = 2e-4 +Δt = 0.5e-4 Δ = max(model.grid.Δy, model.grid.Δz) CFL = Δt / Δ dCFL = (1/Re) * Δt / Δ^2 @show CFL @show dCFL -# while model.clock.time < 1 +while model.clock.time < 5 + # Δt = model.clock.time < 0.3 ? 0.5e-4 : 5e-4 time_step!(model; Δt=Δt, Nt=100, init_with_euler=model.clock.time == 0 ? true : false) @printf("Time: %.4f\n", model.clock.time) @@ -44,7 +64,7 @@ dCFL = (1/Re) * Δt / Δ^2 ζ = reverse(log10.(abs.(ζ)), dims=1) # ax.streamplot(y, z, v, w) - pcolormesh(y, z, ζ, vmin=1e-3, cmap="viridis"); title(); xlabel("y"); ylabel("z"); + pcolormesh(y, z, ζ, vmin=1e-3, cmap="viridis") # fig.colorbar(im, ax=ax) ax.set_title(@sprintf("Lid-driven cavity vorticity log₁₀(ζ): Re=%d, t=%.4f", Re, model.clock.time)) @@ -52,4 +72,4 @@ dCFL = (1/Re) * Δt / Δ^2 ax.set_xlim([0, 1]); ax.set_ylim([-1, 0]); ax.set_aspect(1) gcf(); -# end +end From fd113ad34a99d3f4d7adb3f31b9e0839d782d184 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 18:07:03 -0500 Subject: [PATCH 04/31] Update lid-driven cavity to v0.17.0 syntax and switch to Plots.jl --- .../lid_driven_cavity/lid_driven_cavity.jl | 57 +++++++++---------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 34ce646da3..a3b7160fed 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,5 +1,7 @@ -using Oceananigans, PyPlot, Printf -using Oceananigans: NoPenetrationBC, NonDimensionalModel +using Oceananigans +using Plots, Printf + +using Oceananigans: NoPenetrationBC #### #### Data from tables 1 and 2 of Ghia et al. (1982). @@ -33,43 +35,38 @@ wbcs = ChannelBCs(top = NoPenetrationBC(), bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) Re = 400 -model = NonDimensionalModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), Re=Re, Ri=0, Pr=Inf, Ro=Inf, boundary_conditions=bcs) - -nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) -push!(model.diagnostics, nan_checker) -fig, ax = subplots(nrows=1, ncols=1, figsize=(10, 10)) +model = NonDimensionalModel(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)), + Re=Re, Pr=Inf, Ro=Inf, tracers=nothing, buoyancy=nothing, boundary_conditions=bcs) -Δt = 0.5e-4 Δ = max(model.grid.Δy, model.grid.Δz) -CFL = Δt / Δ -dCFL = (1/Re) * Δt / Δ^2 -@show CFL -@show dCFL - -while model.clock.time < 5 - # Δt = model.clock.time < 0.3 ? 0.5e-4 : 5e-4 - time_step!(model; Δt=Δt, Nt=100, init_with_euler=model.clock.time == 0 ? true : false) - @printf("Time: %.4f\n", model.clock.time) - - y = collect(model.grid.yC) - z = collect(model.grid.zC) + +y = collect(model.grid.yC) +z = collect(model.grid.zC) +# p = heatmap(y, z, zeros(Ny, Nz), color=:viridis, clims=(1e-3, 1), show=true) + +Δt = 0.1e-4 + +# wizard = TimeStepWizard(cfl=0.2, Δt=0.5e-4, max_change=1.1, max_Δt=5.0) + +while model.clock.time < 1e-2 + time_step!(model; Δt=Δt, Nt=1, init_with_euler=model.clock.time == 0 ? true : false) + v = model.velocities.v.data[1, :, :] w = model.velocities.w.data[1, :, :] Δy, Δz = model.grid.Δy, model.grid.Δz - dvdz = (v[1:Ny, 2:Nz+1] - v[1:Ny, 1:Nz])/ Δz - dwdy = (w[1:Ny, 1:Nz] - w[2:Ny+1, 1:Nz])/ Δy + dvdz = (v[1:Ny, 2:Nz+1] - v[1:Ny, 1:Nz]) / Δz + dwdy = (w[2:Ny+1, 1:Nz] - w[1:Ny, 1:Nz]) / Δy ζ = dwdy - dvdz ζ = reverse(log10.(abs.(ζ)), dims=1) - # ax.streamplot(y, z, v, w) - pcolormesh(y, z, ζ, vmin=1e-3, cmap="viridis") - # fig.colorbar(im, ax=ax) + # heatmap!(p, y, z, ζ, color=:viridis, clims=(-3, 1), show=true) - ax.set_title(@sprintf("Lid-driven cavity vorticity log₁₀(ζ): Re=%d, t=%.4f", Re, model.clock.time)) - ax.set_xlabel("\$y\$"); ax.set_ylabel("\$z\$"); - ax.set_xlim([0, 1]); ax.set_ylim([-1, 0]); - ax.set_aspect(1) - gcf(); + u, v, w = model.velocities + u_max = max(maximum(v.data), maximum(w.data)) + CFL = u_max * Δt / Δ + dCFL = (1/Re) * Δt / Δ^2 + @printf("Time: %.4f, CFL: %.3g, dCFL: %.3g, max (v, w, ζ): %.2g, %.2g, %.2g\n", + model.clock.time, CFL, dCFL, maximum(v.data.parent), maximum(w.data.parent), maximum(ζ)) end From aeae31e9467a54481d20974c45928d94d1622fc7 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 18:38:36 -0500 Subject: [PATCH 05/31] Adaptive time-stepping did not help. --- .../lid_driven_cavity/lid_driven_cavity.jl | 40 ++++++++++++------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index a3b7160fed..7ecaebc4f5 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -2,6 +2,12 @@ using Oceananigans using Plots, Printf using Oceananigans: NoPenetrationBC +using Oceananigans.Diagnostics: AdvectiveCFL + +# Workaround for plotting many frames. +# See: https://github.com/JuliaPlots/Plots.jl/issues/1723 +# import GR +# GR.inline("png") #### #### Data from tables 1 and 2 of Ghia et al. (1982). @@ -37,20 +43,24 @@ bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) Re = 400 model = NonDimensionalModel(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)), - Re=Re, Pr=Inf, Ro=Inf, tracers=nothing, buoyancy=nothing, boundary_conditions=bcs) + Re=Re, Pr=Inf, Ro=Inf, + coriolis=nothing, tracers=nothing, buoyancy=nothing, boundary_conditions=bcs) Δ = max(model.grid.Δy, model.grid.Δz) y = collect(model.grid.yC) z = collect(model.grid.zC) -# p = heatmap(y, z, zeros(Ny, Nz), color=:viridis, clims=(1e-3, 1), show=true) +# p = heatmap(y, z, zeros(Ny, Nz), color=:viridis, show=true) + +# Δt = 0.5e-4 -Δt = 0.1e-4 +wizard = TimeStepWizard(cfl=0.2, Δt=1e-5, max_change=1.1, max_Δt=1e-2) +cfl = AdvectiveCFL(wizard) -# wizard = TimeStepWizard(cfl=0.2, Δt=0.5e-4, max_change=1.1, max_Δt=5.0) +while model.clock.time < 5e-3 + update_Δt!(wizard, model) -while model.clock.time < 1e-2 - time_step!(model; Δt=Δt, Nt=1, init_with_euler=model.clock.time == 0 ? true : false) + time_step!(model; Δt=wizard.Δt, Nt=1, init_with_euler=model.clock.time == 0 ? true : false) v = model.velocities.v.data[1, :, :] w = model.velocities.w.data[1, :, :] @@ -59,14 +69,16 @@ while model.clock.time < 1e-2 dvdz = (v[1:Ny, 2:Nz+1] - v[1:Ny, 1:Nz]) / Δz dwdy = (w[2:Ny+1, 1:Nz] - w[1:Ny, 1:Nz]) / Δy ζ = dwdy - dvdz - ζ = reverse(log10.(abs.(ζ)), dims=1) - - # heatmap!(p, y, z, ζ, color=:viridis, clims=(-3, 1), show=true) + ζ = log10.(abs.(ζ)) u, v, w = model.velocities - u_max = max(maximum(v.data), maximum(w.data)) - CFL = u_max * Δt / Δ - dCFL = (1/Re) * Δt / Δ^2 - @printf("Time: %.4f, CFL: %.3g, dCFL: %.3g, max (v, w, ζ): %.2g, %.2g, %.2g\n", - model.clock.time, CFL, dCFL, maximum(v.data.parent), maximum(w.data.parent), maximum(ζ)) + + # heatmap!(p, y, z, ζ, color=:viridis, show=true) + # heatmap!(p, y, z, interior(v)[1, :, :], color=:viridis, show=true) + + u_max = max(maximum(interior(v)), maximum(interior(w))) + CFL = cfl(model) + dCFL = (1/Re) * wizard.Δt / Δ^2 + @printf("Time: %.4g, Δt: %.4g CFL: %.3g, dCFL: %.3g, max (v, w, ζ): %.2g, %.2g, %.2g\n", + model.clock.time, wizard.Δt, CFL, dCFL, maximum(v.data.parent), maximum(w.data.parent), 10^maximum(ζ)) end From 20cc66b15d1c5bf7e3c8dd46ee7bd1e1e17ac5a1 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 21:21:56 -0500 Subject: [PATCH 06/31] Vortical initial condition did not help --- .../lid_driven_cavity/lid_driven_cavity.jl | 35 +++++++++++-------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 7ecaebc4f5..e3e5791f59 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -28,23 +28,28 @@ ũ = Dict( Nx, Ny, Nz = 1, 128, 128 Lx, Ly, Lz = 1, 1, 1 -vbcs = ChannelBCs(top = BoundaryCondition(Value, 1), - bottom = BoundaryCondition(Value, 0), +vbcs = ChannelBCs(top = BoundaryCondition(Value, 0.0), + bottom = BoundaryCondition(Value, 0.0), north = NoPenetrationBC(), south = NoPenetrationBC()) wbcs = ChannelBCs(top = NoPenetrationBC(), bottom = NoPenetrationBC(), - north = BoundaryCondition(Value, 0), - south = BoundaryCondition(Value, 0)) + north = BoundaryCondition(Value, 0.0), + south = BoundaryCondition(Value, 0.0)) bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) -Re = 400 +grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz)) -model = NonDimensionalModel(grid=RegularCartesianGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz)), - Re=Re, Pr=Inf, Ro=Inf, - coriolis=nothing, tracers=nothing, buoyancy=nothing, boundary_conditions=bcs) +Re = 100 +model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, + coriolis=nothing, tracers=nothing, buoyancy=nothing, + boundary_conditions=bcs) + +v₀(x, y, z) = -1 + 2z +w₀(x, y, z) = 1 - 2y +set!(model, v=v₀, w=w₀) Δ = max(model.grid.Δy, model.grid.Δz) @@ -54,13 +59,13 @@ z = collect(model.grid.zC) # Δt = 0.5e-4 -wizard = TimeStepWizard(cfl=0.2, Δt=1e-5, max_change=1.1, max_Δt=1e-2) +wizard = TimeStepWizard(cfl=0.01, Δt=1e-7, max_change=1.1, max_Δt=1e-6) cfl = AdvectiveCFL(wizard) -while model.clock.time < 5e-3 +while model.clock.time < 4e-4 update_Δt!(wizard, model) - time_step!(model; Δt=wizard.Δt, Nt=1, init_with_euler=model.clock.time == 0 ? true : false) + time_step!(model; Δt=wizard.Δt, Nt=10, init_with_euler=model.clock.time == 0 ? true : false) v = model.velocities.v.data[1, :, :] w = model.velocities.w.data[1, :, :] @@ -76,9 +81,11 @@ while model.clock.time < 5e-3 # heatmap!(p, y, z, ζ, color=:viridis, show=true) # heatmap!(p, y, z, interior(v)[1, :, :], color=:viridis, show=true) - u_max = max(maximum(interior(v)), maximum(interior(w))) + v_max = maximum(abs, interior(v)) + w_max = maximum(abs, interior(w)) CFL = cfl(model) dCFL = (1/Re) * wizard.Δt / Δ^2 - @printf("Time: %.4g, Δt: %.4g CFL: %.3g, dCFL: %.3g, max (v, w, ζ): %.2g, %.2g, %.2g\n", - model.clock.time, wizard.Δt, CFL, dCFL, maximum(v.data.parent), maximum(w.data.parent), 10^maximum(ζ)) + @printf("Time: %1.2e, Δt: %1.2e CFL: %1.2e, dCFL: %1.2e, max (v, w, ζ): %1.2e, %1.2e, %1.2e\n", + model.clock.time, wizard.Δt, CFL, dCFL, v_max, w_max, 10^maximum(ζ)) + # @printf("(u, v) @ corner: %.2e, %.2e\n", ) end From 1d318e2d42b7c51c55b34cbe0ed4101dbb45dab1 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 21:36:03 -0500 Subject: [PATCH 07/31] Slowly ramping up v_top still blows up --- .../lid_driven_cavity/lid_driven_cavity.jl | 24 ++++++++++++------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index e3e5791f59..dba01eb65f 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -2,7 +2,7 @@ using Oceananigans using Plots, Printf using Oceananigans: NoPenetrationBC -using Oceananigans.Diagnostics: AdvectiveCFL +using Oceananigans.Diagnostics # Workaround for plotting many frames. # See: https://github.com/JuliaPlots/Plots.jl/issues/1723 @@ -28,7 +28,7 @@ ũ = Dict( Nx, Ny, Nz = 1, 128, 128 Lx, Ly, Lz = 1, 1, 1 -vbcs = ChannelBCs(top = BoundaryCondition(Value, 0.0), +vbcs = ChannelBCs(top = BoundaryCondition(Value, 1.0), bottom = BoundaryCondition(Value, 0.0), north = NoPenetrationBC(), south = NoPenetrationBC()) @@ -47,9 +47,8 @@ model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, coriolis=nothing, tracers=nothing, buoyancy=nothing, boundary_conditions=bcs) -v₀(x, y, z) = -1 + 2z -w₀(x, y, z) = 1 - 2y -set!(model, v=v₀, w=w₀) +nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) +push!(model.diagnostics, nan_checker) Δ = max(model.grid.Δy, model.grid.Δz) @@ -59,13 +58,19 @@ z = collect(model.grid.zC) # Δt = 0.5e-4 -wizard = TimeStepWizard(cfl=0.01, Δt=1e-7, max_change=1.1, max_Δt=1e-6) +wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=1e-4) cfl = AdvectiveCFL(wizard) -while model.clock.time < 4e-4 +v_top(t) = min(1, t) + +while model.clock.time < 1e-3 + t = model.clock.time + update_Δt!(wizard, model) - time_step!(model; Δt=wizard.Δt, Nt=10, init_with_euler=model.clock.time == 0 ? true : false) + model.boundary_conditions.solution.v.z.top = BoundaryCondition(Value, v_top(t)) + + time_step!(model; Δt=wizard.Δt, Nt=10, init_with_euler = t == 0 ? true : false) v = model.velocities.v.data[1, :, :] w = model.velocities.w.data[1, :, :] @@ -87,5 +92,6 @@ while model.clock.time < 4e-4 dCFL = (1/Re) * wizard.Δt / Δ^2 @printf("Time: %1.2e, Δt: %1.2e CFL: %1.2e, dCFL: %1.2e, max (v, w, ζ): %1.2e, %1.2e, %1.2e\n", model.clock.time, wizard.Δt, CFL, dCFL, v_max, w_max, 10^maximum(ζ)) - # @printf("(u, v) @ corner: %.2e, %.2e\n", ) + + @show model.boundary_conditions.solution.v.z.top.condition end From 780a86455b925822770ec389d9553eb39305146e Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 21:53:03 -0500 Subject: [PATCH 08/31] Varying relaxation timescale did not help --- .../lid_driven_cavity/lid_driven_cavity.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index dba01eb65f..df5467d777 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -40,12 +40,28 @@ wbcs = ChannelBCs(top = NoPenetrationBC(), bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) +@inline function Fv(i, j, k, grid, time, U, C, p) + if k == Nz + @inbounds return - 2*1e-5/grid.Δz^2 * (U.v[i, j, Nz] - 1) + else + return 0 + end +end + +# @inline Fv(i, j, k, grid, time, U, C, p) = +# @inbounds ifelse(k == 1, - 2/grid.Δz^2 * U.v[i, j, 1], 0) + +@inline Fw(i, j, k, grid, time, U, C, p) = + @inbounds ifelse(j == 1 || j == Nz, - 2/grid.Δy^2 * U.w[i, j, 1], 0) + +forcing = ModelForcing(v=Fv) + grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz)) Re = 100 model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, coriolis=nothing, tracers=nothing, buoyancy=nothing, - boundary_conditions=bcs) + boundary_conditions=bcs, forcing=forcing) nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) push!(model.diagnostics, nan_checker) From 376bbfedcb7f56e04bf9133d285b2c3a52ae5d6e Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 22:04:41 -0500 Subject: [PATCH 09/31] Strong relaxation for no-slip BC did not help --- .../lid_driven_cavity/lid_driven_cavity.jl | 27 ++++++++++++------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index df5467d777..bfb533f43d 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -40,21 +40,28 @@ wbcs = ChannelBCs(top = NoPenetrationBC(), bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) +# @inline Fv(i, j, k, grid, time, U, C, p) = +# @inbounds ifelse(k == 1, - 2/grid.Δz^2 * U.v[i, j, 1], 0) + @inline function Fv(i, j, k, grid, time, U, C, p) - if k == Nz - @inbounds return - 2*1e-5/grid.Δz^2 * (U.v[i, j, Nz] - 1) + if k == 1 + return @inbounds - 2/grid.Δz^2 * U.v[i, j, 1] else return 0 end end -# @inline Fv(i, j, k, grid, time, U, C, p) = -# @inbounds ifelse(k == 1, - 2/grid.Δz^2 * U.v[i, j, 1], 0) - -@inline Fw(i, j, k, grid, time, U, C, p) = - @inbounds ifelse(j == 1 || j == Nz, - 2/grid.Δy^2 * U.w[i, j, 1], 0) +@inline function Fw(i, j, k, grid, time, U, C, p) + if j == 1 + return @inbounds - 2/grid.Δy^2 * U.w[i, 1, k] + elseif j == grid.Ny + return @inbounds - 2/grid.Δy^2 * U.w[i, grid.Ny, k] + else + return 0 + end +end -forcing = ModelForcing(v=Fv) +forcing = ModelForcing(v=Fv, w=Fw) grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz)) @@ -74,12 +81,12 @@ z = collect(model.grid.zC) # Δt = 0.5e-4 -wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=1e-4) +wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=1e-5) cfl = AdvectiveCFL(wizard) v_top(t) = min(1, t) -while model.clock.time < 1e-3 +while model.clock.time < 4e-3 t = model.clock.time update_Δt!(wizard, model) From c881fea182a19b7385528c91989b9df6148bf9d0 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 13 Dec 2019 22:05:19 -0500 Subject: [PATCH 10/31] Noise did not help --- verification/lid_driven_cavity/lid_driven_cavity.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index bfb533f43d..3151b70300 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -73,6 +73,9 @@ model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) push!(model.diagnostics, nan_checker) +ε(x, y, z) = 1e-2 * randn() +set!(model, v=ε, w=ε) + Δ = max(model.grid.Δy, model.grid.Δz) y = collect(model.grid.yC) From 8eec0a99ac88e48a4ccb968c706145d12bf4836e Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sat, 14 Dec 2019 09:17:21 -0500 Subject: [PATCH 11/31] Clean up lid-driven cavity and switch to vorticity computation --- .../lid_driven_cavity/lid_driven_cavity.jl | 59 ++++--------------- 1 file changed, 13 insertions(+), 46 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 3151b70300..f7dfc89914 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,8 +1,10 @@ -using Oceananigans using Plots, Printf -using Oceananigans: NoPenetrationBC +using Oceananigans using Oceananigans.Diagnostics +using Oceananigans.AbstractOperations + +using Oceananigans: Face, Cell, NoPenetrationBC # Workaround for plotting many frames. # See: https://github.com/JuliaPlots/Plots.jl/issues/1723 @@ -40,29 +42,6 @@ wbcs = ChannelBCs(top = NoPenetrationBC(), bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) -# @inline Fv(i, j, k, grid, time, U, C, p) = -# @inbounds ifelse(k == 1, - 2/grid.Δz^2 * U.v[i, j, 1], 0) - -@inline function Fv(i, j, k, grid, time, U, C, p) - if k == 1 - return @inbounds - 2/grid.Δz^2 * U.v[i, j, 1] - else - return 0 - end -end - -@inline function Fw(i, j, k, grid, time, U, C, p) - if j == 1 - return @inbounds - 2/grid.Δy^2 * U.w[i, 1, k] - elseif j == grid.Ny - return @inbounds - 2/grid.Δy^2 * U.w[i, grid.Ny, k] - else - return 0 - end -end - -forcing = ModelForcing(v=Fv, w=Fw) - grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz)) Re = 100 @@ -73,10 +52,11 @@ model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) push!(model.diagnostics, nan_checker) -ε(x, y, z) = 1e-2 * randn() -set!(model, v=ε, w=ε) +u, v, w = model.velocities +ζ_op = ∂y(w) - ∂z(v) +ζ = Field(Face, Face, Cell, model.architecture, model.grid) +ζ_comp = Computation(ζ_op, ζ) -Δ = max(model.grid.Δy, model.grid.Δz) y = collect(model.grid.yC) z = collect(model.grid.zC) @@ -87,37 +67,24 @@ z = collect(model.grid.zC) wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=1e-5) cfl = AdvectiveCFL(wizard) -v_top(t) = min(1, t) - while model.clock.time < 4e-3 t = model.clock.time update_Δt!(wizard, model) - model.boundary_conditions.solution.v.z.top = BoundaryCondition(Value, v_top(t)) - time_step!(model; Δt=wizard.Δt, Nt=10, init_with_euler = t == 0 ? true : false) - v = model.velocities.v.data[1, :, :] - w = model.velocities.w.data[1, :, :] - - Δy, Δz = model.grid.Δy, model.grid.Δz - dvdz = (v[1:Ny, 2:Nz+1] - v[1:Ny, 1:Nz]) / Δz - dwdy = (w[2:Ny+1, 1:Nz] - w[1:Ny, 1:Nz]) / Δy - ζ = dwdy - dvdz - ζ = log10.(abs.(ζ)) + compute!(ζ_comp) - u, v, w = model.velocities - - # heatmap!(p, y, z, ζ, color=:viridis, show=true) + # heatmap!(p, y, z, interior(ζ)[1, :, :], color=:viridis, show=true) # heatmap!(p, y, z, interior(v)[1, :, :], color=:viridis, show=true) v_max = maximum(abs, interior(v)) w_max = maximum(abs, interior(w)) + ζ_max = maximum(abs, interior(ζ)) CFL = cfl(model) dCFL = (1/Re) * wizard.Δt / Δ^2 - @printf("Time: %1.2e, Δt: %1.2e CFL: %1.2e, dCFL: %1.2e, max (v, w, ζ): %1.2e, %1.2e, %1.2e\n", - model.clock.time, wizard.Δt, CFL, dCFL, v_max, w_max, 10^maximum(ζ)) - @show model.boundary_conditions.solution.v.z.top.condition + @printf("Time: %1.2e, Δt: %1.2e CFL: %1.2e, dCFL: %1.2e, max (v, w, ζ): %1.2e, %1.2e, %1.2e\n", + model.clock.time, wizard.Δt, CFL, dCFL, v_max, w_max, ζ_max) end From c8819fb9c8659da43d02ba37c2055fdb4fc111fd Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sun, 8 Mar 2020 23:31:58 -0400 Subject: [PATCH 12/31] Update Lid-driven cavity script to Oceananigans v0.25 --- .../lid_driven_cavity/lid_driven_cavity.jl | 115 +++++++++--------- 1 file changed, 57 insertions(+), 58 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index f7dfc89914..5c08253a23 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,19 +1,13 @@ -using Plots, Printf - +using Printf using Oceananigans +using Oceananigans: Face, Cell using Oceananigans.Diagnostics +using Oceananigans.OutputWriters using Oceananigans.AbstractOperations -using Oceananigans: Face, Cell, NoPenetrationBC - -# Workaround for plotting many frames. -# See: https://github.com/JuliaPlots/Plots.jl/issues/1723 -# import GR -# GR.inline("png") - -#### -#### Data from tables 1 and 2 of Ghia et al. (1982). -#### +##### +##### Data from tables 1 and 2 of Ghia et al. (1982). +##### j̃ = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 156, 126, 129] ỹ = [0.0, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0] @@ -23,68 +17,73 @@ ũ = Dict( 400 => [0.0, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0] ) -#### -#### Model setup -#### - -Nx, Ny, Nz = 1, 128, 128 -Lx, Ly, Lz = 1, 1, 1 +##### +##### Model setup +##### -vbcs = ChannelBCs(top = BoundaryCondition(Value, 1.0), - bottom = BoundaryCondition(Value, 0.0), - north = NoPenetrationBC(), - south = NoPenetrationBC()) +Ny = Nz = 128 +Ly = Lz = 1.0 -wbcs = ChannelBCs(top = NoPenetrationBC(), - bottom = NoPenetrationBC(), - north = BoundaryCondition(Value, 0.0), - south = BoundaryCondition(Value, 0.0)) +topology = (Flat, Bounded, Bounded) +domain = (x=(0, 1), y=(0, Ly), z=(0, Lz)) +grid = RegularCartesianGrid(topology=topology, size=(1, Ny, Nz); domain...) -bcs = ChannelSolutionBCs(v=vbcs, w=wbcs) +v_bcs = VVelocityBoundaryConditions(grid, + top = ValueBoundaryCondition(1.0), + bottom = ValueBoundaryCondition(0.0) +) -grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(0, Lz)) +w_bcs = WVelocityBoundaryConditions(grid, + north = ValueBoundaryCondition(0.0), + south = ValueBoundaryCondition(0.0) +) -Re = 100 -model = NonDimensionalModel(grid=grid, Re=Re, Pr=Inf, Ro=Inf, - coriolis=nothing, tracers=nothing, buoyancy=nothing, - boundary_conditions=bcs, forcing=forcing) +Re = 100 # Reynolds number -nan_checker = NaNChecker(model; frequency=10, fields=Dict(:w => model.velocities.w)) -push!(model.diagnostics, nan_checker) +model = IncompressibleModel( + grid=grid, + buoyancy=nothing, + tracers=nothing, + coriolis=nothing, + boundary_conditions = (v=v_bcs, w=w_bcs), + closure = ConstantIsotropicDiffusivity(ν=1/Re) +) u, v, w = model.velocities ζ_op = ∂y(w) - ∂z(v) -ζ = Field(Face, Face, Cell, model.architecture, model.grid) +ζ = Field(Cell, Face, Face, model.architecture, model.grid, TracerBoundaryConditions(grid)) ζ_comp = Computation(ζ_op, ζ) - -y = collect(model.grid.yC) -z = collect(model.grid.zC) -# p = heatmap(y, z, zeros(Ny, Nz), color=:viridis, show=true) - -# Δt = 0.5e-4 - -wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=1e-5) +max_Δt = 0.25 * model.grid.Δy^2 * Re / 2 +wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=max_Δt) cfl = AdvectiveCFL(wizard) +dcfl = DiffusiveCFL(wizard) -while model.clock.time < 4e-3 - t = model.clock.time +function print_progress(simulation) + model = simulation.model - update_Δt!(wizard, model) + # Calculate simulation progress in %. + progress = 100 * (model.clock.time / simulation.stop_time) - time_step!(model; Δt=wizard.Δt, Nt=10, init_with_euler = t == 0 ? true : false) + # Find maximum velocities. + vmax = maximum(abs, model.velocities.v.data.parent) + wmax = maximum(abs, model.velocities.w.data.parent) - compute!(ζ_comp) + i, t = model.clock.iteration, model.clock.time + @printf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e\n", + progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt) +end - # heatmap!(p, y, z, interior(ζ)[1, :, :], color=:viridis, show=true) - # heatmap!(p, y, z, interior(v)[1, :, :], color=:viridis, show=true) +simulation = Simulation(model, Δt=wizard, stop_time=10, progress=print_progress, progress_frequency=20) - v_max = maximum(abs, interior(v)) - w_max = maximum(abs, interior(w)) - ζ_max = maximum(abs, interior(ζ)) - CFL = cfl(model) - dCFL = (1/Re) * wizard.Δt / Δ^2 +fields = Dict("v" => model.velocities.v, "w" => model.velocities.w, "ζ" => model -> ζ_comp(model)) +dims = Dict("ζ" => ("xC", "yF", "zF")) +global_attributes = Dict("Re" => Re) +output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) - @printf("Time: %1.2e, Δt: %1.2e CFL: %1.2e, dCFL: %1.2e, max (v, w, ζ): %1.2e, %1.2e, %1.2e\n", - model.clock.time, wizard.Δt, CFL, dCFL, v_max, w_max, ζ_max) -end +simulation.output_writers[:fields] = + NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.1, + global_attributes=global_attributes, output_attributes=output_attributes, + dimensions=dims) + +run!(simulation) From aea2a233066e861d864b2df6da78ee285c472a65 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sun, 8 Mar 2020 23:32:11 -0400 Subject: [PATCH 13/31] Plotting lid-driven cavity vorticity --- .../lid_driven_cavity/plot_vorticity.py | 31 +++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 verification/lid_driven_cavity/plot_vorticity.py diff --git a/verification/lid_driven_cavity/plot_vorticity.py b/verification/lid_driven_cavity/plot_vorticity.py new file mode 100644 index 0000000000..38bf8adf55 --- /dev/null +++ b/verification/lid_driven_cavity/plot_vorticity.py @@ -0,0 +1,31 @@ +import joblib +import xarray as xr +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import cmocean +import ffmpeg + +def plot_lid_driven_cavity_vorticity(Re, n): + ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + + fig, ax = plt.subplots(figsize=(16, 9), dpi=200) + fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) + + # v = ds.v.isel(time=n).squeeze() + # v.plot.pcolormesh(ax=ax, vmin=-1, vmax=1, cmap=cmocean.cm.balance, extend="both") + + ζ = ds.ζ.isel(time=n).squeeze() + ζ.plot.pcolormesh(ax=ax, cmap=cmocean.cm.curl, extend="both", + norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + + ax.set_title("") + ax.set_xlabel("y") + ax.set_ylabel("z") + ax.set_aspect("equal") + + print(f"Saving frame {n}/{ds.time.size-1}...") + plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") + plt.close("all") + +ds = xr.open_dataset(f"lid_driven_cavity_Re100.nc") +plot_lid_driven_cavity_vorticity(100, ds.time.size-1) From a28ae76f89a04779a564f71958bed7353f5b2123 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 9 Mar 2020 21:37:00 -0400 Subject: [PATCH 14/31] Cleanup --- .../lid_driven_cavity/lid_driven_cavity.jl | 18 +--------- .../lid_driven_cavity/plot_vorticity.py | 34 ++++++++++++++----- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 5c08253a23..d33e2c82fb 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -5,22 +5,6 @@ using Oceananigans.Diagnostics using Oceananigans.OutputWriters using Oceananigans.AbstractOperations -##### -##### Data from tables 1 and 2 of Ghia et al. (1982). -##### - -j̃ = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 156, 126, 129] -ỹ = [0.0, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0] - -ũ = Dict( - 100 => [0.0, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0], - 400 => [0.0, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0] -) - -##### -##### Model setup -##### - Ny = Nz = 128 Ly = Lz = 1.0 @@ -74,7 +58,7 @@ function print_progress(simulation) progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt) end -simulation = Simulation(model, Δt=wizard, stop_time=10, progress=print_progress, progress_frequency=20) +simulation = Simulation(model, Δt=wizard, stop_time=100, progress=print_progress, progress_frequency=20) fields = Dict("v" => model.velocities.v, "w" => model.velocities.w, "ζ" => model -> ζ_comp(model)) dims = Dict("ζ" => ("xC", "yF", "zF")) diff --git a/verification/lid_driven_cavity/plot_vorticity.py b/verification/lid_driven_cavity/plot_vorticity.py index 38bf8adf55..479b9f1bcf 100644 --- a/verification/lid_driven_cavity/plot_vorticity.py +++ b/verification/lid_driven_cavity/plot_vorticity.py @@ -1,10 +1,23 @@ import joblib +import numpy as np import xarray as xr import matplotlib.pyplot as plt import matplotlib.colors as colors import cmocean import ffmpeg +##### +##### Data from tables 1 and 2 of Ghia et al. (1982). +##### + +j = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] +y = [0.0, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0] + +V = { + 100: [0.0, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0], + 400: [0.0, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0] +} + def plot_lid_driven_cavity_vorticity(Re, n): ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") @@ -14,18 +27,23 @@ def plot_lid_driven_cavity_vorticity(Re, n): # v = ds.v.isel(time=n).squeeze() # v.plot.pcolormesh(ax=ax, vmin=-1, vmax=1, cmap=cmocean.cm.balance, extend="both") - ζ = ds.ζ.isel(time=n).squeeze() - ζ.plot.pcolormesh(ax=ax, cmap=cmocean.cm.curl, extend="both", - norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + # ζ = ds.ζ.isel(time=n).squeeze() + # ζ.plot.pcolormesh(ax=ax, cmap=cmocean.cm.curl, extend="both", + # norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + + # ax.set_title("") + # ax.set_xlabel("y") + # ax.set_ylabel("z") + # ax.set_aspect("equal") + + plt.plot(y, V[Re], linestyle="", marker="o") - ax.set_title("") - ax.set_xlabel("y") - ax.set_ylabel("z") - ax.set_aspect("equal") + v = ds.v.isel(time=n, yF=64) + plt.plot(ds.zC, v.values.flatten()) print(f"Saving frame {n}/{ds.time.size-1}...") plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -ds = xr.open_dataset(f"lid_driven_cavity_Re100.nc") +ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") plot_lid_driven_cavity_vorticity(100, ds.time.size-1) From c519a3fdd59420c1ebfd0dbf3f8244a56d51115f Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 9 Mar 2020 22:33:23 -0400 Subject: [PATCH 15/31] Sweet velocity plots --- .../lid_driven_cavity/plot_vorticity.py | 65 ++++++++++++++----- 1 file changed, 47 insertions(+), 18 deletions(-) diff --git a/verification/lid_driven_cavity/plot_vorticity.py b/verification/lid_driven_cavity/plot_vorticity.py index 479b9f1bcf..96d9a9aefb 100644 --- a/verification/lid_driven_cavity/plot_vorticity.py +++ b/verification/lid_driven_cavity/plot_vorticity.py @@ -10,40 +10,69 @@ ##### Data from tables 1 and 2 of Ghia et al. (1982). ##### -j = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] -y = [0.0, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0] +j_Ghia = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] +y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000] -V = { - 100: [0.0, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0], - 400: [0.0, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0] +v_Ghia = { + 100: [0.0000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0000], + 400: [0.0000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0000] +} + +k_Ghia = [1, 9, 10, 11, 13, 21, 30, 31, 65, 104, 111, 117, 122, 123, 124, 125, 129] +z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000] + +w_Ghia = { + 100: [0.0000, 0.09233, 0.10091, 0.1089, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000] } def plot_lid_driven_cavity_vorticity(Re, n): ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") - fig, ax = plt.subplots(figsize=(16, 9), dpi=200) + fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 9), dpi=200) + plt.subplots_adjust(hspace=0.25) fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) - # v = ds.v.isel(time=n).squeeze() - # v.plot.pcolormesh(ax=ax, vmin=-1, vmax=1, cmap=cmocean.cm.balance, extend="both") + ax_v, ax_V = axes[0, 0], axes[0, 1] + ax_w, ax_W = axes[1, 0], axes[1, 1] - # ζ = ds.ζ.isel(time=n).squeeze() - # ζ.plot.pcolormesh(ax=ax, cmap=cmocean.cm.curl, extend="both", - # norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + v_line = ds.v.isel(time=n, yF=64) + ax_v.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none") + ax_v.plot(ds.zC, v_line.values.flatten(), label="Oceananigans.jl", color="tab:blue") + ax_v.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_v.set_xlabel("z") + ax_v.set_ylabel("v") + ax_v.set_xlim([0, 1]) - # ax.set_title("") - # ax.set_xlabel("y") - # ax.set_ylabel("z") - # ax.set_aspect("equal") + w_line = ds.w.isel(time=n, zF=64) + ax_w.plot(z_Ghia, w_Ghia[Re], label="Ghia et al. (1982)", color="tab:orange", linestyle="", marker="o", fillstyle="none") + ax_w.plot(ds.yC, w_line.values.flatten(), label="Oceananigans.jl", color="tab:orange") + ax_w.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_w.set_xlabel("y") + ax_w.set_ylabel("w") + ax_w.set_xlim([0, 1]) - plt.plot(y, V[Re], linestyle="", marker="o") + v = ds.v.isel(time=n).squeeze() + img_v = v.plot.pcolormesh(ax=ax_V, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) + fig.colorbar(img_v, ax=ax_V, extend="both") + ax_V.axvline(x=0.5, color="tab:blue", alpha=0.5) + ax_V.set_title("v-velocity") + ax_V.set_xlabel("y") + ax_V.set_ylabel("z") + ax_V.set_aspect("equal") - v = ds.v.isel(time=n, yF=64) - plt.plot(ds.zC, v.values.flatten()) + w = ds.w.isel(time=n).squeeze() + img_w = w.plot.pcolormesh(ax=ax_W, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) + fig.colorbar(img_w, ax=ax_W, extend="both") + ax_W.axhline(y=0.5, color="tab:orange", alpha=0.5) + ax_W.set_title("w-velocity") + ax_W.set_xlabel("y") + ax_W.set_ylabel("z") + ax_W.set_aspect("equal") print(f"Saving frame {n}/{ds.time.size-1}...") plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") +Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") plot_lid_driven_cavity_vorticity(100, ds.time.size-1) From 48bec1a0aa58211801b0c8272ea2c94e7ef021ce Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 9 Mar 2020 22:43:58 -0400 Subject: [PATCH 16/31] Cleanup --- .../plot_lid_driven_cavity.py | 86 +++++++++++++++++++ .../lid_driven_cavity/plot_vorticity.py | 78 ----------------- 2 files changed, 86 insertions(+), 78 deletions(-) create mode 100644 verification/lid_driven_cavity/plot_lid_driven_cavity.py delete mode 100644 verification/lid_driven_cavity/plot_vorticity.py diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py new file mode 100644 index 0000000000..8395fa8e86 --- /dev/null +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -0,0 +1,86 @@ +import joblib +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import cmocean +import ffmpeg + +from numpy import nan + +##### +##### Data from tables 1 and 2 of Ghia et al. (1982). +##### + +j_Ghia = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] +y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000] + +v_Ghia = { + 100: [0.0000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0000], + 400: [0.0000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0000] +} + +k_Ghia = [1, 9, 10, 11, 13, 21, 30, 31, 65, 104, 111, 117, 122, 123, 124, 125, 129] +z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000] + +w_Ghia = { + 100: [0.0000, 0.09233, 0.10091, 0.1089, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000] +} + +y_ζ_Ghia = [0.0000, 0.0625, 0.1250, 0.1875, 0.2500, 0.3125, 0.3750, 0.4375, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8125, 0.8750, 0.9375, 1.0000] + +ζ_Ghia = { + 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan] +} + +def plot_velocity_frame(Re, n): + ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + + fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 9), dpi=200) + plt.subplots_adjust(hspace=0.25) + fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) + + ax_v_line, ax_v_mesh = axes[0, 0], axes[0, 1] + ax_w_line, ax_w_mesh = axes[1, 0], axes[1, 1] + + v_line = ds.v.isel(time=n, yF=64) + ax_v_line.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none") + ax_v_line.plot(ds.zC, v_line.values.flatten(), label="Oceananigans.jl", color="tab:blue") + ax_v_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_v_line.set_xlabel("z") + ax_v_line.set_ylabel("v") + ax_v_line.set_xlim([0, 1]) + + w_line = ds.w.isel(time=n, zF=64) + ax_w_line.plot(z_Ghia, w_Ghia[Re], label="Ghia et al. (1982)", color="tab:orange", linestyle="", marker="o", fillstyle="none") + ax_w_line.plot(ds.yC, w_line.values.flatten(), label="Oceananigans.jl", color="tab:orange") + ax_w_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_w_line.set_xlabel("y") + ax_w_line.set_ylabel("w") + ax_w_line.set_xlim([0, 1]) + + v = ds.v.isel(time=n).squeeze() + img_v = v.plot.pcolormesh(ax=ax_v_mesh, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) + fig.colorbar(img_v, ax=ax_v_mesh, extend="both") + ax_v_mesh.axvline(x=0.5, color="tab:blue", alpha=0.5) + ax_v_mesh.set_title("v-velocity") + ax_v_mesh.set_xlabel("y") + ax_v_mesh.set_ylabel("z") + ax_v_mesh.set_aspect("equal") + + w = ds.w.isel(time=n).squeeze() + img_w = w.plot.pcolormesh(ax=ax_w_mesh, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) + fig.colorbar(img_w, ax=ax_w_mesh, extend="both") + ax_w_mesh.axhline(y=0.5, color="tab:orange", alpha=0.5) + ax_w_mesh.set_title("w-velocity") + ax_w_mesh.set_xlabel("y") + ax_w_mesh.set_ylabel("z") + ax_w_mesh.set_aspect("equal") + + print(f"Saving velocity frame {n}/{ds.time.size-1}...") + plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") + plt.close("all") + +Re = 100 +ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") +plot_lid_driven_cavity_vorticity(100, ds.time.size-1) diff --git a/verification/lid_driven_cavity/plot_vorticity.py b/verification/lid_driven_cavity/plot_vorticity.py deleted file mode 100644 index 96d9a9aefb..0000000000 --- a/verification/lid_driven_cavity/plot_vorticity.py +++ /dev/null @@ -1,78 +0,0 @@ -import joblib -import numpy as np -import xarray as xr -import matplotlib.pyplot as plt -import matplotlib.colors as colors -import cmocean -import ffmpeg - -##### -##### Data from tables 1 and 2 of Ghia et al. (1982). -##### - -j_Ghia = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] -y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000] - -v_Ghia = { - 100: [0.0000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0000], - 400: [0.0000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0000] -} - -k_Ghia = [1, 9, 10, 11, 13, 21, 30, 31, 65, 104, 111, 117, 122, 123, 124, 125, 129] -z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000] - -w_Ghia = { - 100: [0.0000, 0.09233, 0.10091, 0.1089, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000] -} - -def plot_lid_driven_cavity_vorticity(Re, n): - ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") - - fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 9), dpi=200) - plt.subplots_adjust(hspace=0.25) - fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) - - ax_v, ax_V = axes[0, 0], axes[0, 1] - ax_w, ax_W = axes[1, 0], axes[1, 1] - - v_line = ds.v.isel(time=n, yF=64) - ax_v.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none") - ax_v.plot(ds.zC, v_line.values.flatten(), label="Oceananigans.jl", color="tab:blue") - ax_v.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) - ax_v.set_xlabel("z") - ax_v.set_ylabel("v") - ax_v.set_xlim([0, 1]) - - w_line = ds.w.isel(time=n, zF=64) - ax_w.plot(z_Ghia, w_Ghia[Re], label="Ghia et al. (1982)", color="tab:orange", linestyle="", marker="o", fillstyle="none") - ax_w.plot(ds.yC, w_line.values.flatten(), label="Oceananigans.jl", color="tab:orange") - ax_w.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) - ax_w.set_xlabel("y") - ax_w.set_ylabel("w") - ax_w.set_xlim([0, 1]) - - v = ds.v.isel(time=n).squeeze() - img_v = v.plot.pcolormesh(ax=ax_V, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) - fig.colorbar(img_v, ax=ax_V, extend="both") - ax_V.axvline(x=0.5, color="tab:blue", alpha=0.5) - ax_V.set_title("v-velocity") - ax_V.set_xlabel("y") - ax_V.set_ylabel("z") - ax_V.set_aspect("equal") - - w = ds.w.isel(time=n).squeeze() - img_w = w.plot.pcolormesh(ax=ax_W, vmin=-1, vmax=1, cmap=cmocean.cm.balance, add_colorbar=False) - fig.colorbar(img_w, ax=ax_W, extend="both") - ax_W.axhline(y=0.5, color="tab:orange", alpha=0.5) - ax_W.set_title("w-velocity") - ax_W.set_xlabel("y") - ax_W.set_ylabel("z") - ax_W.set_aspect("equal") - - print(f"Saving frame {n}/{ds.time.size-1}...") - plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") - plt.close("all") - -Re = 100 -ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") -plot_lid_driven_cavity_vorticity(100, ds.time.size-1) From 6930d1ba56deb6873cfc22bdeb3630852938cc04 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 9 Mar 2020 22:55:11 -0400 Subject: [PATCH 17/31] Plotting vorticity --- .../plot_lid_driven_cavity.py | 36 +++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 8395fa8e86..11fab4e24c 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -78,9 +78,41 @@ def plot_velocity_frame(Re, n): ax_w_mesh.set_aspect("equal") print(f"Saving velocity frame {n}/{ds.time.size-1}...") - plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") + plt.savefig(f"lid_driven_cavity_velocity_Re{Re}_{n:05d}.png") + plt.close("all") + +def plot_vorticity_frame(Re, n): + ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + + fig, axes = plt.subplots(ncols=2, figsize=(16, 9), dpi=200) + fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) + + ax_line, ax_mesh = axes[0], axes[1] + + ζ_line = -ds.ζ.isel(time=n, zF=127) + ax_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") + ax_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:red") + ax_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_line.set_xlabel("y") + ax_line.set_ylabel("vorticity $\zeta$") + ax_line.set_xlim([0, 1]) + ax_line.set_ylim(bottom=0) + + ζ = ds.ζ.isel(yF=slice(0, -2), zF=slice(0, -2), time=n).squeeze() + img_ζ = ζ.plot.pcolormesh(ax=ax_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, + norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + fig.colorbar(img_ζ, ax=ax_mesh, extend="both") + ax_mesh.axhline(y=ds.zF[127], color="tab:orange", alpha=0.5) + ax_mesh.set_title("vorticity") + ax_mesh.set_xlabel("y") + ax_mesh.set_ylabel("z") + ax_mesh.set_aspect("equal") + + print(f"Saving vorticity frame {n}/{ds.time.size-1}...") + plt.savefig(f"lid_driven_cavity_vorticity_Re{Re}_{n:05d}.png") plt.close("all") Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") -plot_lid_driven_cavity_vorticity(100, ds.time.size-1) +# plot_velocity_frame(Re, ds.time.size-1) +plot_vorticity_frame(Re, ds.time.size-1) From 297f90a3b3400be864f75f37e568e15d450f9152 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Mon, 9 Mar 2020 23:03:51 -0400 Subject: [PATCH 18/31] Parallelize plotting and make movies of velocity and vorticity --- .../plot_lid_driven_cavity.py | 40 ++++++++++++++++--- 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 11fab4e24c..86d4f55eb7 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -35,6 +35,8 @@ def plot_velocity_frame(Re, n): ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + Ny = ds.yC.size + Nz = ds.zC.size fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 9), dpi=200) plt.subplots_adjust(hspace=0.25) @@ -43,7 +45,7 @@ def plot_velocity_frame(Re, n): ax_v_line, ax_v_mesh = axes[0, 0], axes[0, 1] ax_w_line, ax_w_mesh = axes[1, 0], axes[1, 1] - v_line = ds.v.isel(time=n, yF=64) + v_line = ds.v.isel(time=n, yF=Ny//2) ax_v_line.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none") ax_v_line.plot(ds.zC, v_line.values.flatten(), label="Oceananigans.jl", color="tab:blue") ax_v_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) @@ -51,7 +53,7 @@ def plot_velocity_frame(Re, n): ax_v_line.set_ylabel("v") ax_v_line.set_xlim([0, 1]) - w_line = ds.w.isel(time=n, zF=64) + w_line = ds.w.isel(time=n, zF=Nz//2) ax_w_line.plot(z_Ghia, w_Ghia[Re], label="Ghia et al. (1982)", color="tab:orange", linestyle="", marker="o", fillstyle="none") ax_w_line.plot(ds.yC, w_line.values.flatten(), label="Oceananigans.jl", color="tab:orange") ax_w_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) @@ -83,13 +85,15 @@ def plot_velocity_frame(Re, n): def plot_vorticity_frame(Re, n): ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + Ny = ds.yC.size + Nz = ds.zC.size fig, axes = plt.subplots(ncols=2, figsize=(16, 9), dpi=200) fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) ax_line, ax_mesh = axes[0], axes[1] - ζ_line = -ds.ζ.isel(time=n, zF=127) + ζ_line = -ds.ζ.isel(time=n, zF=Nz-1) ax_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") ax_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:red") ax_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) @@ -102,7 +106,7 @@ def plot_vorticity_frame(Re, n): img_ζ = ζ.plot.pcolormesh(ax=ax_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) fig.colorbar(img_ζ, ax=ax_mesh, extend="both") - ax_mesh.axhline(y=ds.zF[127], color="tab:orange", alpha=0.5) + ax_mesh.axhline(y=ds.zF[Nz-1], color="tab:red", alpha=0.5) ax_mesh.set_title("vorticity") ax_mesh.set_xlabel("y") ax_mesh.set_ylabel("z") @@ -114,5 +118,29 @@ def plot_vorticity_frame(Re, n): Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") -# plot_velocity_frame(Re, ds.time.size-1) -plot_vorticity_frame(Re, ds.time.size-1) + +joblib.Parallel(n_jobs=-1)( + joblib.delayed(plot_velocity_frame)(Re, n) + for n in range(ds.time.size) +) + +joblib.Parallel(n_jobs=-1)( + joblib.delayed(plot_vorticity_frame)(Re, n) + for n in range(ds.time.size) +) + +( + ffmpeg + .input(f"lid_driven_cavity_velocity_Re{Re}_%05d.png", framerate=30) + .output(f"lid_driven_cavity_velocity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') + .overwrite_output() + .run() +) + +( + ffmpeg + .input(f"lid_driven_cavity_vorticity_Re{Re}_%05d.png", framerate=30) + .output(f"lid_driven_cavity_vorticity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') + .overwrite_output() + .run() +) From 24f368ff221cdbd9411007b1732eba5d1cc99634 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Wed, 11 Mar 2020 08:11:28 -0400 Subject: [PATCH 19/31] Modularize lid-driven cavity verification --- .../lid_driven_cavity/lid_driven_cavity.jl | 120 ++++++++++-------- 1 file changed, 67 insertions(+), 53 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index d33e2c82fb..a6876ace4c 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -5,69 +5,83 @@ using Oceananigans.Diagnostics using Oceananigans.OutputWriters using Oceananigans.AbstractOperations -Ny = Nz = 128 -Ly = Lz = 1.0 - -topology = (Flat, Bounded, Bounded) -domain = (x=(0, 1), y=(0, Ly), z=(0, Lz)) -grid = RegularCartesianGrid(topology=topology, size=(1, Ny, Nz); domain...) - -v_bcs = VVelocityBoundaryConditions(grid, - top = ValueBoundaryCondition(1.0), - bottom = ValueBoundaryCondition(0.0) -) - -w_bcs = WVelocityBoundaryConditions(grid, - north = ValueBoundaryCondition(0.0), - south = ValueBoundaryCondition(0.0) -) - -Re = 100 # Reynolds number - -model = IncompressibleModel( - grid=grid, - buoyancy=nothing, - tracers=nothing, - coriolis=nothing, - boundary_conditions = (v=v_bcs, w=w_bcs), - closure = ConstantIsotropicDiffusivity(ν=1/Re) -) - -u, v, w = model.velocities -ζ_op = ∂y(w) - ∂z(v) -ζ = Field(Cell, Face, Face, model.architecture, model.grid, TracerBoundaryConditions(grid)) -ζ_comp = Computation(ζ_op, ζ) - -max_Δt = 0.25 * model.grid.Δy^2 * Re / 2 -wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=max_Δt) -cfl = AdvectiveCFL(wizard) -dcfl = DiffusiveCFL(wizard) +function simulate_lid_driven_cavity(; Re, end_time=10.0) + Ny = Nz = 128 + Ly = Lz = 1.0 + + topology = (Flat, Bounded, Bounded) + domain = (x=(0, 1), y=(0, Ly), z=(0, Lz)) + grid = RegularCartesianGrid(topology=topology, size=(1, Ny, Nz); domain...) + + v_bcs = VVelocityBoundaryConditions(grid, + top = ValueBoundaryCondition(1.0), + bottom = ValueBoundaryCondition(0.0) + ) + + w_bcs = WVelocityBoundaryConditions(grid, + north = ValueBoundaryCondition(0.0), + south = ValueBoundaryCondition(0.0) + ) + + model = IncompressibleModel( + grid = grid, + buoyancy = nothing, + tracers = nothing, + coriolis = nothing, + boundary_conditions = (v=v_bcs, w=w_bcs), + closure = ConstantIsotropicDiffusivity(ν=1/Re) + ) + + u, v, w = model.velocities + ζ_op = ∂y(w) - ∂z(v) + ζ = Field(Cell, Face, Face, model.architecture, model.grid, TracerBoundaryConditions(grid)) + ζ_computation = Computation(ζ_op, ζ) + + fields = Dict( + "v" => model.velocities.v, + "w" => model.velocities.w, + "ζ" => model -> ζ_computation(model) + ) + + dims = Dict("ζ" => ("xC", "yF", "zF")) + global_attributes = Dict("Re" => Re) + output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) + + field_output_writer = + NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.01, + global_attributes=global_attributes, output_attributes=output_attributes, + dimensions=dims) + + max_Δt = 0.25 * model.grid.Δy^2 * Re / 2 # Make sure not to violate diffusive CFL. + wizard = TimeStepWizard(cfl=0.1, Δt=1e-6, max_change=1.1, max_Δt=max_Δt) + + cfl = AdvectiveCFL(wizard) + dcfl = DiffusiveCFL(wizard) + + simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=print_progress, + progress_frequency=20, parameters=(cfl=cfl, dcfl=dcfl)) + + simulation.output_writers[:fields] = field_output_writer + + run!(simulation) + + return simulation +end function print_progress(simulation) model = simulation.model + cfl, dcfl = simulation.parameters # Calculate simulation progress in %. progress = 100 * (model.clock.time / simulation.stop_time) # Find maximum velocities. - vmax = maximum(abs, model.velocities.v.data.parent) - wmax = maximum(abs, model.velocities.w.data.parent) + vmax = maximum(abs, interior(model.velocities.v)) + wmax = maximum(abs, interior(model.velocities.w)) i, t = model.clock.iteration, model.clock.time @printf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e\n", progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt) -end - -simulation = Simulation(model, Δt=wizard, stop_time=100, progress=print_progress, progress_frequency=20) -fields = Dict("v" => model.velocities.v, "w" => model.velocities.w, "ζ" => model -> ζ_comp(model)) -dims = Dict("ζ" => ("xC", "yF", "zF")) -global_attributes = Dict("Re" => Re) -output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) - -simulation.output_writers[:fields] = - NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.1, - global_attributes=global_attributes, output_attributes=output_attributes, - dimensions=dims) - -run!(simulation) + return nothing +end From c2a6fefe5ed2b681a6647e61ceb0f27cbc9573dd Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Wed, 11 Mar 2020 09:35:12 -0400 Subject: [PATCH 20/31] Combine velocity and vorticity plots --- .../plot_lid_driven_cavity.py | 72 ++++++------------- 1 file changed, 23 insertions(+), 49 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 86d4f55eb7..8cdba6b50d 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -33,17 +33,18 @@ 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan] } -def plot_velocity_frame(Re, n): +def plot_lid_driven_cavity_frame(Re, n): ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") Ny = ds.yC.size Nz = ds.zC.size - fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16, 9), dpi=200) + fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(16, 16), dpi=200) plt.subplots_adjust(hspace=0.25) fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) ax_v_line, ax_v_mesh = axes[0, 0], axes[0, 1] ax_w_line, ax_w_mesh = axes[1, 0], axes[1, 1] + ax_ζ_line, ax_ζ_mesh = axes[2, 0], axes[2, 1] v_line = ds.v.isel(time=n, yF=Ny//2) ax_v_line.plot(y_Ghia, v_Ghia[Re], label="Ghia et al. (1982)", color="tab:blue", linestyle="", marker="o", fillstyle="none") @@ -79,68 +80,41 @@ def plot_velocity_frame(Re, n): ax_w_mesh.set_ylabel("z") ax_w_mesh.set_aspect("equal") - print(f"Saving velocity frame {n}/{ds.time.size-1}...") - plt.savefig(f"lid_driven_cavity_velocity_Re{Re}_{n:05d}.png") - plt.close("all") - -def plot_vorticity_frame(Re, n): - ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") - Ny = ds.yC.size - Nz = ds.zC.size - - fig, axes = plt.subplots(ncols=2, figsize=(16, 9), dpi=200) - fig.suptitle(f"Lid-driven cavity, Re = {Re}, t = {ds.time[n].values:.2f}", fontsize=16) - - ax_line, ax_mesh = axes[0], axes[1] - ζ_line = -ds.ζ.isel(time=n, zF=Nz-1) - ax_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") - ax_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:red") - ax_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) - ax_line.set_xlabel("y") - ax_line.set_ylabel("vorticity $\zeta$") - ax_line.set_xlim([0, 1]) - ax_line.set_ylim(bottom=0) + ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") + ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:red") + ax_ζ_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) + ax_ζ_line.set_xlabel("y") + ax_ζ_line.set_ylabel("vorticity $\zeta$") + ax_ζ_line.set_xlim([0, 1]) + ax_ζ_line.set_ylim(bottom=0) ζ = ds.ζ.isel(yF=slice(0, -2), zF=slice(0, -2), time=n).squeeze() - img_ζ = ζ.plot.pcolormesh(ax=ax_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, + img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) - fig.colorbar(img_ζ, ax=ax_mesh, extend="both") - ax_mesh.axhline(y=ds.zF[Nz-1], color="tab:red", alpha=0.5) - ax_mesh.set_title("vorticity") - ax_mesh.set_xlabel("y") - ax_mesh.set_ylabel("z") - ax_mesh.set_aspect("equal") - - print(f"Saving vorticity frame {n}/{ds.time.size-1}...") - plt.savefig(f"lid_driven_cavity_vorticity_Re{Re}_{n:05d}.png") + fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both") + ax_ζ_mesh.axhline(y=ds.zF[Nz-1], color="tab:red", alpha=0.5) + ax_ζ_mesh.set_title("vorticity") + ax_ζ_mesh.set_xlabel("y") + ax_ζ_mesh.set_ylabel("z") + ax_ζ_mesh.set_aspect("equal") + + print(f"Saving lid-driven cavity Re={Re} frame {n}/{ds.time.size-1}...") + plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") joblib.Parallel(n_jobs=-1)( - joblib.delayed(plot_velocity_frame)(Re, n) + joblib.delayed(plot_lid_driven_cavity_frame)(Re, n) for n in range(ds.time.size) ) -joblib.Parallel(n_jobs=-1)( - joblib.delayed(plot_vorticity_frame)(Re, n) - for n in range(ds.time.size) -) - -( - ffmpeg - .input(f"lid_driven_cavity_velocity_Re{Re}_%05d.png", framerate=30) - .output(f"lid_driven_cavity_velocity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') - .overwrite_output() - .run() -) - ( ffmpeg - .input(f"lid_driven_cavity_vorticity_Re{Re}_%05d.png", framerate=30) - .output(f"lid_driven_cavity_vorticity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') + .input(f"lid_driven_cavity_Re{Re}_%05d.png", framerate=30) + .output(f"lid_driven_cavity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') .overwrite_output() .run() ) From 2e6ae8517e2f4462ac9c07c0d30ae7e552f76046 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Thu, 23 Jul 2020 07:51:00 -0400 Subject: [PATCH 21/31] Update output plots a bit --- verification/lid_driven_cavity/lid_driven_cavity.jl | 2 +- verification/lid_driven_cavity/plot_lid_driven_cavity.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index a6876ace4c..7953385396 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -48,7 +48,7 @@ function simulate_lid_driven_cavity(; Re, end_time=10.0) output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) field_output_writer = - NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.01, + NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.1, global_attributes=global_attributes, output_attributes=output_attributes, dimensions=dims) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 8cdba6b50d..d5acd5737b 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -91,7 +91,7 @@ def plot_lid_driven_cavity_frame(Re, n): ζ = ds.ζ.isel(yF=slice(0, -2), zF=slice(0, -2), time=n).squeeze() img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, - norm=colors.SymLogNorm(linthresh=1e-2, vmin=-1e2, vmax=1e2)) + norm=colors.SymLogNorm(base=10, linthresh=1e-2, vmin=-1e2, vmax=1e2)) fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both") ax_ζ_mesh.axhline(y=ds.zF[Nz-1], color="tab:red", alpha=0.5) ax_ζ_mesh.set_title("vorticity") From a8cda4716d5ac655bba60e4c0336ccec270b76a4 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Thu, 23 Jul 2020 14:14:37 -0400 Subject: [PATCH 22/31] Lid-driven cavity at Re=400 --- verification/lid_driven_cavity/lid_driven_cavity.jl | 4 ++++ verification/lid_driven_cavity/plot_lid_driven_cavity.py | 9 ++++++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 7953385396..fbda6a2c9f 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -85,3 +85,7 @@ function print_progress(simulation) return nothing end + +# simulate_lid_driven_cavity(Re=100, end_time=25.0) +# simulate_lid_driven_cavity(Re=400, end_time=50.0) + diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index d5acd5737b..9550658d4b 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -24,13 +24,15 @@ z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000] w_Ghia = { - 100: [0.0000, 0.09233, 0.10091, 0.1089, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000] + 100: [0.0000, 0.09233, 0.10091, 0.10890, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000], + 400: [0.0000, 0.18360, 0.19713, 0.20920, 0.22965, 0.28124, 0.30203, 0.30174, 0.05186, -0.38598, -0.44993, -0.23827, -0.22847, -0.19254, -0.15663, -0.12146, 0.0000] } y_ζ_Ghia = [0.0000, 0.0625, 0.1250, 0.1875, 0.2500, 0.3125, 0.3750, 0.4375, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8125, 0.8750, 0.9375, 1.0000] ζ_Ghia = { - 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan] + 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan], + 400: [nan, 53.6863, 34.6351, 26.5825, 21.0985, 16.8900, 13.7040, 11.4537, 10.0545, 9.38889, 9.34599, 9.88979, 11.2018, 13.9068, 19.6859, 35.0773, nan] } def plot_lid_driven_cavity_frame(Re, n): @@ -103,7 +105,7 @@ def plot_lid_driven_cavity_frame(Re, n): plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -Re = 100 +Re = 400 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") joblib.Parallel(n_jobs=-1)( @@ -118,3 +120,4 @@ def plot_lid_driven_cavity_frame(Re, n): .overwrite_output() .run() ) + From 5250349dcec16622de6e851475f0845046a1ea08 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Thu, 23 Jul 2020 21:54:14 -0400 Subject: [PATCH 23/31] Tedious data entry --- .../plot_lid_driven_cavity.py | 35 +++++++++++++------ 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 9550658d4b..ce77059ddc 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -9,30 +9,45 @@ from numpy import nan ##### -##### Data from tables 1 and 2 of Ghia et al. (1982). +##### Data from tables 1, 2, and 4 of Ghia et al. (1982). ##### j_Ghia = [1, 8, 9, 10, 14, 23, 37, 59, 65, 80, 95, 110, 123, 124, 125, 126, 129] y_Ghia = [0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531, 0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000] v_Ghia = { - 100: [0.0000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.0000], - 400: [0.0000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.0000] + 100: [0.00000, -0.03717, -0.04192, -0.04775, -0.06434, -0.10150, -0.15662, -0.21090, -0.20581, -0.13641, 0.00332, 0.23151, 0.68717, 0.73722, 0.78871, 0.84123, 1.00000], + 400: [0.00000, -0.08186, -0.09266, -0.10338, -0.14612, -0.24299, -0.32726, -0.17119, -0.11477, 0.02135, 0.16256, 0.29093, 0.55892, 0.61756, 0.68439, 0.75837, 1.00000], + 1000: [0.00000, -0.18109, -0.20196, -0.22220, -0.29730, -0.38289, -0.27805, -0.10648, -0.06080, 0.05702, 0.18719, 0.33304, 0.46604, 0.51117, 0.57492, 0.65928, 1.00000], + 3200: [0.00000, -0.32407, -0.35344, -0.37827, -0.41933, -0.34323, -0.24427, -0.86636, -0.04272, 0.07156, 0.19791, 0.34682, 0.46101, 0.46547, 0.48296, 0.53236, 1.00000], + 5000: [0.00000, -0.41165, -0.42901, -0.43643, -0.40435, -0.33050, -0.22855, -0.07404, -0.03039, 0.08183, 0.20087, 0.33556, 0.46036, 0.45992, 0.46120, 0.48223, 1.00000], + 7500: [0.00000, -0.43154, -0.43590, -0.43025, -0.38324, -0.32393, -0.23176, -0.07503, -0.03800, 0.08342, 0.20591, 0.34228, 0.47167, 0.47323, 0.47048, 0.47244, 1.00000], + 10000: [0.00000, -0.42735, -0.42537, -0.41657, -0.38000, -0.32709, -0.23186, -0.07540, 0.03111, 0.08344, 0.20673, 0.34635, 0.47804, 0.48070, 0.47783, 0.47221, 1.00000] } k_Ghia = [1, 9, 10, 11, 13, 21, 30, 31, 65, 104, 111, 117, 122, 123, 124, 125, 129] z_Ghia = [0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344, 0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000] w_Ghia = { - 100: [0.0000, 0.09233, 0.10091, 0.10890, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.0000], - 400: [0.0000, 0.18360, 0.19713, 0.20920, 0.22965, 0.28124, 0.30203, 0.30174, 0.05186, -0.38598, -0.44993, -0.23827, -0.22847, -0.19254, -0.15663, -0.12146, 0.0000] + 100: [0.00000, 0.09233, 0.10091, 0.10890, 0.12317, 0.16077, 0.17507, 0.17527, 0.05454, -0.24533, -0.22445, -0.16914, -0.10313, -0.08864, -0.07391, -0.05906, 0.00000], + 400: [0.00000, 0.18360, 0.19713, 0.20920, 0.22965, 0.28124, 0.30203, 0.30174, 0.05186, -0.38598, -0.44993, -0.23827, -0.22847, -0.19254, -0.15663, -0.12146, 0.00000], + 1000: [0.00000, 0.27485, 0.29012, 0.30353, 0.32627, 0.37095, 0.33075, 0.32235, 0.02526, -0.31966, -0.42665, -0.51550, -0.39188, -0.33714, -0.27669, -0.21388, 0.00000], + 3200: [0.00000, 0.39560, 0.40917, 0.41906, 0.42768, 0.37119, 0.29030, 0.28188, 0.00999, -0.31184, -0.37401, -0.44307, -0.54053, -0.52357, -0.47425, -0.39017, 0.00000], + 5000: [0.00000, 0.42447, 0.43329, 0.43648, 0.42951, 0.35368, 0.28066, 0.27280, 0.00945, -0.30018, -0.36214, -0.41442, -0.52876, -0.55408, -0.55069, -0.49774, 0.00000], + 7500: [0.00000, 0.43979, 0.44030, 0.43564, 0.41824, 0.35060, 0.28117, 0.27348, 0.00824, -0.30448, -0.36213, -0.41050, -0.48590, -0.52347, -0.55216, -0.53858, 0.00000], + 10000: [0.00000, 0.43983, 0.43733, 0.43124, 0.41487, 0.35070, 0.28003, 0.27224, 0.00831, -0.30719, -0.36737, -0.41496, -0.45863, -0.49099, -0.52987, -0.54302, 0.00000] } y_ζ_Ghia = [0.0000, 0.0625, 0.1250, 0.1875, 0.2500, 0.3125, 0.3750, 0.4375, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8125, 0.8750, 0.9375, 1.0000] ζ_Ghia = { - 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan], - 400: [nan, 53.6863, 34.6351, 26.5825, 21.0985, 16.8900, 13.7040, 11.4537, 10.0545, 9.38889, 9.34599, 9.88979, 11.2018, 13.9068, 19.6859, 35.0773, nan] + 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan], + 400: [nan, 53.6863, 34.6351, 26.5825, 21.0985, 16.8900, 13.7040, 11.4537, 10.0545, 9.38889, 9.34599, 9.88979, 11.2018, 13.9068, 19.6859, 35.0773, nan], + 1000: [nan, 42.1124, 23.8707, 18.3120, 16.0458, 14.8061, 14.1374, 14.0928, 14.8901, 16.8350, 20.2666, 25.4341, 32.2953, 40.5437, 51.0557, 75.5980, nan], + 3200: [nan, 49.9664, 34.2327, 30.4779, 27.9514, 25.8572, 24.4639, 24.1457, 25.3889, 28.9413, 35.8795, 47.1443, 61.7864, 75.6401, 89.3391, 126.670, nan], + 5000: [nan, 56.7091, 41.3394, 38.0436, 35.3504, 33.0486, 31.5791, 31.3793, 33.0115, 37.3609, 45.8622, 60.0065, 77.9509, 91.5682, 103.436, 146.702, nan], + 7500: [nan, 61.4046, 50.0769, 46.8901, 43.5641, 40.6123, 38.6951, 38.3834, 40.3982, 45.9128, 56.9345, 75.6334, 98.2364, 111.115, 125.131, 180.927, nan], + 10000: [nan, 66.0352, 57.7756, 54.3725, 50.3792, 46.8672, 44.6303, 44.3287, 46.8271, 53.5905, 67.1400, 90.0231, 116.275, 127.928, 145.073, 209.452, nan] } def plot_lid_driven_cavity_frame(Re, n): @@ -82,9 +97,9 @@ def plot_lid_driven_cavity_frame(Re, n): ax_w_mesh.set_ylabel("z") ax_w_mesh.set_aspect("equal") - ζ_line = -ds.ζ.isel(time=n, zF=Nz-1) + ζ_line = ds.ζ.isel(time=n, zF=Nz-1) ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") - ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:red") + ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:purple") ax_ζ_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) ax_ζ_line.set_xlabel("y") ax_ζ_line.set_ylabel("vorticity $\zeta$") @@ -95,7 +110,7 @@ def plot_lid_driven_cavity_frame(Re, n): img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, norm=colors.SymLogNorm(base=10, linthresh=1e-2, vmin=-1e2, vmax=1e2)) fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both") - ax_ζ_mesh.axhline(y=ds.zF[Nz-1], color="tab:red", alpha=0.5) + ax_ζ_mesh.axhline(y=ds.zF[Nz-2], color="tab:purple", alpha=0.5) ax_ζ_mesh.set_title("vorticity") ax_ζ_mesh.set_xlabel("y") ax_ζ_mesh.set_ylabel("z") From ed4f4d5435849297b562fa4e34d535c9630a076d Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Thu, 23 Jul 2020 21:57:33 -0400 Subject: [PATCH 24/31] Run all cases --- verification/lid_driven_cavity/lid_driven_cavity.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index fbda6a2c9f..40536d3c64 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -86,6 +86,11 @@ function print_progress(simulation) return nothing end -# simulate_lid_driven_cavity(Re=100, end_time=25.0) -# simulate_lid_driven_cavity(Re=400, end_time=50.0) + simulate_lid_driven_cavity(Re=100, end_time=15) + simulate_lid_driven_cavity(Re=400, end_time=20) + simulate_lid_driven_cavity(Re=1000, end_time=25) + simulate_lid_driven_cavity(Re=3200, end_time=50) + simulate_lid_driven_cavity(Re=5000, end_time=50) + simulate_lid_driven_cavity(Re=7500, end_time=75) + simulate_lid_driven_cavity(Re=10000, end_time=100) From d8967352045141757682e91c00427dcb8fe00006 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sun, 26 Jul 2020 09:05:34 -0400 Subject: [PATCH 25/31] Better vorticity plots --- .../lid_driven_cavity/plot_lid_driven_cavity.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index ce77059ddc..378b176784 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -7,6 +7,7 @@ import ffmpeg from numpy import nan +from xarray.ufuncs import fabs ##### ##### Data from tables 1, 2, and 4 of Ghia et al. (1982). @@ -97,20 +98,21 @@ def plot_lid_driven_cavity_frame(Re, n): ax_w_mesh.set_ylabel("z") ax_w_mesh.set_aspect("equal") - ζ_line = ds.ζ.isel(time=n, zF=Nz-1) - ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:red", linestyle="", marker="o", fillstyle="none") + ζ_line = fabs(ds.ζ.isel(time=n, zF=Nz-1)) + ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:purple", linestyle="", marker="o", fillstyle="none") ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:purple") ax_ζ_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) ax_ζ_line.set_xlabel("y") - ax_ζ_line.set_ylabel("vorticity $\zeta$") + ax_ζ_line.set_ylabel("vorticity $|\zeta$|") ax_ζ_line.set_xlim([0, 1]) ax_ζ_line.set_ylim(bottom=0) ζ = ds.ζ.isel(yF=slice(0, -2), zF=slice(0, -2), time=n).squeeze() img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, - norm=colors.SymLogNorm(base=10, linthresh=1e-2, vmin=-1e2, vmax=1e2)) + # vmin=-100, vmax=100) + norm=colors.SymLogNorm(base=10, linthresh=1, vmin=-1e2, vmax=1e2)) fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both") - ax_ζ_mesh.axhline(y=ds.zF[Nz-2], color="tab:purple", alpha=0.5) + ax_ζ_mesh.axhline(y=ds.zF[Nz-2], color="tab:purple", alpha=1.0) ax_ζ_mesh.set_title("vorticity") ax_ζ_mesh.set_xlabel("y") ax_ζ_mesh.set_ylabel("z") @@ -120,7 +122,7 @@ def plot_lid_driven_cavity_frame(Re, n): plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -Re = 400 +Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") joblib.Parallel(n_jobs=-1)( From 539cda3a0f44310006973a46345d6800184fda3f Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sun, 26 Jul 2020 11:08:27 -0400 Subject: [PATCH 26/31] Manual data entry gone wrong --- .../lid_driven_cavity/plot_lid_driven_cavity.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index 378b176784..fe5be26c2c 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -44,11 +44,11 @@ ζ_Ghia = { 100: [nan, 40.0110, 22.5378, 16.2862, 12.7844, 10.4199, 8.69628, 7.43218, 6.57451, 6.13973, 6.18946, 6.82674, 8.22110, 10.7414, 15.6591, 30.7923, nan], 400: [nan, 53.6863, 34.6351, 26.5825, 21.0985, 16.8900, 13.7040, 11.4537, 10.0545, 9.38889, 9.34599, 9.88979, 11.2018, 13.9068, 19.6859, 35.0773, nan], - 1000: [nan, 42.1124, 23.8707, 18.3120, 16.0458, 14.8061, 14.1374, 14.0928, 14.8901, 16.8350, 20.2666, 25.4341, 32.2953, 40.5437, 51.0557, 75.5980, nan], - 3200: [nan, 49.9664, 34.2327, 30.4779, 27.9514, 25.8572, 24.4639, 24.1457, 25.3889, 28.9413, 35.8795, 47.1443, 61.7864, 75.6401, 89.3391, 126.670, nan], - 5000: [nan, 56.7091, 41.3394, 38.0436, 35.3504, 33.0486, 31.5791, 31.3793, 33.0115, 37.3609, 45.8622, 60.0065, 77.9509, 91.5682, 103.436, 146.702, nan], - 7500: [nan, 61.4046, 50.0769, 46.8901, 43.5641, 40.6123, 38.6951, 38.3834, 40.3982, 45.9128, 56.9345, 75.6334, 98.2364, 111.115, 125.131, 180.927, nan], - 10000: [nan, 66.0352, 57.7756, 54.3725, 50.3792, 46.8672, 44.6303, 44.3287, 46.8271, 53.5905, 67.1400, 90.0231, 116.275, 127.928, 145.073, 209.452, nan] + 1000: [nan, 42.1124, 23.8707, 18.3120, 16.0458, 14.8061, 14.1374, 14.0928, 14.8901, 16.8350, 20.2666, 25.4341, 32.2953, 40.5437, 51.0557, 75.5980, nan][::-1], + 3200: [nan, 49.9664, 34.2327, 30.4779, 27.9514, 25.8572, 24.4639, 24.1457, 25.3889, 28.9413, 35.8795, 47.1443, 61.7864, 75.6401, 89.3391, 126.670, nan][::-1], + 5000: [nan, 56.7091, 41.3394, 38.0436, 35.3504, 33.0486, 31.5791, 31.3793, 33.0115, 37.3609, 45.8622, 60.0065, 77.9509, 91.5682, 103.436, 146.702, nan][::-1], + 7500: [nan, 61.4046, 50.0769, 46.8901, 43.5641, 40.6123, 38.6951, 38.3834, 40.3982, 45.9128, 56.9345, 75.6334, 98.2364, 111.115, 125.131, 180.927, nan][::-1], + 10000: [nan, 66.0352, 57.7756, 54.3725, 50.3792, 46.8672, 44.6303, 44.3287, 46.8271, 53.5905, 67.1400, 90.0231, 116.275, 127.928, 145.073, 209.452, nan][::-1] } def plot_lid_driven_cavity_frame(Re, n): @@ -122,7 +122,7 @@ def plot_lid_driven_cavity_frame(Re, n): plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -Re = 100 +Re = 10000 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") joblib.Parallel(n_jobs=-1)( From d060385aa7f4c498819a133d5ade7c546177fcac Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 11 Sep 2020 22:11:40 -0400 Subject: [PATCH 27/31] Update `lid_driven_cavity.jl` --- verification/lid_driven_cavity/lid_driven_cavity.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 40536d3c64..ec0efa6e68 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -5,12 +5,11 @@ using Oceananigans.Diagnostics using Oceananigans.OutputWriters using Oceananigans.AbstractOperations -function simulate_lid_driven_cavity(; Re, end_time=10.0) +function simulate_lid_driven_cavity(; Re, end_time) Ny = Nz = 128 - Ly = Lz = 1.0 topology = (Flat, Bounded, Bounded) - domain = (x=(0, 1), y=(0, Ly), z=(0, Lz)) + domain = (x=(0, 1), y=(0, 1), z=(0, 1)) grid = RegularCartesianGrid(topology=topology, size=(1, Ny, Nz); domain...) v_bcs = VVelocityBoundaryConditions(grid, @@ -29,7 +28,7 @@ function simulate_lid_driven_cavity(; Re, end_time=10.0) tracers = nothing, coriolis = nothing, boundary_conditions = (v=v_bcs, w=w_bcs), - closure = ConstantIsotropicDiffusivity(ν=1/Re) + closure = IsotropicDiffusivity(ν=1/Re) ) u, v, w = model.velocities @@ -48,7 +47,7 @@ function simulate_lid_driven_cavity(; Re, end_time=10.0) output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) field_output_writer = - NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", interval=0.1, + NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", time_interval=0.1, global_attributes=global_attributes, output_attributes=output_attributes, dimensions=dims) @@ -59,7 +58,7 @@ function simulate_lid_driven_cavity(; Re, end_time=10.0) dcfl = DiffusiveCFL(wizard) simulation = Simulation(model, Δt=wizard, stop_time=end_time, progress=print_progress, - progress_frequency=20, parameters=(cfl=cfl, dcfl=dcfl)) + iteration_interval=20, parameters=(cfl=cfl, dcfl=dcfl)) simulation.output_writers[:fields] = field_output_writer From 05c1f7741b4c41d5095ed24de8840349ef16d6e6 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 11 Sep 2020 22:30:34 -0400 Subject: [PATCH 28/31] Update vorticity plots --- verification/lid_driven_cavity/plot_lid_driven_cavity.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index fe5be26c2c..ea6fd8db61 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -98,7 +98,7 @@ def plot_lid_driven_cavity_frame(Re, n): ax_w_mesh.set_ylabel("z") ax_w_mesh.set_aspect("equal") - ζ_line = fabs(ds.ζ.isel(time=n, zF=Nz-1)) + ζ_line = fabs(ds.ζ.isel(time=n, zF=Nz)) ax_ζ_line.plot(y_ζ_Ghia, ζ_Ghia[Re], label="Ghia et al. (1982)", color="tab:purple", linestyle="", marker="o", fillstyle="none") ax_ζ_line.plot(ds.yF, ζ_line.values.flatten(), label="Oceananigans.jl", color="tab:purple") ax_ζ_line.legend(loc="lower left", bbox_to_anchor=(0, 1.01, 1, 0.2), ncol=2, frameon=False) @@ -107,12 +107,11 @@ def plot_lid_driven_cavity_frame(Re, n): ax_ζ_line.set_xlim([0, 1]) ax_ζ_line.set_ylim(bottom=0) - ζ = ds.ζ.isel(yF=slice(0, -2), zF=slice(0, -2), time=n).squeeze() + ζ = ds.ζ.isel(time=n).squeeze() img_ζ = ζ.plot.pcolormesh(ax=ax_ζ_mesh, cmap=cmocean.cm.curl, extend="both", add_colorbar=False, - # vmin=-100, vmax=100) norm=colors.SymLogNorm(base=10, linthresh=1, vmin=-1e2, vmax=1e2)) fig.colorbar(img_ζ, ax=ax_ζ_mesh, extend="both") - ax_ζ_mesh.axhline(y=ds.zF[Nz-2], color="tab:purple", alpha=1.0) + ax_ζ_mesh.axhline(y=1, color="tab:purple", alpha=1.0) ax_ζ_mesh.set_title("vorticity") ax_ζ_mesh.set_xlabel("y") ax_ζ_mesh.set_ylabel("z") @@ -122,7 +121,7 @@ def plot_lid_driven_cavity_frame(Re, n): plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -Re = 10000 +Re = 100 ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") joblib.Parallel(n_jobs=-1)( From eb4002b5c81846dc7827633880a5bc917b9cf391 Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Fri, 11 Sep 2020 22:34:31 -0400 Subject: [PATCH 29/31] Automate everything --- .../lid_driven_cavity/lid_driven_cavity.jl | 20 ++++++------- .../plot_lid_driven_cavity.py | 29 +++++++++---------- 2 files changed, 23 insertions(+), 26 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index ec0efa6e68..beda690e08 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -5,12 +5,10 @@ using Oceananigans.Diagnostics using Oceananigans.OutputWriters using Oceananigans.AbstractOperations -function simulate_lid_driven_cavity(; Re, end_time) - Ny = Nz = 128 - +function simulate_lid_driven_cavity(; Re, N, end_time) topology = (Flat, Bounded, Bounded) domain = (x=(0, 1), y=(0, 1), z=(0, 1)) - grid = RegularCartesianGrid(topology=topology, size=(1, Ny, Nz); domain...) + grid = RegularCartesianGrid(topology=topology, size=(1, N, N); domain...) v_bcs = VVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(1.0), @@ -85,11 +83,11 @@ function print_progress(simulation) return nothing end - simulate_lid_driven_cavity(Re=100, end_time=15) - simulate_lid_driven_cavity(Re=400, end_time=20) - simulate_lid_driven_cavity(Re=1000, end_time=25) - simulate_lid_driven_cavity(Re=3200, end_time=50) - simulate_lid_driven_cavity(Re=5000, end_time=50) - simulate_lid_driven_cavity(Re=7500, end_time=75) - simulate_lid_driven_cavity(Re=10000, end_time=100) + simulate_lid_driven_cavity(Re=100, N=128, end_time=15) + simulate_lid_driven_cavity(Re=400, N=128, end_time=20) + simulate_lid_driven_cavity(Re=1000, N=128, end_time=25) + simulate_lid_driven_cavity(Re=3200, N=128, end_time=50) + simulate_lid_driven_cavity(Re=5000, N=256, end_time=50) + simulate_lid_driven_cavity(Re=7500, N=256, end_time=75) + simulate_lid_driven_cavity(Re=10000, N=256, end_time=100) diff --git a/verification/lid_driven_cavity/plot_lid_driven_cavity.py b/verification/lid_driven_cavity/plot_lid_driven_cavity.py index ea6fd8db61..4256b25d69 100644 --- a/verification/lid_driven_cavity/plot_lid_driven_cavity.py +++ b/verification/lid_driven_cavity/plot_lid_driven_cavity.py @@ -121,19 +121,18 @@ def plot_lid_driven_cavity_frame(Re, n): plt.savefig(f"lid_driven_cavity_Re{Re}_{n:05d}.png") plt.close("all") -Re = 100 -ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") - -joblib.Parallel(n_jobs=-1)( - joblib.delayed(plot_lid_driven_cavity_frame)(Re, n) - for n in range(ds.time.size) -) - -( - ffmpeg - .input(f"lid_driven_cavity_Re{Re}_%05d.png", framerate=30) - .output(f"lid_driven_cavity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') - .overwrite_output() - .run() -) +for Re in [100, 400, 1000, 3200, 5000, 7500, 10000]: + ds = xr.open_dataset(f"lid_driven_cavity_Re{Re}.nc") + joblib.Parallel(n_jobs=-1)( + joblib.delayed(plot_lid_driven_cavity_frame)(Re, n) + for n in range(ds.time.size) + ) + + ( + ffmpeg + .input(f"lid_driven_cavity_Re{Re}_%05d.png", framerate=30) + .output(f"lid_driven_cavity_Re{Re}.mp4", crf=15, pix_fmt='yuv420p') + .overwrite_output() + .run() + ) From 038d803e610d3141f38d00d40f247475d77207da Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sat, 12 Sep 2020 08:47:58 -0400 Subject: [PATCH 30/31] Better logging --- verification/lid_driven_cavity/lid_driven_cavity.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index beda690e08..2a1654a900 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -1,10 +1,14 @@ using Printf +using Logging + using Oceananigans using Oceananigans: Face, Cell using Oceananigans.Diagnostics using Oceananigans.OutputWriters using Oceananigans.AbstractOperations +Logging.global_logger(OceananigansLogger()) + function simulate_lid_driven_cavity(; Re, N, end_time) topology = (Flat, Bounded, Bounded) domain = (x=(0, 1), y=(0, 1), z=(0, 1)) @@ -77,8 +81,8 @@ function print_progress(simulation) wmax = maximum(abs, interior(model.velocities.w)) i, t = model.clock.iteration, model.clock.time - @printf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e\n", - progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt) + @info @sprintf("[%06.2f%%] i: %d, t: %.3f, U_max: (%.2e, %.2e), CFL: %.2e, dCFL: %.2e, next Δt: %.2e", + progress, i, t, vmax, wmax, cfl(model), dcfl(model), simulation.Δt.Δt) return nothing end From 08d5cb5ce844ede1266e431f0992ba4ff611b79e Mon Sep 17 00:00:00 2001 From: ali-ramadhan Date: Sun, 20 Sep 2020 11:17:34 -0400 Subject: [PATCH 31/31] Embed YouTube videos for lid-driven cavity --- docs/make.jl | 6 +- docs/src/verification/lid_driven_cavity.md | 90 +++++++++++++++++++- docs/src/verification/taylor_green_vortex.md | 13 --- 3 files changed, 92 insertions(+), 17 deletions(-) delete mode 100644 docs/src/verification/taylor_green_vortex.md diff --git a/docs/make.jl b/docs/make.jl index e8adfd0958..2b1b627f51 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -99,9 +99,9 @@ numerical_pages = [ ] verification_pages = [ - "Taylor-Green vortex" => "verification/taylor_green_vortex.md", - "Stratified Couette flow" => "verification/stratified_couette_flow.md", - "Convergence tests" => "verification/convergence_tests.md" + "Convergence tests" => "verification/convergence_tests.md", + "Lid-driven cavity" => "verification/lid_driven_cavity.md", + "Stratified Couette flow" => "verification/stratified_couette_flow.md" ] appendix_pages = [ diff --git a/docs/src/verification/lid_driven_cavity.md b/docs/src/verification/lid_driven_cavity.md index 338f1efab0..d0c3d51e00 100644 --- a/docs/src/verification/lid_driven_cavity.md +++ b/docs/src/verification/lid_driven_cavity.md @@ -2,8 +2,96 @@ The lid-driven cavity test problem has been used for a long time as a simple verification test for computational fluid dynamics codes. First described by [Burggraf66](@cite), the fluid is contained in a square cavity with Dirchlet boundary -conditions on all four sides. The top wall moves with some velocity ``U`` while the other three walls are stationary. +conditions on all four sides. The top wall moves with velocity ``U = 1`` while the other three walls are stationary. The solution reaches a laminar steady-state whose properties can be compared with a huge amount of existing data. The canonical database is given by [Ghia82](@cite) who report detailed information on the velocity fields as well as the streamline and vorticity contours at various Reynolds numbers. More accurate data is reported by [Botella98](@cite), [Erturk05](@cite), and [Bruneau06](@cite). + +Below are lid-driven cavity simulation results from Oceananigans.jl compared with [Ghia82](@cite) for Re = 100, 400, +1000, 3200, 5000, 7500, and 10000. + +## Re = 100 + +```@raw html + +
+
+ +
+
+``` + +## Re = 400 + +```@raw html +
+
+ +
+
+``` + +## Re = 1,000 + +```@raw html +
+
+ +
+
+``` + +## Re = 3,200 + +```@raw html +
+
+ +
+
+``` + +## Re = 5,000 + +```@raw html +
+
+ +
+
+``` + +## Re = 7,500 + +```@raw html +
+
+ +
+
+``` + +## Re = 10,000 + +```@raw html +
+
+ +
+
+``` diff --git a/docs/src/verification/taylor_green_vortex.md b/docs/src/verification/taylor_green_vortex.md deleted file mode 100644 index 79e32ce2f3..0000000000 --- a/docs/src/verification/taylor_green_vortex.md +++ /dev/null @@ -1,13 +0,0 @@ -# Taylor-Green vortex -An exact solution to the two-dimensional incompressible Navier-Stokes equations is given by [Taylor37](@cite) describing -the unsteady flow of a vortex decaying under viscous dissipation. The viscous terms balance the time derivatives while -the nonlinear advection terms balance the pressure gradient term. We use the doubly-periodic solution described by -(p. 310) [Hesthaven07](@cite) - -```math -\begin{aligned} - u(x, y, t) &= -\sin(2\pi y) e^{-4\pi^2\nu t} \\ - v(x, y, t) &= \sin(2\pi x) e^{-4\pi^2\nu t} \\ - p(x, y, t) &= -\cos(2\pi x) \cos(2\pi y) e^{-8\pi^2\nu t} -\end{aligned} -```