advectiondiffusion.jl

### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ afcb23ea-f189-11ea-3b4c-1bcd4dc1e1cb
using DiffEqOperators, DifferentialEquations, Plots, PlutoUI, LinearAlgebra, Statistics, Measurements

# ╔═╡ c1975882-f19a-11ea-3ac1-b32fd1b87d2d
begin
	#pyplot()
	gr()
	ICbounds = [-5.0, 0.0]
	ICregion(x) = (ICbounds[1] ≤ x < ICbounds[2])
	nx = 1000
	t0 = 0.0
	tf = 5.0
	T(difo, advo, DIFF, ADV) = -advo * ADV + difo * DIFF
	dt = 0.01
	default(palette=palette(:Set1_9))
	md"IC/BC region and other shared parameters"
end

# ╔═╡ 90a08d04-f196-11ea-3d0c-77c2e58f093c
md"""
# Transition snapshot or equilibrium?
	
This notebook is about comparing the typical "spread" of a tracer coming from (but equivalently going to) a given region, with two distinct approaches:
- The first more "standard" approach is to **step a tracer in time** with a given initial condition (IC) and let it move around (and maybe decay) until a given time $\color{red}t$.
- The second approach, often adopted by OCIM users, is to impose the initial condition as a boundary condition instead (BC), let the tracer decay with a timescale $\color{blue}\tau$ everywhere else, and **directly solve for the steady-state**.  
	
This is what this notebook simulates. In the plot below, you can see the result of using both approaches (and a few other things). Use the sliders to change the parameters, the time *t* when the snapshot is taken, the decay timescale τ, and the advection and diffusion parameters.

#### Q: If $\color{red}t \color{black}= \color{blue}\tau$, which approach shows more "travel"?
"""

# ╔═╡ 6b64381e-f18d-11ea-05ae-19c6205e3999
md"""
t $(@bind t Slider(t0:dt:tf, default=0.0, show_value=true))\
τ $(@bind τ Slider(t0:dt:tf, default=1.0, show_value=true))

advection $(@bind advo Slider(0:100, default=25.0, show_value=true)) | diffusion $(@bind difo Slider(0.5:0.5:50, default=10.0, show_value=true))

show: advection $(@bind showadvo CheckBox()) | diffusion  $(@bind showdifo CheckBox())
"""

# ╔═╡ 4074da66-f189-11ea-0018-bd624fcd0bc6
begin	
	xbounds = (-10.0, 100.0)
	x1 = range(xbounds..., length=nx)
	h1 = step(x1)
	
	# Initial condition
	IC1 = ICregion.(x1)
	
	# transport (with BC)
	DIFF1 = CenteredDifference(2, 2, h1, nx)
	ADV1 = CenteredDifference(1, 2, h1, nx)
	BC1 = RobinBC((0.0,1.0,0.0), (0.0,1.0,0.0), h1, 1)
	transport1 = T(difo, advo, DIFF1, ADV1) * BC1
	# decay
	decay = DiffEqScalar(-1/τ)
	
	# Full operator, problem, and solution
	OP0 = AffineDiffEqOperator{Float64}((transport1,), (), 1.0IC1)
	OP1 = AffineDiffEqOperator{Float64}((transport1 + decay,), (), 1.0IC1)
	prob0 = ODEProblem(OP0, IC1, (t0, tf))
	prob1 = ODEProblem(OP1, IC1, (t0, tf))
	sol0 = solve(prob0, Tsit5())
	#sol1 = solve(prob1, Tsit5())
	
	# fast relaxation in IC region
	A_relax = DiffEqArrayOperator(-Diagonal(IC1/0.01))
	b_relax = IC1/0.01
	
	# Full operator, problem, and solution
	OP3 = AffineDiffEqOperator{Float64}((transport1 + decay, A_relax), (b_relax,), 1.0IC1)
	prob3 = ODEProblem(OP3, IC1, (t0, tf))
	#sol3 = solve(prob3, Tsit5())
	
	prob3SS = SteadyStateProblem(prob3)
	sol3SS = solve(prob3SS)
		
	md"Setup and solve the time-stepped problem"
end

# ╔═╡ 8c9a0a9c-f19a-11ea-045e-93340744a138
begin
	x2bounds = (0.0, 100.0)
	x2 = range(x2bounds..., length=nx)
	h2 = step(x2)
	
	# Initial condition 2 (because the x grid is different)
	IC2 = ICregion.(x2)

	# transport (with BC)
	DIFF2 = CenteredDifference(2, 2, h2, nx)
	ADV2 = CenteredDifference(1, 2, h2, nx)
	BC2 = RobinBC((1.0,0.0,1.0), (0.0,1.0,0.0), h2, 1)
	transport2 = T(difo, advo, DIFF2, ADV2) * BC2
	
	OP2 = AffineDiffEqOperator{Float64}((transport2 + decay,), (), 1.0IC2)
	
	prob2 = ODEProblem(OP2, IC2, (t0, tf))
	#sol2 = solve(prob2, Tsit5())
	
	prob2SS = SteadyStateProblem(prob2)
	#sol2SS = solve(prob2SS)
	
	md"Setup and solve the steady-state problem"
end

# ╔═╡ 9197d9d2-f18d-11ea-218e-ef1e770b0d76
begin
	# IC and BC
	xlims = (-10,110)
	plt = plot(x1, IC1, lab="Initial condition (IC)", c=:black, ylims=(-0.25,1.25), xlims=xlims, lw=1, yl="s(x, t=$t)", xl="x")
	#scatter!(plt, [x2[1]], [1.0], lab="Boundary condition (BC)", c=:black, m=:circle, ms=5)
	vspan!(plt, collect(ICbounds), α=0.1, c=:black, lab="IC region")
	# Annotations (t, τ, adv, and diff)
	#annotate!(plt, 70, 0.55, "t=$t")
	annotate!(plt, 25, 1, text("τ=$τ", 10))
	#annotate!(plt, 90, 0.55, "adv=$advo")
	#annotate!(plt, 90, 0.4, "diff=$difo")
	# Solution plots
	plot!(plt, x1, sol0(t), lab="s(t) (time-stepped solution)", lw=5, c=1)
	#plot!(plt, x1, sol1(t), lab="", lw=1, c=1, α=0.5)
	plot!(plt, x1, sol3SS, lab="s* (Steady-state solution)", c=2, lw=5)
	#plot!(plt, x2, sol2(t), lab="IC with BC (thick) or fast relaxation (thin)", α=0.5, lw=5, c=3)
	#plot!(plt, x1, sol3(t), lab="", lw=2, c=3)
	# t and τ * adv
	#vspan!(plt, advo * t .+ [0, -5], lab="(t = $t) × (adv = $advo)", c=1, α=0.1)
	#vspan!(plt, advo * τ .+ [0, -5], lab="(τ = $τ) × (adv = $advo)", c=2, α=0.1)
	showadvo && plot!(plt, mean(ICbounds) .+ [0, advo * t], 0.25/3 * [1,1], arrow=true, c=3, lab="advection displacement") 
	showdifo && plot!(plt, [mean(ICbounds)+ advo * t ± sqrt(2 * t * difo)], [2*0.25/3], lab="diffusion, √(2Dt)", c=5, msc=5) 
	plt2 = heatmap(x1, [0,1], (sol0(t) * [1, 1]')', axis=nothing, yl="s(t)", colorbar=false, xlims=xlims)
	plt3 = heatmap(x1, [0,1], (sol3SS.u * [1, 1]')', axis=nothing, yl="s*", yaxis=nothing, colorbar=false, xlims=xlims)
	showadvo && plot!(plt3, mean(ICbounds) .+ [0, advo * t], [1,1]/2, arrow=true, c=:white, lw=2, lab="")
	l = @layout [a{0.1h}
                 b{0.1h} 
				 c{0.8h} ]
	plot(plt2, plt3, plt, layout=l, link=:x)
end

# ╔═╡ Cell order:
# ╟─afcb23ea-f189-11ea-3b4c-1bcd4dc1e1cb
# ╟─c1975882-f19a-11ea-3ac1-b32fd1b87d2d
# ╟─4074da66-f189-11ea-0018-bd624fcd0bc6
# ╟─8c9a0a9c-f19a-11ea-045e-93340744a138
# ╟─90a08d04-f196-11ea-3d0c-77c2e58f093c
# ╟─9197d9d2-f18d-11ea-218e-ef1e770b0d76
# ╟─6b64381e-f18d-11ea-05ae-19c6205e3999