DiffeqOperators for non-uniform grid

[JianghuiDu]

I have a problem of using DiffeqOperators.jl on non-uniform grid. Here is an example of a simple diffusion equation D\nabla^2C+1=0 between [0,50] with Robin boundary conditions on both side. It is solved using centered spatial difference and Method of Line time-stepping towards steady state. I did uniform grid first and then a non-uniform grid later. Compared to the analytical solution. The uniform grid results are right but not the non-uniform one. It seems the upper boundary condition is wrong when using non-uniform grid.

# grid set up
L = 50.0 # domain [0,L]

N = 1000 # number of grid points

# a uniform grid
h = L/(N+1) # spatial steps
x = collect(range(h, step=h, length=N))
dx = vcat(x,L).-vcat(0,x)

# a non-uniform grid (finer mesh at the boundaries)
x = L*sin.(collect(range(h, step=h, length=N))/L*π/2).^2
dx = vcat(x,L).-vcat(0,x)

# a simple diffusion equation
D=10.0
Δ2 = CenteredDifference(2, 2,dx, N);
bc = RobinBC((0,-D,10.0),(1.0,0.0,0.0),dx,1);

function du_operator(du,u,parm,t)
    du .=  D.*(Δ2*bc*u) .+1
end

prob_operator = ODEProblem(du_operator,zeros(N),(0.0,10e3));
cb = TerminateSteadyState(1e-12, 1e-9,DiffEqCallbacks.allDerivPass);
@time sol = solve(prob_operator,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1),
            reltol=1e-9,abstol=1e-12,callback=cb,tstops=10.0);

# analytical solution
analytical_sol =  175 .- x .- x.^2/20

[briochemc]

EDIT: I just saw that there is already an issue on DiffEqOperators.jl, so apologies for the noise!

---

@ChrisRackauckas Could you take a look at this? Can the boundary condition operators not handle non-uniform grids? Maybe there is some issue here with the location of the ghost points? Not sure…

Anyway, it seems @JianghuiDu is correct that there is an issue here, so I rewrote the same problem below as a self-contained MWE and added plots of the solutions, maybe that will help someone identify the issue?

Here is a similar MWE:

using DifferentialEquations, DiffEqOperators, Sundials, Plots
L = 50.0 # domain [0,L]
N = 1000 # number of grid points
# a uniform grid
h = L / (N + 1) # spatial steps
x1 = collect(range(h, step=h, length=N))
# a non-uniform grid (finer mesh at the boundaries)
x2 = L * sin.(collect(range(h, step=h, length=N)) / L * π / 2) .^ 2
# dx vector
dx(x) = [x; L] - [0.0; x]
# a simple diffusion equation
D = 10.0
Δ(x) = CenteredDifference(2, 2, dx(x), N)
# Boundary condition
bc(x) = RobinBC((0.0, -D, 10.0), (1.0, 0.0, 0.0), dx(x), 1)
# full operator is an affine operator so let's use that
op(x) = AffineDiffEqOperator{Float64}((D * Δ(x) * bc(x),), (1.0,), zeros(N))
# problem
prob(x) = ODEProblem(op(x), ones(N), (0.0,10e3))
# callback
cb = TerminateSteadyState(1e-12, 1e-9, DiffEqCallbacks.allDerivPass)
# solution as a function of the x-grid
sol(x) = solve(prob(x), CVODE_BDF(linear_solver=:Band, jac_upper=1, jac_lower=1),
            reltol=1e-9, abstol=1e-12, callback=cb, tstops=10.0)
# compute both solutions
sol1 = sol(x1)
sol2 = sol(x2)
# analytical solution
analytical_sol(x) = 175 - x - x^2 / 20
# plot
plot(x1, analytical_sol, lw=10, lab="analytical solution", xlab="x", bg=:black, α=0.5)
plot!(x1, sol1(sol1.t[end]), lw=3, lab="uniform-grid solution")
plot!(x2, sol2(sol2.t[end]), lw=3, lab="non-uniform-grid solution")

and this is the plot I get:

[ChrisRackauckas]

No worries, thanks for bringing it up. Yes, we can follow up on the repo. I think it might be a boundary condition implementation issue since that is what would cause the left to drop. Non-uniforms grids must not be properly taken into account at the boundary, making it leak.