Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Catch up with Oceananigans.jl #78

Merged
merged 10 commits into from
Oct 9, 2020
563 changes: 380 additions & 183 deletions Manifest.toml

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion benchmarks/benchmark_static_atmosphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ using JULES
const timer = TimerOutput()
const Δt = 1

grid = RegularCartesianGrid(size=(32, 32, 32), length=(1, 1, 1), halo=(2, 2, 2))
grid = RegularCartesianGrid(size=(32, 32, 32), extent=(1, 1, 1), halo=(2, 2, 2))
model = CompressibleModel(grid=grid, thermodynamic_variable=Energy())
time_step!(model, Δt) # warmup to compile

Expand Down
2 changes: 1 addition & 1 deletion src/JULES.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using Oceananigans
using Oceananigans: AbstractGrid

using Oceananigans.TurbulenceClosures:
ConstantIsotropicDiffusivity, DiffusivityFields, with_tracers
IsotropicDiffusivity, DiffusivityFields, with_tracers

export
Entropy, Energy,
Expand Down
1 change: 1 addition & 0 deletions src/Operators/Operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ using Oceananigans.Grids
using Oceananigans.Operators

using Oceananigans: AbstractGrid
using Oceananigans.TurbulenceClosures: IsotropicDiffusivity

export
kinetic_energy,
Expand Down
20 changes: 10 additions & 10 deletions src/Operators/compressible_operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -102,15 +102,15 @@
* ∂zᵃᵃᶠ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))
end

@inline function ∂ⱼDᶜⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tracer_index, ρ, ρc, args...)
@inline function ∂ⱼDᶜⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tracer_index, ρ, ρc, args...)
@inbounds κ = closure.κ[tracer_index]
return (1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, κ, ρ, ρc)
+ δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, κ, ρ, ρc)
+ δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, κ, ρ, ρc)))
end

@inline function ∂ⱼtᶜDᶜⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tvar_diag,
@inline function ∂ⱼtᶜDᶜⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tvar_diag,
tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc, args...)

@inbounds κ = closure.κ[tracer_index]
Expand All @@ -120,7 +120,7 @@ end
+ δzᵃᵃᶜ(i, j, k, grid, diffusive_tvar_flux_z, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)))
end

@inline function ∂ⱼDᵖⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tracer_index,
@inline function ∂ⱼDᵖⱼ(i, j, k, grid, closure::IsotropicDiffusivity, tracer_index,
p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃, args ...)

@inbounds κ = closure.κ[tracer_index]
Expand Down Expand Up @@ -189,34 +189,34 @@ end
@inline viscous_flux_wx(i, j, k, grid, νᶠᶜᶠ, ρ, ρw, ρu) = ℑxzᶠᵃᶠ(i, j, k, grid, ρ) * νᶠᶜᶠ * Axᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wx(i, j, k, grid, ρ, ρw, ρu)
@inline viscous_flux_wy(i, j, k, grid, νᶜᶠᶠ, ρ, ρw, ρv) = ℑyzᵃᶠᶠ(i, j, k, grid, ρ) * νᶜᶠᶠ * Ayᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wy(i, j, k, grid, ρ, ρw, ρv)

@inline ∂ⱼτ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline ∂ⱼτ₁ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, closure.ν, ρ, ρũ.ρu)
+ δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, closure.ν, ρ, ρũ.ρu, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, closure.ν, ρ, ρũ.ρu, ρũ.ρw)))

@inline ∂ⱼτ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline ∂ⱼτ₂ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, closure.ν, ρ, ρũ.ρv, ρũ.ρu)
+ δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, closure.ν, ρ, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, closure.ν, ρ, ρũ.ρv, ρũ.ρw)))

@inline ∂ⱼτ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline ∂ⱼτ₃ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶠ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, closure.ν, ρ, ρũ.ρw, ρũ.ρu)
+ δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, closure.ν, ρ, ρũ.ρw, ρũ.ρv)
+ δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, closure.ν, ρ, ρũ.ρw)))

@inline u₁∂ⱼτ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline u₁∂ⱼτ₁ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
U_over_ρ(i, j, k, grid, ρ, ρũ.ρu) * ∂ⱼτ₁ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline u₂∂ⱼτ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline u₂∂ⱼτ₂ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
V_over_ρ(i, j, k, grid, ρ, ρũ.ρv) * ∂ⱼτ₂ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline u₃∂ⱼτ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline u₃∂ⱼτ₃ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
W_over_ρ(i, j, k, grid, ρ, ρũ.ρw) * ∂ⱼτ₃ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)

@inline Q_dissipation(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
@inline Q_dissipation(i, j, k, grid, closure::IsotropicDiffusivity, ρ, ρũ, args...) =
( ℑxᶜᵃᵃ(i, j, k, grid, u₁∂ⱼτ₁ⱼ, closure, ρ, ρũ, args...)
+ ℑyᵃᶜᵃ(i, j, k, grid, u₂∂ⱼτ₂ⱼ, closure, ρ, ρũ, args...)
+ ℑzᵃᵃᶜ(i, j, k, grid, u₃∂ⱼτ₃ⱼ, closure, ρ, ρũ, args...))
11 changes: 6 additions & 5 deletions src/compressible_model.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
using Oceananigans
using Oceananigans.Models: AbstractModel, Clock, tracernames
using Oceananigans.Forcing: zeroforcing
using Oceananigans.Forcings: model_forcing

#####
##### Definition of a compressible model
Expand Down Expand Up @@ -48,9 +48,10 @@ function CompressibleModel(;
extra_tracers = nothing,
tracernames = collect_tracers(thermodynamic_variable, gases, microphysics, extra_tracers),
coriolis = nothing,
closure = ConstantIsotropicDiffusivity(float_type, ν=0.5, κ=0.5),
diffusivities = DiffusivityFields(architecture, grid, tracernames, closure),
forcing = ModelForcing(),
closure = IsotropicDiffusivity(float_type, ν=0.5, κ=0.5),
boundary_conditions = NamedTuple(),
diffusivities = DiffusivityFields(architecture, grid, tracernames, boundary_conditions, closure),
forcing = NamedTuple(),
gravity = g_Earth,
slow_forcings = ForcingFields(architecture, grid, tracernames),
right_hand_sides = RightHandSideFields(architecture, grid, tracernames),
Expand All @@ -59,7 +60,7 @@ function CompressibleModel(;

gravity = float_type(gravity)
tracers = TracerFields(architecture, grid, tracernames)
forcing = ModelForcing(tracernames, forcing)
forcing = model_forcing(tracernames; forcing...)
closure = with_tracers(tracernames, closure)
total_density = CellField(architecture, grid)

Expand Down
24 changes: 17 additions & 7 deletions src/time_stepping.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,9 @@ using JULES.Operators

using Oceananigans: datatuple
using Oceananigans.BoundaryConditions
import Oceananigans: time_step!

import Oceananigans.TimeSteppers: time_step!
import Oceananigans.Simulations: ab2_or_rk3_time_step!

#####
##### Utilities for time stepping
Expand All @@ -18,8 +20,10 @@ end
##### Time-stepping algorithm
#####

ab2_or_rk3_time_step!(model::CompressibleModel, Δt; euler) = time_step!(model, Δt)

# Adding kwargs... so this time_step! can work with Oceananigans.Simulation
function time_step!(model::CompressibleModel, Δt; kwargs...)
function time_step!(model::CompressibleModel, Δt)
ρ = model.total_density
ρũ = model.momenta
ρc̃ = model.tracers
Expand All @@ -43,13 +47,16 @@ function time_step!(model::CompressibleModel, Δt; kwargs...)

@debug "Computing slow forcings..."
update_total_density!(ρ, model.grid, model.gases, ρc̃)
fill_halo_regions!(merge((Σρ=ρ,), ρũ, ρc̃), model.architecture)
fill_halo_regions!(merge((Σρ=ρ,), ρũ, ρc̃), model.architecture, model.clock, nothing)

fill_halo_regions!(ρũ.ρw, model.architecture, model.clock, nothing)
fill_halo_regions!(IV_ρũ.ρw, model.architecture, model.clock, nothing)

compute_slow_forcings!(
F̃, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
model.coriolis, model.closure, ρ, ρũ, ρc̃, K̃, model.forcing, model.clock.time)
model.coriolis, model.closure, ρ, ρũ, ρc̃, K̃, model.forcing, model.clock)

fill_halo_regions!(F̃.ρw, model.architecture)
fill_halo_regions!(F̃.ρw, model.architecture, model.clock, nothing)

for rk3_iter in 1:3
@debug "RK3 step #$rk3_iter..."
Expand All @@ -60,15 +67,18 @@ function time_step!(model::CompressibleModel, Δt; kwargs...)
model.gases, model.gravity, ρ, ρũ, ρc̃, F̃)

update_total_density!(ρ, model.grid, model.gases, ρc̃)
fill_halo_regions!(merge((Σρ=ρ,), ρũ, ρc̃), model.architecture)
fill_halo_regions!(merge((Σρ=ρ,), ρũ, ρc̃), model.architecture, model.clock, nothing)
else
compute_rhs_args = (R̃, model.grid, model.thermodynamic_variable,
model.gases, model.gravity, ρ, IV_ρũ, IV_ρc̃, F̃)

update_total_density!(ρ, model.grid, model.gases, IV_ρc̃)
fill_halo_regions!(merge((Σρ=ρ,), IV_ρũ, IV_ρc̃), model.architecture)
fill_halo_regions!(merge((Σρ=ρ,), IV_ρũ, IV_ρc̃), model.architecture, model.clock, nothing)
end

fill_halo_regions!(ρũ.ρw, model.architecture, model.clock, nothing)
fill_halo_regions!(IV_ρũ.ρw, model.architecture, model.clock, nothing)

compute_right_hand_sides!(compute_rhs_args...)

@debug " Advancing variables..."
Expand Down
8 changes: 4 additions & 4 deletions src/time_stepping_kernels.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,12 @@ update_total_density!(model) =
"""
Slow forcings include viscous dissipation, diffusion, and Coriolis terms.
"""
function compute_slow_forcings!(F̃, grid, tvar, gases, gravity, coriolis, closure, ρ, ρũ, ρc̃, K̃, forcing, time)
function compute_slow_forcings!(F̃, grid, tvar, gases, gravity, coriolis, closure, ρ, ρũ, ρc̃, K̃, forcing, clock)
@inbounds begin
for k in 1:grid.Nz, j in 1:grid.Ny, i in 1:grid.Nx
F̃.ρu[i, j, k] = FU(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.u(i, j, k, grid, time, ρũ, ρc̃, nothing)
F̃.ρv[i, j, k] = FV(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.v(i, j, k, grid, time, ρũ, ρc̃, nothing)
F̃.ρw[i, j, k] = FW(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.w(i, j, k, grid, time, ρũ, ρc̃, nothing)
F̃.ρu[i, j, k] = FU(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.u(i, j, k, grid, clock, nothing)
F̃.ρv[i, j, k] = FV(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.v(i, j, k, grid, clock, nothing)
F̃.ρw[i, j, k] = FW(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) + forcing.w(i, j, k, grid, clock, nothing)
end

for (tracer_index, ρc_name) in enumerate(propertynames(ρc̃))
Expand Down
4 changes: 2 additions & 2 deletions test/test_lazy_fields.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ using JULES: LazyVelocityFields, LazyTracerFields
@testset "Lazy fields" begin
@info "Testing lazy fields..."

grid = RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
model = CompressibleModel(grid=grid, gases=DryEarth(),
thermodynamic_variable=Energy())

Expand All @@ -22,7 +22,7 @@ using JULES: LazyVelocityFields, LazyTracerFields

@test primitive_tracers.e[10, 11, 12] == 4/π

grid = RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
model = CompressibleModel(grid=grid, gases=DryEarth(),
thermodynamic_variable=Entropy())

Expand Down
4 changes: 2 additions & 2 deletions test/test_models.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
@testset "Energy thermodynamic variable" begin
@info " Testing model construction with energy thermodynamic variable..."

grid = RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
model = CompressibleModel(grid = grid, gases = DryEarth(),
thermodynamic_variable = Energy())
@test model isa CompressibleModel
Expand All @@ -13,7 +13,7 @@
@testset "Entropy thermodynamic variable" begin
@info " Testing model construction with entropy thermodynamic variable..."

grid = RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
model = CompressibleModel(grid = grid, gases = DryEarth(),
thermodynamic_variable = Entropy())
@test model isa CompressibleModel
Expand Down
21 changes: 14 additions & 7 deletions test/test_regression.jl
Original file line number Diff line number Diff line change
@@ -1,17 +1,20 @@
using Oceananigans.Fields: interiorparent

function summarize_regression_test(field_names, fields, correct_fields)
for (field_name, φ, φ_c) in zip(field_names, fields, correct_fields)
Δ = Array(φ) .- φ_c
function summarize_regression_test(fields, correct_fields)
for (field_name, φ, φ_c) in zip(keys(fields), fields, correct_fields)
Δ = φ .- φ_c

Δ_min = minimum(Δ)
Δ_max = maximum(Δ)
Δ_mean = mean(Δ)
Δ_abs_mean = mean(abs, Δ)
Δ_std = std(Δ)

@info(@sprintf("Δ%s: min=%.6e, max=%.6e, mean=%.6e, absmean=%.6e, std=%.6e",
field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std))
matching = sum(φ .≈ φ_c)
grid_points = length(φ_c)

@info @sprintf("Δ%s: min=%+.6e, max=%+.6e, mean=%+.6e, absmean=%+.6e, std=%+.6e (%d/%d matching grid points)",
field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std, matching, grid_points)
end
end

Expand Down Expand Up @@ -61,7 +64,9 @@ end
push!(correct_fields, file["ρs"])
end

summarize_regression_test(field_names, fields, correct_fields)
fields = NamedTuple{Tuple(field_names)}(Tuple(fields))
correct_fields = NamedTuple{Tuple(field_names)}(Tuple(correct_fields))
summarize_regression_test(fields, correct_fields)

@test all(isapprox.(interior(model.momenta.ρu), file["ρu"], atol=1e-12))

Expand Down Expand Up @@ -125,7 +130,9 @@ end
push!(correct_fields, file["ρs"])
end

summarize_regression_test(field_names, fields, correct_fields)
fields = NamedTuple{Tuple(field_names)}(Tuple(fields))
correct_fields = NamedTuple{Tuple(field_names)}(Tuple(correct_fields))
summarize_regression_test(fields, correct_fields)

@test all(isapprox.(interior(model.momenta.ρu), file["ρu"], atol=1e-12))

Expand Down
4 changes: 2 additions & 2 deletions test/test_time_stepping.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
@testset "Energy thermodynamic variable" begin
@info " Testing time stepping with energy thermodynamic variable..."

grid = RegularCartesianGrid(size=(16, 16, 16), halo=(2, 2, 2), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), halo=(2, 2, 2), extent=(1, 1, 1))
model = CompressibleModel(grid = grid, gases = DryEarth(),
thermodynamic_variable = Energy())
time_step!(model, 1)
Expand All @@ -14,7 +14,7 @@
@testset "Entropy thermodynamic variable" begin
@info " Testing time stepping with entropy thermodynamic variable..."

grid = RegularCartesianGrid(size=(16, 16, 16), halo=(2, 2, 2), length=(1, 1, 1))
grid = RegularCartesianGrid(size=(16, 16, 16), halo=(2, 2, 2), extent=(1, 1, 1))
model = CompressibleModel(grid = grid, gases = DryEarth(),
thermodynamic_variable = Entropy())
time_step!(model, 1)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ function simulate_dry_rising_thermal_bubble(; thermodynamic_variable, end_time=1
grid = grid,
gases = DryEarth(),
thermodynamic_variable = tvar,
closure = ConstantIsotropicDiffusivity(ν=75.0, κ=75.0)
closure = IsotropicDiffusivity(ν=75.0, κ=75.0)
)

#####
Expand Down Expand Up @@ -130,7 +130,7 @@ function simulate_dry_rising_thermal_bubble(; thermodynamic_variable, end_time=1
sim_parameters = (make_plots=make_plots, ρʰᵈ=ρʰᵈ, ρ̄ᵢ=ρ̄ᵢ, ρ̄s̄ᵢ=ρ̄s̄ᵢ)
end

simulation = Simulation(model, Δt=0.1, stop_time=end_time, progress_frequency=50,
simulation = Simulation(model, Δt=0.1, stop_time=end_time, iteration_interval=50,
progress=print_progress_and_make_plots, parameters=sim_parameters)
run!(simulation)

Expand Down Expand Up @@ -209,7 +209,7 @@ function print_progress_and_make_plots(simulation)

p = plot(u_plot, w_plot, ρ_plot, tvar_plot, layout=(2, 2), dpi=200, show=true)

n = Int(model.clock.iteration / simulation.progress_frequency)
n = Int(model.clock.iteration / simulation.iteration_interval)
n == 1 && !isdir("frames") && mkdir("frames")
savefig(p, @sprintf("frames/thermal_bubble_%s_%03d.png", typeof(tvar), n))
end
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ model = CompressibleModel(
grid = grid,
gases = DryEarth(),
thermodynamic_variable = tvar,
closure = ConstantIsotropicDiffusivity(ν=75.0, κ=75.0)
closure = IsotropicDiffusivity(ν=75.0, κ=75.0)
)

#####
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ function simulate_three_gas_dry_rising_thermal_bubble(;
grid = grid,
gases = DryEarth3(),
thermodynamic_variable = tvar,
closure = ConstantIsotropicDiffusivity(ν=75.0, κ=75.0)
closure = IsotropicDiffusivity(ν=75.0, κ=75.0)
)

#####
Expand Down Expand Up @@ -164,7 +164,7 @@ function simulate_three_gas_dry_rising_thermal_bubble(;
sim_parameters = (make_plots=make_plots, ρʰᵈ=ρʰᵈ, ρ̄ᵢ=ρ̄ᵢ, ρ̄s̄ᵢ=ρ̄s̄ᵢ)
end

simulation = Simulation(model, Δt=0.1, stop_time=end_time, progress_frequency=50,
simulation = Simulation(model, Δt=0.1, stop_time=end_time, iteration_interval=50,
progress=print_progress_and_make_plots, parameters=sim_parameters)
run!(simulation)

Expand Down Expand Up @@ -253,7 +253,7 @@ function print_progress_and_make_plots(simulation)
p = plot(ρ₁_plot, ρ₂_plot, ρ₃_plot, ρ_plot, ρ′_plot, tvar_plot,
layout=(2, 3), show=true, dpi=200)

n = Int(model.clock.iteration / simulation.progress_frequency)
n = Int(model.clock.iteration / simulation.iteration_interval)
n == 1 && !isdir("frames") && mkdir("frames")
savefig(p, @sprintf("frames/three_gas_thermal_bubble_%s_%03d.png", typeof(tvar), n))
end
Expand Down