System of PDEs `solve` error

[henry2004y]

Describe the bug 🐞

For a PDE system with 3 unknowns as a function of 1D space x and time t, instability is detected while solving. The equations to be solved are the simplest case for Langmuir waves in plasma physics.

Expected behavior

Here is a simple finite difference solver for the problem:

$$ \begin{aligned} \frac{\partial n}{\partial t} &= -n_0\frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial t} &= \frac{e}{m}\frac{\partial \phi}{\partial x} \\ \frac{\partial^2 \phi}{\partial x^2} &= \frac{e}{\epsilon_0}n \end{aligned} $$

in a domain $x\in [0, 2\pi], t\in [0, 1]$ with a periodic BC and initial conditions

$$ \begin{aligned} n &= 0.02 \cos(x) \\ u &= 0 \\ \phi &= 0 \end{aligned} $$

I expect MethodOfLines.jl should produce similar results.

function solve_poisson(source, dx)
   N = length(source)

   # Construct the finite difference matrix (periodic BC)
   A = zeros(N, N)
   for j in axes(A, 2), i in axes(A, 1)
      if i == j
         A[i,j] = -2 / dx^2
      elseif abs(i - j) == 1
         A[i,j] = 1 / dx^2
      elseif (i == 1 && j == N) || (i == N && j == 1)
         A[i,j] = 1 / dx^2
      end
   end

   # Solve the linear system
   potential = A \ source

   return potential
end

# Constants
n0 = 1.0         # Background density
e = 1.0          # Electron charge 
m = 1.0          # Electron mass
ε₀ = 1.0         # Permittivity of free space

# Discretization
L = 2π           # Domain length
N = 100          # Number of grid points
dx = L / N
x = range(0, L, length=N)

tspan = (0.0, 1.0)
dt = 0.01
trange = range(tspan[1], tspan[2], step=dt)

# Initial conditions
n = @. 0.02*cos(x)  # Density perturbation
u = zeros(N)        # Zero initial velocity
E = zeros(N)        # Zero initial electric field
source = similar(n) # Source term for Poisson's equation

# Time stepping loop
for t in trange
   # Update electric field (Poisson's equation)
   @. source = e/ε₀ * n
   potential = solve_poisson(source, dx)

   for i in eachindex(E)
      if i == 1
         E[i] = (potential[end] - potential[2]) / 2dx
      elseif i == N
         E[i] = (potential[end-1] - potential[1]) / 2dx
      else
         E[i] = (potential[i-1] - potential[i+1]) / 2dx
      end
   end

   # Update velocity (momentum equation)
   @. u += -dt * e/m * E

   # Update density (continuity equation)
   for i in eachindex(n)
      if i == 1
         n[i] += - dt * n0 * (u[2] - u[end]) / 2dx
      elseif i == N
         n[i] += - dt * n0 * (u[1] - u[end-1]) / 2dx
      else
         n[i] += - dt * n0 * (u[i+1] - u[i-1]) / 2dx
      end
   end
end

# Visualization
#=
using PyPlot

plot(x, n, label="Density")
plot(x, u, label="Velocity")
plot(x, E, label="Electric field") 
xlabel("x")
ylabel("n, E")
legend()
=#

Minimal Reproducible Example 👇

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x t
@variables n(..) u(..) ϕ(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

t_min = 0.0
t_max = 1.0
x_min = 0.0
x_max = 2pi

n0 = 1.0
e = 1.0
m = 1.0
ε₀ = 1.0
	
eq = [
   Dt(n(x,t)) ~ - n0 * Dx(u(x,t)),
   Dt(u(x,t)) ~ e / m * Dx(ϕ(x,t)),
   Dxx(ϕ(x,t)) ~ e / ε₀ * n(x,t)
]

domains = [
   x ∈ Interval(x_min, x_max),
   t ∈ Interval(t_min, t_max)
]

bcs = [
   n(x, t_min) ~ 0.02*cos(x),
   n(x_min, t) ~ n(x_max, t),
       
   u(x, t_min) ~ 0,
   u(x_min, t) ~ u(x_max, t),

   ϕ(x, t_min) ~ 0,
   ϕ(x_min, t) ~ ϕ(x_max, t)
]

@named pdesys = PDESystem(eq, bcs, domains, [x, t], [n(x,t), u(x,t), ϕ(x,t)])

N = 16
dx = (x_max - x_min) / N

discretization = MOLFiniteDifference([x=>dx], t, approx_order=2)

@time prob = discretize(pdesys, discretization)

@time sol = solve(prob, ROS3P(), saveat=0.1)

Error & Stacktrace ⚠️

┌ Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.
└ @ MethodOfLines C:\Users\hyzho\.julia\packages\MethodOfLines\xyn4D\src\system_parsing\pde_system_transformation.jl:43
 13.934843 seconds (20.65 M allocations: 1.293 GiB, 5.27% gc time, 92.48% compilation time: <1% of which was recompilation)
┌ Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\hyzho\.julia\packages\SciMLBase\NkxGm\src\integrator_interface.jl:623
┌ Warning: Solution has length 1 in dimension x. Interpolation will not be possible for variable u(x, t). Solution return code is Unstable.
...

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `D:\Computer\MOL\Project.toml`
⌅ [5b8099bc] DomainSets v0.6.7
  [94925ecb] MethodOfLines v0.11.0
  [961ee093] ModelingToolkit v9.12.0
  [1dea7af3] OrdinaryDiffEq v6.74.1

• Output of versioninfo()

Julia Version 1.10.2
Commit bd47eca2c8 (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 12 × Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)