# DiffeqOperators for non-uniform grid

**URL:** <https://discourse.julialang.org/t/diffeqoperators-for-non-uniform-grid/46802>\
**Category:** General Usage\
**Tags:** diffeq, pde\
**Created:** [September 17, 2020, 8:22pm UTC](https://discourse.julialang.org/t/diffeqoperators-for-non-uniform-grid/46802 "2020-09-17T20:22:53Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![JianghuiDu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jianghuidu/32/8028_2.png) [@JianghuiDu](https://discourse.julialang.org/u/JianghuiDu)\
**Post date:** [September 17, 2020, 8:22pm UTC](https://discourse.julialang.org/t/diffeqoperators-for-non-uniform-grid/46802/1 "2020-09-17T20:22:53Z")

</div>

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.

```julia
# 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

```

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [September 29, 2020, 1:57am UTC](https://discourse.julialang.org/t/diffeqoperators-for-non-uniform-grid/46802/2 "2020-09-29T01:57:11Z")

</div>

**EDIT: I just saw that there is already an issue on [DiffEqOperators.jl](https://github.com/SciML/DiffEqOperators.jl/issues/267), 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:

```julia
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:

 ![Screen Shot 2020-09-29 at 11.53.56 am](https://global.discourse-cdn.com/julialang/original/3X/9/c/9c6b305f08821ddfe0a5f83b220e8992d5c9cb19.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 29, 2020, 2:31pm UTC](https://discourse.julialang.org/t/diffeqoperators-for-non-uniform-grid/46802/3 "2020-09-29T14:31:43Z")

</div>

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.