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