Merge pull request #85 from thabbott/ar/higher-order-advection

Higher-order advection schemes

src/Operators/Operators.jl

@@ -10,8 +10,8 @@ using Oceananigans.TurbulenceClosures: IsotropicDiffusivity
 export
     kinetic_energy,
     hdivᶜᶜᵃ, divᶜᶜᶜ, ∇²,
-    div_uc, ∂ⱼpuⱼ, ∂ⱼDᶜⱼ, ∂ⱼtᶜDᶜⱼ, ∂ⱼDᵖⱼ,
-    div_ρuũ, div_ρvũ, div_ρwũ,
+    div_ρuũ, div_ρvũ, div_ρwũ, div_ρUc,
+    ∂ⱼpuⱼ, ∂ⱼDᶜⱼ, ∂ⱼtᶜDᶜⱼ, ∂ⱼDᵖⱼ,
     ∂ⱼτ₁ⱼ, ∂ⱼτ₂ⱼ, ∂ⱼτ₃ⱼ, Q_dissipation
 
 include("compressible_operators.jl")

src/Operators/compressible_operators.jl

@@ -1,3 +1,5 @@
+using Oceananigans.Advection
+
 #####
 ##### Legend:
 ##### ρ  -> density fields
@@ -52,15 +54,15 @@
 ##### Tracer advection
 #####
 
-@inline advective_tracer_flux_x(i, j, k, grid, ρ, ρu, ρc) = Ax_ψᵃᵃᶠ(i, j, k, grid, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρc) / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
-@inline advective_tracer_flux_y(i, j, k, grid, ρ, ρv, ρc) = Ay_ψᵃᵃᶠ(i, j, k, grid, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρc) / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
-@inline advective_tracer_flux_z(i, j, k, grid, ρ, ρw, ρc) = Az_ψᵃᵃᵃ(i, j, k, grid, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρc) / ℑzᵃᵃᶠ(i, j, k, grid, ρ)
+@inline tracer_flux_ρuc(i, j, k, grid, scheme, ρ, ρu, ρc) = advective_tracer_flux_x(i, j, k, grid, scheme, ρu, ρc) / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
+@inline tracer_flux_ρvc(i, j, k, grid, scheme, ρ, ρv, ρc) = advective_tracer_flux_y(i, j, k, grid, scheme, ρv, ρc) / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
+@inline tracer_flux_ρwc(i, j, k, grid, scheme, ρ, ρw, ρc) = advective_tracer_flux_z(i, j, k, grid, scheme, ρw, ρc) / ℑzᵃᵃᶠ(i, j, k, grid, ρ)
 
-@inline div_uc(i, j, k, grid, ρ, ρũ, ρc) =
+@inline div_ρUc(i, j, k, grid, scheme, density, momenta, ρc) =
     (1/Vᵃᵃᶜ(i, j, k, grid)
-        * (  δxᶜᵃᵃ(i, j, k, grid, advective_tracer_flux_x, ρ, ρũ.ρu, ρc)
-           + δyᵃᶜᵃ(i, j, k, grid, advective_tracer_flux_y, ρ, ρũ.ρv, ρc)
-           + δzᵃᵃᶜ(i, j, k, grid, advective_tracer_flux_z, ρ, ρũ.ρw, ρc)))
+        * (  δxᶜᵃᵃ(i, j, k, grid, tracer_flux_ρuc, scheme, density, momenta.ρu, ρc)
+           + δyᵃᶜᵃ(i, j, k, grid, tracer_flux_ρvc, scheme, density, momenta.ρv, ρc)
+           + δzᵃᵃᶜ(i, j, k, grid, tracer_flux_ρwc, scheme, density, momenta.ρw, ρc)))
 
 #####
 ##### Diffusion
@@ -134,34 +136,34 @@ end
 ##### Momentum advection
 #####
 
-@inline momentum_flux_ρuu(i, j, k, grid, ρ, ρu) = @inbounds ℑxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶜᵃᵃ(i, j, k, grid, ρu) / ρ[i, j, k]
-@inline momentum_flux_ρvv(i, j, k, grid, ρ, ρv) = @inbounds ℑyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶜᵃ(i, j, k, grid, ρv) / ρ[i, j, k]
-@inline momentum_flux_ρww(i, j, k, grid, ρ, ρw) = @inbounds ℑzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶜ(i, j, k, grid, ρw) / ρ[i, j, k]
+@inline momentum_flux_ρuu(i, j, k, grid, scheme, ρ, ρu) = @inbounds momentum_flux_uu(i, j, k, grid, scheme, ρu) / ρ[i, j, k]
+@inline momentum_flux_ρvv(i, j, k, grid, scheme, ρ, ρv) = @inbounds momentum_flux_vv(i, j, k, grid, scheme, ρv) / ρ[i, j, k]
+@inline momentum_flux_ρww(i, j, k, grid, scheme, ρ, ρw) = @inbounds momentum_flux_ww(i, j, k, grid, scheme, ρw) / ρ[i, j, k]
 
-@inline momentum_flux_ρuv(i, j, k, grid, ρ, ρu, ρv) = ℑxᶠᵃᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρu) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
-@inline momentum_flux_ρuw(i, j, k, grid, ρ, ρu, ρw) = ℑxᶠᵃᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρu) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
-@inline momentum_flux_ρvu(i, j, k, grid, ρ, ρu, ρv) = ℑyᵃᶠᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρv) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
-@inline momentum_flux_ρvw(i, j, k, grid, ρ, ρv, ρw) = ℑyᵃᶠᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρv) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
-@inline momentum_flux_ρwu(i, j, k, grid, ρ, ρu, ρw) = ℑzᵃᵃᶠ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρw) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
-@inline momentum_flux_ρwv(i, j, k, grid, ρ, ρv, ρw) = ℑzᵃᵃᶠ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρw) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
+@inline momentum_flux_ρuv(i, j, k, grid, scheme, ρ, ρu, ρv) = momentum_flux_uv(i, j, k, grid, scheme, ρu, ρv) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
+@inline momentum_flux_ρuw(i, j, k, grid, scheme, ρ, ρu, ρw) = momentum_flux_uw(i, j, k, grid, scheme, ρu, ρw) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
+@inline momentum_flux_ρvu(i, j, k, grid, scheme, ρ, ρu, ρv) = momentum_flux_vu(i, j, k, grid, scheme, ρu, ρv) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
+@inline momentum_flux_ρvw(i, j, k, grid, scheme, ρ, ρv, ρw) = momentum_flux_vw(i, j, k, grid, scheme, ρv, ρw) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
+@inline momentum_flux_ρwu(i, j, k, grid, scheme, ρ, ρu, ρw) = momentum_flux_wu(i, j, k, grid, scheme, ρu, ρw) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
+@inline momentum_flux_ρwv(i, j, k, grid, scheme, ρ, ρv, ρw) = momentum_flux_wv(i, j, k, grid, scheme, ρv, ρw) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
 
-@inline div_ρuũ(i, j, k, grid, ρ, ρũ) =
+@inline div_ρuũ(i, j, k, grid, scheme, ρ, ρũ) =
     (1/Vᵃᵃᶜ(i, j, k, grid)
-        * (  δxᶠᵃᵃ(i, j, k, grid, momentum_flux_ρuu, ρ, ρũ.ρu)
-           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρuv, ρ, ρũ.ρu, ρũ.ρv)
-           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρuw, ρ, ρũ.ρu, ρũ.ρw)))
+        * (  δxᶠᵃᵃ(i, j, k, grid, momentum_flux_ρuu, scheme, ρ, ρũ.ρu)
+           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρuv, scheme, ρ, ρũ.ρu, ρũ.ρv)
+           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρuw, scheme, ρ, ρũ.ρu, ρũ.ρw)))
 
-@inline div_ρvũ(i, j, k, grid, ρ, ρũ) =
+@inline div_ρvũ(i, j, k, grid, scheme, ρ, ρũ) =
     (1/Vᵃᵃᶜ(i, j, k, grid)
-        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρvu, ρ, ρũ.ρu, ρũ.ρv)
-           + δyᵃᶠᵃ(i, j, k, grid, momentum_flux_ρvv, ρ, ρũ.ρv)
-           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρvw, ρ, ρũ.ρv, ρũ.ρw)))
+        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρvu, scheme, ρ, ρũ.ρu, ρũ.ρv)
+           + δyᵃᶠᵃ(i, j, k, grid, momentum_flux_ρvv, scheme, ρ, ρũ.ρv)
+           + δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρvw, scheme, ρ, ρũ.ρv, ρũ.ρw)))
 
-@inline div_ρwũ(i, j, k, grid, ρ, ρũ) =
+@inline div_ρwũ(i, j, k, grid, scheme, ρ, ρũ) =
     (1/Vᵃᵃᶠ(i, j, k, grid)
-        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρwu, ρ, ρũ.ρu, ρũ.ρw)
-           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρwv, ρ, ρũ.ρv, ρũ.ρw)
-           + δzᵃᵃᶠ(i, j, k, grid, momentum_flux_ρww, ρ, ρũ.ρw)))
+        * (  δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρwu, scheme, ρ, ρũ.ρu, ρũ.ρw)
+           + δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρwv, scheme, ρ, ρũ.ρv, ρũ.ρw)
+           + δzᵃᵃᶠ(i, j, k, grid, momentum_flux_ρww, scheme, ρ, ρũ.ρw)))
 
 #####
 ##### Viscous dissipation

src/compressible_model.jl

@@ -1,15 +1,17 @@
 using Oceananigans
+using Oceananigans.Advection: CenteredSecondOrder
 using Oceananigans.Models: AbstractModel, Clock, tracernames
 using Oceananigans.Forcings: model_forcing
 
 #####
 ##### Definition of a compressible model
 #####
 
-mutable struct CompressibleModel{A, FT, Ω, D, M, V, T, L, K, Θ, G, X, C, F, S} <: AbstractModel
+mutable struct CompressibleModel{A, FT, Ω, ∇, D, M, V, T, L, K, Θ, G, X, C, F, S} <: AbstractModel
               architecture :: A
                       grid :: Ω
                      clock :: Clock{FT}
+                 advection :: ∇
              total_density :: D
                    momenta :: M
                 velocities :: V
@@ -37,6 +39,7 @@ function CompressibleModel(;
               architecture = CPU(),
                 float_type = Float64,
                      clock = Clock{float_type}(0, 0),
+                 advection = CenteredSecondOrder(),
                    momenta = MomentumFields(architecture, grid),
                      gases = DryEarth(float_type),
     thermodynamic_variable = Energy(),
@@ -59,7 +62,7 @@ function CompressibleModel(;
     lazy_tracers = LazyTracerFields(architecture, grid, total_density, tracers)
 
     return CompressibleModel(
-        architecture, grid, clock, total_density, momenta, velocities, tracers,
+        architecture, grid, clock, advection, total_density, momenta, velocities, tracers,
         lazy_tracers, diffusivities, thermodynamic_variable, gases, gravity,
         coriolis, closure, forcing, time_stepper)
 end

src/source_terms.jl

@@ -38,38 +38,38 @@ using Oceananigans.Coriolis
     end
 end
 
-@inline function ρu_slow_source_term(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃, FU)
+@inline function ρu_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FU)
     @inbounds begin
-        return (- div_ρuũ(i, j, k, grid, ρ, ρũ)
+        return (- div_ρuũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                 - ∂p∂x(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                 + FU[i, j, k])
     end
 end
 
-@inline function ρv_slow_source_term(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃, FV)
+@inline function ρv_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FV)
     @inbounds begin
-        return (- div_ρvũ(i, j, k, grid, ρ, ρũ)
+        return (- div_ρvũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                 - ∂p∂y(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                 + FV[i, j, k])
     end
 end
 
-@inline function ρw_slow_source_term(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃, FW)
+@inline function ρw_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FW)
     @inbounds begin
-        return (- div_ρwũ(i, j, k, grid, ρ, ρũ)
+        return (- div_ρwũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                 - ∂p∂z(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                 - gravity * ℑzᵃᵃᶠ(i, j, k, grid, ρ)
                 + FW[i, j, k])
     end
 end
 
-@inline function ρc_slow_source_term(i, j, k, grid, ρ, ρũ, ρc, FC)
+@inline function ρc_fast_source_term(i, j, k, grid, advection_scheme, ρ, ρũ, ρc, FC)
     @inbounds begin
-        return -div_uc(i, j, k, grid, ρ, ρũ, ρc) + FC[i, j, k]
+        return -div_ρUc(i, j, k, grid, advection_scheme, ρ, ρũ, ρc) + FC[i, j, k]
     end
 end
 
-@inline ρc_slow_source_term(i, j, k, grid::AbstractGrid{T}, tvar::Entropy, gases, gravity, ρ, ρũ, ρc̃) where T = zero(T)
+@inline ρt_fast_source_term(i, j, k, grid::AbstractGrid{T}, tvar::Entropy, gases, gravity, ρ, ρũ, ρc̃) where T = zero(T)
 
-@inline ρc_slow_source_term(i, j, k, grid, tvar::Energy, gases, gravity, ρ, ρũ, ρc̃) =
+@inline ρt_fast_source_term(i, j, k, grid, tvar::Energy, gases, gravity, ρ, ρũ, ρc̃) =
     -∂ⱼpuⱼ(i, j, k, grid, diagnose_pressure, tvar, gases, gravity, ρ, ρũ, ρc̃)

src/time_stepping.jl

@@ -58,8 +58,8 @@ function time_step!(model::CompressibleModel, Δt)
     fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
     fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)
 
-    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable,
-                               model.gases, model.gravity, total_density, momenta, tracers, slow_source_terms)
+    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
+                               model.advection, total_density, momenta, tracers, slow_source_terms)
 
     advance_variables!(intermediate_fields, model.grid, momenta, tracers, fast_source_terms, Δt=first_stage_Δt)
 
@@ -75,8 +75,8 @@ function time_step!(model::CompressibleModel, Δt)
     fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
     fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)
 
-    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable,
-                               model.gases, model.gravity, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)
+    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
+                               model.advection, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)
 
     advance_variables!(intermediate_fields, model.grid, momenta, tracers, fast_source_terms, Δt=second_stage_Δt)
 
@@ -94,8 +94,8 @@ function time_step!(model::CompressibleModel, Δt)
     fill_halo_regions!(momenta.ρw, model.architecture, model.clock, nothing)
     fill_halo_regions!(intermediate_momenta.ρw, model.architecture, model.clock, nothing)
 
-    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable,
-                               model.gases, model.gravity, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)
+    compute_fast_source_terms!(fast_source_terms, model.grid, model.thermodynamic_variable, model.gases, model.gravity,
+                               model.advection, total_density, intermediate_momenta, intermediate_tracers, slow_source_terms)
 
     state_variables = (momenta..., tracers = tracers)
     advance_variables!(state_variables, model.grid, momenta, tracers, fast_source_terms, Δt=third_stage_Δt)

src/time_stepping_kernels.jl

@@ -42,12 +42,12 @@ end
 """
 Fast forcings include advection, pressure gradient, and buoyancy terms.
 """
-function compute_fast_source_terms!(fast_source_terms, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, slow_source_terms)
+function compute_fast_source_terms!(fast_source_terms, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms)
     @inbounds begin
         for k in 1:grid.Nz, j in 1:grid.Ny, i in 1:grid.Nx
-            fast_source_terms.ρu[i, j, k] = ρu_slow_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, slow_source_terms.ρu)
-            fast_source_terms.ρv[i, j, k] = ρv_slow_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, slow_source_terms.ρv)
-            fast_source_terms.ρw[i, j, k] = ρw_slow_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers, slow_source_terms.ρw)
+            fast_source_terms.ρu[i, j, k] = ρu_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρu)
+            fast_source_terms.ρv[i, j, k] = ρv_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρv)
+            fast_source_terms.ρw[i, j, k] = ρw_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, advection_scheme, total_density, momenta, tracers, slow_source_terms.ρw)
         end
 
         for ρc_name in propertynames(tracers)
@@ -56,14 +56,13 @@ function compute_fast_source_terms!(fast_source_terms, grid, thermodynamic_varia
             S_ρc = getproperty(slow_source_terms.tracers, ρc_name)
 
             for k in 1:grid.Nz, j in 1:grid.Ny, i in 1:grid.Nx
-                F_ρc[i, j, k] = ρc_slow_source_term(i, j, k, grid, total_density, momenta, ρc, S_ρc)
+                F_ρc[i, j, k] = ρc_fast_source_term(i, j, k, grid, advection_scheme, total_density, momenta, ρc, S_ρc)
             end
         end
 
         for k in 1:grid.Nz, j in 1:grid.Ny, i in 1:grid.Nx
-            fast_source_terms.tracers[1].data[i, j, k] += ρc_slow_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
+            fast_source_terms.tracers[1].data[i, j, k] += ρt_fast_source_term(i, j, k, grid, thermodynamic_variable, gases, gravity, total_density, momenta, tracers)
         end
-
     end
 end
 

verification/periodic_advection/periodic_advection.jl

@@ -2,24 +2,29 @@ using Printf
 using Plots
 using JULES
 using Oceananigans
+using Oceananigans.Advection
 
+using JULES: IdealGas
 using Oceananigans.Grids: Cell, xnodes
 
 ENV["GKSwstype"] = "100"
 
 N = 64
 L = 1
-T = 1000
-Δt = 1e1
+T = 2
+Δt = 1e-2
 Nt = Int(T/Δt)
 
+ideal_gas = IdealGas(Float64, JULES.R⁰, 1; T₀=1, p₀=1, s₀=1, u₀=0)
+
 topo = (Periodic, Periodic, Periodic)
 domain = (x=(-L/2, L/2), y=(0, 1), z=(0, 1))
-grid = RegularCartesianGrid(topology=topo, size=(N, 1, 1), halo=(2, 2, 2); domain...)
+grid = RegularCartesianGrid(topology=topo, size=(N, 1, 1), halo=(4, 4, 4); domain...)
 
 model = CompressibleModel(
                       grid = grid,
-                     gases = DryEarth(),
+                 advection = WENO5(),
+                     gases = (ρ=ideal_gas,),
     thermodynamic_variable = Entropy(),
                    closure = IsotropicDiffusivity(ν=0, κ=0),
                    gravity = 0.0,
@@ -43,7 +48,7 @@ update_total_density!(model)
 @inline x′(x, t, L, U) = mod(x + L/2 - U * t, L) - L/2
 @inline ϕ_Gaussian(x, t; L, U, a=1, c=1/8) = a * exp(-x′(x, t, L, U)^2 / (2c^2))
 @inline ϕ_Square(x, t; L, U, w=0.15) = -w <= x′(x, t, L, U) <= w ? 1.0 : 0.0
-ρc₀(x, y, z) = ϕ_Gaussian(x, 0, L=L, U=1)
+ρc₀(x, y, z) = ϕ_Square(x, 0, L=L, U=1)
 set!(model.tracers.ρc, ρc₀)
 
 anim = @animate for n in 1:Nt
@@ -54,7 +59,7 @@ anim = @animate for n in 1:Nt
 
     x = xnodes(Cell, grid)
     ρc = interior(model.tracers.ρc)[:]
-    plot(x, ρc, lw=2, label="", title=title, xlims=(-L/2, L/2), ylims=(0, 1), dpi=200)
+    plot(x, ρc, lw=2, label="", title=title, xlims=(-L/2, L/2), ylims=(-0.2, 1.1), dpi=200)
 end
 
 mp4(anim, "periodic_advection.mp4", fps=15)