# This script is a (currently unsuccessful) attempt to implement a partial differential 
# equation describing changes in concentration with respect to longitudinal distance along 
# a stream channel. Changes in concentration result from both advection and a first-order 
# reaction in the hyporheic zone, which exchanges with the channel in such a way that 
# the ages of water returning to the channel are distributed as a power law. 

##### ----- Define Parameters, Initial Conditions, Boundary Conditions ----- #####
τ₋ = 1.0 # minimum residence time in seconds
τ₊ = 1000.0 # maximum residence time in seconds
α = 1.6 # power-law exponent
k = 0.0 # first-order rate constant in seconds^-1
S = 1.0 # volume of water in hyporheic zone
# Calculate rate of hyporheic exchange (volume per time)
q = qCalc_powerLaw(S, τ₋, τ₊, α)
v = 1 # Velocity
p = [τ₋, τ₊, α, k, q, v] # parameter vector
tspan = (0.0, 10.0) # simulation timespan
xₙ = 10 # physical length of simulation
N = 999 # number of interior points
Δx = xₙ/(N+1) # grid spacing along x dimension (longitudinal)
x = range(0, step = Δx, length = (N+2)) # x-points of solution, including boundaries
h(p,t) = zeros(N) # history function describing solute concentrations before the start of the simulation
C₀ = zeros(N) # initial condition
for i in 1:N
    C₀[i] = 1.0 - i*Δx/xₙ
end

##### ----- Finite Difference Operator Method ----- #####
ΔUpwind = DiffEqOperators.UpwindDifference(1,4,Δx,xₙ) # first derivative operator, fourth-order approximation
Dirichletbc = DiffEqOperators.DirichletBC(1.0,0.0) # boundary condition operator

@advectReact(ΔUpwind,Dirichletbc)

# Build the DDE problem
prob = DiffEqBase.ODEProblem(advectionReaction!, C₀, tspan, p)

# Solve problem
solution = solve(prob,Rosenbrock23())

# Plot results 
t = solution.t
u = solution.u

plt = plot()
for i in 1:length(t)
    plot!(x, Dirichletbc*u[[i]][1] ,label="t=$(t[i])")
end
display(plt)