# Using DifferentialEquations.jl with ghost cells

**URL:** <https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561>\
**Category:** General Usage\
**Created:** [June 11, 2026, 11:11am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561 "2026-06-11T11:11:31Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lucas](https://avatars.discourse-cdn.com/v4/letter/l/ba8739/32.png) [@Lucas](https://discourse.julialang.org/u/Lucas)\
**Post date:** [June 11, 2026, 11:11am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561/1 "2026-06-11T11:11:32Z")

</div>

I’m trying to solve a few differential equations using the method of lines. I’m using a five-point laplacian and ghost cells. I will add some code below, for the case of an imaginary-time Gross-Pitaevskii equation with chemical potential.

This solver works fine with an explicit method like Tsit5(), but doesn’t with implicit solvers. One issue is that I used OffsetArrays, which the standard automatic differentiation backend doesn’t like. But even using finite differences, it is not working; I suspect because the ghost cells make the Jacobian singular.

How can this problem be avoided? I know that in this case, it is easy to omit the ghost cells and replace them by a separate stencil at the edges, but I would like to implement more complicated boundary conditions in the future, in which the ghost cell approach is most practical.

MWE:

```julia

using DifferentialEquations
using OffsetArrays
using NonlinearSolve: AutoFiniteDiff # finite-difference Jacobian for the stiff solver
using LinearAlgebra: norm
using Printf

# ---- boundary conditions: fill the ghost layer from the interior ------------
function apply_bcs!(u, N, bc)
    if bc === :dirichlet # ψ = 0 on the boundary
        u[0, :] .= 0;
        u[N+1, :] .= 0
        u[:, 0] .= 0;
        u[:, N+1] .= 0
    elseif bc === :neumann # zero normal gradient: ghost = nearest interior
        u[0, :] .= u[1, :]; u[N+1, :] .= u[N, :]
        u[:, 0] .= u[:, 1]; u[:, N+1] .= u[:, N]
    else
        error("unknown bc $bc")
    end
    return nothing
end

# ---- RHS: imaginary-time GPE gradient flow ----------------------------------
function gpe_imag!(du, u, p, t)
    (; N, dx, g, μ, V, bc) = p
    if !(u isa OffsetArray)
        u = OffsetArray(u, 0:N+1, 0:N+1)
        du = OffsetArray(du, 0:N+1, 0:N+1)
    end
    apply_bcs!(u, N, bc)
    @inbounds for j in 1:N, i in 1:N
        Δψ = (u[i-1, j] + u[i+1, j] + u[i, j+1] + u[i, j-1] - 4u[i, j]) / dx^2
        du[i, j] = 0.5Δψ - V[i, j] * u[i, j] - g * u[i, j]^3 + μ * u[i, j]
    end
    # du at the ghost cells is deliberately left untouched (exactly as f_steady! does).
    return nothing
end

# ghost-padded state from an interior field: indices 0:N+1, interior 1:N
function padded(interior, N)
    u = OffsetArray(zeros(N + 2, N + 2), 0:N+1, 0:N+1)
    u[1:N, 1:N] .= interior
    return u
end

# ---- problem setup (no global constants — everything goes into p) -----------
function build_problem(bc; N = 40, L = 12.0, g = 1.0, μ = 1.0, τ_F = 100.0)
    dx = L / (N - 1)
    xs = range(-L/2, L/2, length = N) # interior coordinates
    V = [0.5 * (x^2 + y^2) for x in xs, y in xs] # harmonic trap on the interior
    u0 = padded([exp(-(x^2 + y^2) / 2) for x in xs, y in xs], N) # Gaussian seed
    p = (; N, L, dx, g, μ, V, bc)
    return ODEProblem(gpe_imag!, u0, (0.0, τ_F), p), p
end

# ---- run one solver, report what happens ------------------------------------
function try_solver(label, alg, prob, N)
    println(label); flush(stdout)
    try
        sol = solve(prob, alg; save_everystep = false, abstol = 1e-6, reltol = 1e-6)
        uf = sol.u[end]
        @printf("ran retcode=%-9s ‖ψ_final‖=%.4f\n", sol.retcode, norm(uf[1:N, 1:N]))
    catch e
        # msg = replace(sprint(showerror, e), '\n' => " ")
        # @printf("THREW %s\n%50s%s\n", nameof(typeof(e)), "", first(msg, 200))
        println(sprint(showerror, e))
    end
    flush(stdout)
end

function run_all(bc)
    prob, p = build_problem(bc)
    @printf("\n========== bc = :%s (N=%d, state size %d = (N+2)² ghost-padded, μ=%.1f) ==========\n",
            bc, p.N, (p.N + 2)^2, p.μ)
    println("── implicit (stiff), DEFAULT autodiff = ForwardDiff ──")
    try_solver("Rosenbrock23() [ForwardDiff]", Rosenbrock23(), prob, p.N)
    println("\n── implicit (stiff), autodiff = AutoFiniteDiff() ──")
    try_solver("Rosenbrock23(autodiff = AutoFiniteDiff())", Rosenbrock23(autodiff = AutoFiniteDiff()), prob, p.N)
    println("\n── explicit (no Jacobian) ──")
    try_solver("Tsit5()", Tsit5(), prob, p.N)
end

run_all(:dirichlet)
run_all(:neumann)

```

When using the implicit method with automatic differentiation, I get the error “First call to automatic differentiation for the Jacobian  
failed. This means that the user `f` function is not compatible  
with automatic differentiation …”. Using the finite differences backend gives a warning “dt was forced below floating point epsilon” in case of dirichlet boundary conditions but seems fine here for neumann boundary conditions. The Tsit5() solver always works.

---

<div class="post-metadata">

**Author:** ![lukem12345](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lukem12345/32/221287_2.png) [@lukem12345](https://discourse.julialang.org/u/lukem12345)\
**Post date:** [June 11, 2026, 5:08pm UTC](https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561/2 "2026-06-11T17:08:22Z")

</div>

Setting all of `du` explicitly to 0 right before boundary conditions are applied, i.e., `du .= 0.0`, in `gpe_imag!` seems to resolve the immediate issue and allows `Rosenbrock23(autodiff = AutoFiniteDiff())` to execute:

```julia-auto
── implicit (stiff), autodiff = AutoFiniteDiff() ──
Rosenbrock23(autodiff = AutoFiniteDiff())
ran retcode=Success ‖ψ_final‖=0.7495

── explicit (no Jacobian) ──
Tsit5()
ran retcode=Success ‖ψ_final‖=0.7496

```

I don’t use `OffsetArrays` myself so I can’t speak for whether this unblocks whatever the root issue is, but maybe this helps you get unstuck.

(I see you have a comment `# du at the ghost cells is deliberately left untouched (exactly as f_steady! does). ` but I don’t know what `f_steady!` is.)

---

<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:** [June 11, 2026, 7:42pm UTC](https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561/3 "2026-06-11T19:42:48Z")

</div>

> [@Lucas](#):
>
> This solver works fine with an explicit method like Tsit5(), but doesn’t with implicit solvers. One issue is that I used OffsetArrays, which the standard automatic differentiation backend doesn’t like. But even using finite differences, it is not working; I suspect because the ghost cells make the Jacobian singular.

It’s a DAE. Do not modify `u`, that is bad and can cause a lot of weird stuff with adaptivity. Just define it with algebraic constraints, because that’s all that BCs are anyways.

> **[Differential Algebraic Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/)**
>
> Documentation for DifferentialEquations.jl.

So just add a singular mass matrix and it’ll handle this at high order a lot more effectively.

---

<div class="post-metadata">

**Author:** ![Lucas](https://avatars.discourse-cdn.com/v4/letter/l/ba8739/32.png) [@Lucas](https://discourse.julialang.org/u/Lucas)\
**Post date:** [June 12, 2026, 10:08am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-with-ghost-cells/137561/4 "2026-06-12T10:08:54Z")

</div>

Thanks for the suggestion Chris, I re-implemented it using a mass matrix. (See the code below.) The explicit solver is still considerably faster, as is not unexpected for this problem. Is there a more “correct” way to implement it for this case? If I understand correctly, explicit solvers are not compatible with singular mass matrices.

```julia
using DifferentialEquations
using LinearAlgebra: norm, Diagonal
using SparseArrays: nnz
import Symbolics
using Printf

# ---- RHS: interior = GPE gradient flow, ghosts = BC residual ----------------
function gpe_mm!(du, u, p, t)
    (; N, dx, g, μ, V, bc) = p
    @inbounds for j in 2:N+1, i in 2:N+1
        Δψ = (u[i-1, j] + u[i+1, j] + u[i, j+1] + u[i, j-1] - 4u[i, j]) / dx^2
        du[i, j] = 0.5Δψ - V[i-1, j-1] * u[i, j] - g * u[i, j]^3 + μ * u[i, j]
    end
    apply_boundary!(du, u, N, bc)
    return nothing
end

# ---- ghost-cell constraints: 0 = du[ghost] = BC residual (one loop per edge) ----
function apply_boundary!(du, u, N, bc)
    if bc === :dirichlet # ψ = 0 on the boundary
        @inbounds for k in 1:N+2; du[k, begin] = u[k, begin]; end
        @inbounds for k in 1:N+2; du[k, end] = u[k, end]; end
        @inbounds for k in 2:N+1; du[begin, k] = u[begin, k]; end
        @inbounds for k in 2:N+1; du[end, k] = u[end, k]; end
    elseif bc === :neumann # zero normal gradient: ghost = nearest interior
        @inbounds for k in 2:N+1; du[k, begin] = u[k, begin] - u[k, begin+1]; end
        @inbounds for k in 2:N+1; du[k, end] = u[k, end] - u[k, end-1]; end
        @inbounds for k in 2:N+1; du[begin, k] = u[begin, k] - u[begin+1, k]; end
        @inbounds for k in 2:N+1; du[end, k] = u[end, k] - u[end-1, k]; end
        @inbounds begin # corners pinned to 0
            du[begin, begin] = u[begin, begin]; du[begin, end] = u[begin, end]
            du[end, begin] = u[end, begin]; du[end, end] = u[end, end]
        end
    else
        error("unknown bc $bc")
    end
    return nothing
end

# ---- Jacobian sparsity pattern, detected automatically with Symbolics -------
function jac_sparsity(p)
    N = p.N
    rhs!(du, u) = gpe_mm!(du, u, p, 0.0)
    S = Symbolics.jacobian_sparsity(rhs!, zeros(N + 2, N + 2), zeros(N + 2, N + 2))
    return Float64.(S)
end

# ---- problem setup (no global constants — everything goes into p) -----------
function build_problem(bc; N = 160, L = 32.0, g = 1.0, μ = 1.0, τ_F = 1e4)
    dx = L / (N - 1)
    xs = range(-L/2, L/2, length = N) # interior coordinates
    V = [0.5 * (x^2 + y^2) for x in xs, y in xs] # harmonic trap on the interior

    u0 = zeros(N + 2, N + 2)
    u0[2:N+1, 2:N+1] .= [exp(-(x^2 + y^2) / 2) for x in xs, y in xs] # Gaussian seed
    bc === :neumann && _fill_neumann_ghosts!(u0, N) # consistent IC on the ghosts

    mask = falses(N + 2, N + 2) # singular mass matrix: 1 interior, 0 ghosts
    mask[2:N+1, 2:N+1] .= true
    M = Diagonal(Float64.(vec(mask)))

    p = (; N, L, dx, g, μ, V, bc)
    f = ODEFunction(gpe_mm!, mass_matrix = M, jac_prototype = jac_sparsity(p))
    return ODEProblem(f, u0, (0.0, τ_F), p), p
end

function _fill_neumann_ghosts!(u, N)
    lo, hi = 1, N + 2
    u[lo, 2:N+1] .= u[lo+1, 2:N+1]; u[hi, 2:N+1] .= u[hi-1, 2:N+1]
    u[2:N+1, lo] .= u[2:N+1, lo+1]; u[2:N+1, hi] .= u[2:N+1, hi-1]
    return nothing
end

# ---- warm up (compile out of the timing), then @time the solve --------------
function try_solver(label, prob, N)
    println(label); flush(stdout)
    solve(prob, Rodas5P(); save_everystep = false, abstol = 1e-6, reltol = 1e-6) # warm-up
    sol = @time solve(prob, Rodas5P(); save_everystep = false, abstol = 1e-6, reltol = 1e-6)
    @printf("ran retcode=%-9s ‖ψ_final‖=%.4f\n", sol.retcode, norm(sol.u[end][2:N+1, 2:N+1]))
    flush(stdout)
end

function run_all(bc)
    prob, p = build_problem(bc)
    @printf("\n========== bc = :%s (N=%d, state size %d = (N+2)² with mass-matrix BCs, μ=%.1f) ==========\n",
            bc, p.N, (p.N + 2)^2, p.μ)
    S = prob.f.jac_prototype
    @printf("Jacobian sparsity: %d×%d, %d structural nonzeros (%.3f%% dense)\n",
            size(S, 1), size(S, 2), nnz(S), 100 * nnz(S) / length(S))
    try_solver("Rodas5P()", prob, p.N)
end

run_all(:dirichlet)
run_all(:neumann)

```