Hi everyone,
For one of my problem, I need to solve this ODE backwards:
where X is the position of a marker (could be 1d, 2d or 3d), a(X(t),t) is a velocity vector and X_g is a grid point.
The goal is to find the position of X at the previous timestep n knowing its current position at n+1 following characteristics and knowing a(X(t),t) at n and n+1 (this is for solving an advection problem using a semi-Lagrangian approach).
It is common to solve this equation using the implicit mid-point rule:
X(t^n) = - a(\frac{X(t^{n+1}) + X(t^n)}{2}, t^{n+\frac{1}{2}}) \Delta t + X(t^{n+1})
and iterating until convergence is reached as x^n is on both sides.
What complicates that is that a(X,t) doesn’t have an analytical expression, so a(\frac{X^{n+1} + x^{n}}{2}, t^{n+\frac{1}{2}}) needs to be interpolated from a(X, t^n) and a(X, t^{n+1}) in each iteration. So the final expression is
X(t^n)^{k+1} = - \mathcal{L}(\frac{X_g + X(t^n)^k}{2}, t^{n+\frac{1}{2}}) \Delta t + X_g, \: \textrm{for} \: k \ge 0, X(t^n)^{0} = X_g
with \mathcal{L} being a first order Lagrange interpolant for example.
So, I know how to program that, no problem. My issue is that I don’t know how I could combine that with OrdinaryDiffEq. My goal to do so is to try higher order method such as the implicit Adams–Moulton scheme and that would be nice to not have to implement that myself.
Here is a MWE of the approach “by hand”:
using Plots
using LinearAlgebra: norm
# grid
nx::Int = 51
Lx::Float64 = 1.0
Δx::Float64 = Lx / (nx-1)
x::StepRangeLen = range(0, length=nx, stop= Lx)
# velocity
v::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.5) .+ 1.0
v_timestep::Vector{Float64} = ones(nx) .* cos.(range(0, length=nx, stop= 10*Lx) .- 0.2) .+ 1.0
v_t_half::Vector{Float64} = zeros(nx)
# CFL condition
CFL::Float64 = 1.5
period::Int = 2
tmax::Float64 = period * Lx / maximum(abs.(v))
Δt::Float64 = CFL * Δx / maximum(abs.(v))
X_t0::Vector{Float64} = zeros(nx)
X_t_half::Vector{Float64} = zeros(nx)
function linear_interpolation(x, y, new_x)
idx = searchsortedlast(x, new_x)
if idx == 0
return y[1]
elseif idx == length(x)
return y[end]
else
frac = (new_x - x[idx]) / (x[idx+1] - x[idx])
return y[idx] + frac * (y[idx+1] - y[idx])
end
end
plt = scatter(x, zeros(size(x)))
plot!(plt, x, v, label="v")
plot!(plt, x, v_timestep, label="v_timestep")
v_t_half .= (v .+ v_timestep) / 2
# save the initial position
X_t_old = zeros(size(x))
# iterate 5 times to improve approximation
for i in 1:5
X_t_old = copy(X_t_half)
X_t_half .= x .- v_t_half .* Δt ./ 2
v_t_half .= [linear_interpolation(x, v_t_half, xi) for xi in X_t_half]
display(scatter!(X_t_half[1:10:end], zeros(size(X_t_half[1:10:end])), label=""))
# compute error between X_t_half and X_t_old
error = norm(X_t_half - X_t_old)
println("Error: $error")
end
I fail to convert this code to something compatible with OrdinaryDiffEq.
I guess it should take something in that form:
using OrdinaryDiffEq
p = (v=v, v_timestep=v_timestep)
function system!(du, u, p, t)
v, v_timestep = p
du .= v
end
u0 = collect(x)
tspan = (0.0, -Δt)
prob = ODEProblem(system!, u0, tspan, dt = -Δt)
But I don’t know how I could tell it to interpolate v or if I should maybe provide this as a parameter.
