Advection-diffusion of tracer in one dimension

This is an example demonstrating the advection-diffusion of a passive tracer in one dimension.

Install dependencies

First let's make sure we have all the required packages installed

using Pkg
pkg.add(["PassiveTracerFlows", "CairoMakie", "JLD2"])

Let's begin

First load packages needed to run this example.

using PassiveTracerFlows, CairoMakie, Printf, JLD2, LinearAlgebra

Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)

Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n²
stepper = "RK4"          # timestepper
     dt = 0.02           # timestep
 nsteps = 5000           # total number of time-steps

Physical parameters

L = 2π       # domain size
κ = 0.01     # diffusivity

Flow

We set a constant background flow and pass this to OneDAdvectingFlow with steadyflow = true to indicate the flow is not time dependent.

u(x) = 0.05
advecting_flow = OneDAdvectingFlow(; u, steadyflow = true)

OneDAdvectingFlow
  ├────────── u: Function
  └─ steadyflow: true

Problem setup

We initialize a Problem by providing a set of keyword arguments.

prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; nx=n, Lx=L, κ, dt, stepper)

and define some shortcuts.

sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x = grid.x

-3.141592653589793:0.04908738521234052:3.0925052683774528

Initial condition

We advect-diffuse a concentration field that has an initial concentration set to Gaussian.

gaussian(x, σ) = exp(-x^2 / 2σ^2)

amplitude, spread = 1, 0.15
c₀ = [amplitude * gaussian(x[i], spread) for i in 1:grid.nx]

TracerAdvectionDiffusion.set_c!(prob, c₀)

Saving output

We create the saved output using the Output function from FourierFlows.jl then save the concentration field using the get_concentration function every 50 timesteps.

function get_concentration(prob)
  ldiv!(prob.vars.c, prob.grid.rfftplan, deepcopy(prob.sol))

  return prob.vars.c
end

output = Output(prob, "advection-diffusion1D.jld2",
                (:concentration, get_concentration))

Output
  ├──── prob: FourierFlows.Problem{DataType, Vector{ComplexF64}, Float64, Vector{Float64}}
  ├──── path: advection-diffusion1D_1.jld2
  └── fields: Dict{Symbol, Function}(:concentration => Main.get_concentration)

By calling saveproblem(output) we save information that we will use for plotting later on.

saveproblem(output)

Stepping the problem forward

Now we step the problem forward and save output every 50 timesteps.

save_frequency = 50 # frequency at which output is saved

startwalltime = time()
while clock.step <= nsteps
  if clock.step % save_frequency == 0
    saveoutput(output)
    log = @sprintf("Output saved, step: %04d, t: %.2f, walltime: %.2f min",
                   clock.step, clock.t, (time()-startwalltime) / 60)

    println(log)
  end

  stepforward!(prob)
end

Output saved, step: 0000, t: 0.00, walltime: 0.00 min
Output saved, step: 0050, t: 1.00, walltime: 0.00 min
Output saved, step: 0100, t: 2.00, walltime: 0.00 min
Output saved, step: 0150, t: 3.00, walltime: 0.00 min
Output saved, step: 0200, t: 4.00, walltime: 0.00 min
Output saved, step: 0250, t: 5.00, walltime: 0.00 min
Output saved, step: 0300, t: 6.00, walltime: 0.00 min
Output saved, step: 0350, t: 7.00, walltime: 0.00 min
Output saved, step: 0400, t: 8.00, walltime: 0.00 min
Output saved, step: 0450, t: 9.00, walltime: 0.00 min
Output saved, step: 0500, t: 10.00, walltime: 0.00 min
Output saved, step: 0550, t: 11.00, walltime: 0.00 min
Output saved, step: 0600, t: 12.00, walltime: 0.00 min
Output saved, step: 0650, t: 13.00, walltime: 0.00 min
Output saved, step: 0700, t: 14.00, walltime: 0.00 min
Output saved, step: 0750, t: 15.00, walltime: 0.00 min
Output saved, step: 0800, t: 16.00, walltime: 0.00 min
Output saved, step: 0850, t: 17.00, walltime: 0.00 min
Output saved, step: 0900, t: 18.00, walltime: 0.00 min
Output saved, step: 0950, t: 19.00, walltime: 0.00 min
Output saved, step: 1000, t: 20.00, walltime: 0.00 min
Output saved, step: 1050, t: 21.00, walltime: 0.00 min
Output saved, step: 1100, t: 22.00, walltime: 0.00 min
Output saved, step: 1150, t: 23.00, walltime: 0.00 min
Output saved, step: 1200, t: 24.00, walltime: 0.00 min
Output saved, step: 1250, t: 25.00, walltime: 0.00 min
Output saved, step: 1300, t: 26.00, walltime: 0.00 min
Output saved, step: 1350, t: 27.00, walltime: 0.00 min
Output saved, step: 1400, t: 28.00, walltime: 0.00 min
Output saved, step: 1450, t: 29.00, walltime: 0.00 min
Output saved, step: 1500, t: 30.00, walltime: 0.00 min
Output saved, step: 1550, t: 31.00, walltime: 0.00 min
Output saved, step: 1600, t: 32.00, walltime: 0.00 min
Output saved, step: 1650, t: 33.00, walltime: 0.00 min
Output saved, step: 1700, t: 34.00, walltime: 0.00 min
Output saved, step: 1750, t: 35.00, walltime: 0.00 min
Output saved, step: 1800, t: 36.00, walltime: 0.00 min
Output saved, step: 1850, t: 37.00, walltime: 0.00 min
Output saved, step: 1900, t: 38.00, walltime: 0.00 min
Output saved, step: 1950, t: 39.00, walltime: 0.00 min
Output saved, step: 2000, t: 40.00, walltime: 0.00 min
Output saved, step: 2050, t: 41.00, walltime: 0.00 min
Output saved, step: 2100, t: 42.00, walltime: 0.00 min
Output saved, step: 2150, t: 43.00, walltime: 0.00 min
Output saved, step: 2200, t: 44.00, walltime: 0.00 min
Output saved, step: 2250, t: 45.00, walltime: 0.00 min
Output saved, step: 2300, t: 46.00, walltime: 0.00 min
Output saved, step: 2350, t: 47.00, walltime: 0.00 min
Output saved, step: 2400, t: 48.00, walltime: 0.00 min
Output saved, step: 2450, t: 49.00, walltime: 0.00 min
Output saved, step: 2500, t: 50.00, walltime: 0.00 min
Output saved, step: 2550, t: 51.00, walltime: 0.00 min
Output saved, step: 2600, t: 52.00, walltime: 0.00 min
Output saved, step: 2650, t: 53.00, walltime: 0.00 min
Output saved, step: 2700, t: 54.00, walltime: 0.00 min
Output saved, step: 2750, t: 55.00, walltime: 0.00 min
Output saved, step: 2800, t: 56.00, walltime: 0.00 min
Output saved, step: 2850, t: 57.00, walltime: 0.00 min
Output saved, step: 2900, t: 58.00, walltime: 0.00 min
Output saved, step: 2950, t: 59.00, walltime: 0.00 min
Output saved, step: 3000, t: 60.00, walltime: 0.00 min
Output saved, step: 3050, t: 61.00, walltime: 0.00 min
Output saved, step: 3100, t: 62.00, walltime: 0.00 min
Output saved, step: 3150, t: 63.00, walltime: 0.00 min
Output saved, step: 3200, t: 64.00, walltime: 0.00 min
Output saved, step: 3250, t: 65.00, walltime: 0.00 min
Output saved, step: 3300, t: 66.00, walltime: 0.00 min
Output saved, step: 3350, t: 67.00, walltime: 0.00 min
Output saved, step: 3400, t: 68.00, walltime: 0.00 min
Output saved, step: 3450, t: 69.00, walltime: 0.00 min
Output saved, step: 3500, t: 70.00, walltime: 0.00 min
Output saved, step: 3550, t: 71.00, walltime: 0.00 min
Output saved, step: 3600, t: 72.00, walltime: 0.00 min
Output saved, step: 3650, t: 73.00, walltime: 0.00 min
Output saved, step: 3700, t: 74.00, walltime: 0.00 min
Output saved, step: 3750, t: 75.00, walltime: 0.00 min
Output saved, step: 3800, t: 76.00, walltime: 0.00 min
Output saved, step: 3850, t: 77.00, walltime: 0.00 min
Output saved, step: 3900, t: 78.00, walltime: 0.00 min
Output saved, step: 3950, t: 79.00, walltime: 0.00 min
Output saved, step: 4000, t: 80.00, walltime: 0.00 min
Output saved, step: 4050, t: 81.00, walltime: 0.00 min
Output saved, step: 4100, t: 82.00, walltime: 0.00 min
Output saved, step: 4150, t: 83.00, walltime: 0.00 min
Output saved, step: 4200, t: 84.00, walltime: 0.00 min
Output saved, step: 4250, t: 85.00, walltime: 0.00 min
Output saved, step: 4300, t: 86.00, walltime: 0.00 min
Output saved, step: 4350, t: 87.00, walltime: 0.00 min
Output saved, step: 4400, t: 88.00, walltime: 0.00 min
Output saved, step: 4450, t: 89.00, walltime: 0.00 min
Output saved, step: 4500, t: 90.00, walltime: 0.00 min
Output saved, step: 4550, t: 91.00, walltime: 0.00 min
Output saved, step: 4600, t: 92.00, walltime: 0.00 min
Output saved, step: 4650, t: 93.00, walltime: 0.00 min
Output saved, step: 4700, t: 94.00, walltime: 0.00 min
Output saved, step: 4750, t: 95.00, walltime: 0.00 min
Output saved, step: 4800, t: 96.00, walltime: 0.00 min
Output saved, step: 4850, t: 97.00, walltime: 0.00 min
Output saved, step: 4900, t: 98.00, walltime: 0.00 min
Output saved, step: 4950, t: 99.00, walltime: 0.00 min
Output saved, step: 5000, t: 100.00, walltime: 0.00 min

Visualising the output

We load the .jld2 file and create a timeseries of the concentration field

file = jldopen(output.path)

iterations = parse.(Int, keys(file["snapshots/t"]))

t = [file["snapshots/t/$i"] for i ∈ iterations]
c = [file["snapshots/concentration/$i"] for i ∈ iterations]

Set up the plotting arguments and look at the initial concentration.

x, Lx = file["grid/x"], file["grid/Lx"]

n = Observable(1)
c_anim = @lift Array(c[$n])
title = @lift @sprintf("concentration, t = %s", t[$n])

fig = Figure(size = (600, 600))
ax = Axis(fig[1, 1],
          xlabel = "x",
          ylabel = "c",
          limits = ((-Lx/2, Lx/2), (0, maximum(c[1]))))

lines!(ax, x, c_anim; linewidth = 4)

Lines{Tuple{Vector{Point{2, Float32}}}}

Now, we create a movie of the tracer concentration being advected and diffused.

frames = 1:length(t)
record(fig, "1D_advection-diffusion.mp4", frames, framerate = 18) do i
    n[] = i
end

This page was generated using Literate.jl.