Discretization of semi-coupled wave and diffusion equations fails for step size of pi/(25n)

[NonDairyNeutrino]

Intro

I'm using a semi-coupled system of wave and diffusion equations to simulate a string whose mass/speed of sound diffuses (as a toy to familiarize myself with MethodOfLines.jl for my master's thesis in high-performance computational physics). The initial profiles are both normal distributions, and as such I choose the spatial domain to be of length 4*pi. But I encounter an error in the discretization.

Given a step size dx = pi / n, a warning is issued

Warning: dx for x does not divide domain exactly, adding grid point at x = 6.283185307179586).

followed by a BoundsError such as

ERROR: LoadError: BoundsError: attempt to access 401-element Vector{SymbolicUtils.BasicSymbolic{Real}} at index [CartesianIndex{1}[CartesianIndex(0,), CartesianIndex(1,)]]

any time n is a multiple of 25, as far as I have tested. Additionally, no flags are raised when n = 5. I tested for multiple different dependencies on n, e.g. multiples of 5 not equal to multiples of 25, +-1 around multiples of 25, and all have worked just fine.

I've tried using the VSCode debugger, but in this particular instance, it quickly gets over my head.

Code (does not work)

It appears the BoundsError is raised from within discretize. As such I have omitted the later code to actually get and visualize the solution.

Note: dx = pi / 25.

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets

# independent and dependent variables
@parameters t x        # t - time, x - space
@variables u(..) c(..) # u - wave amplitude, c - speed of sound

# differential operators
Dt  = Differential(t)   # first  order time derivative operator
Dtt = Differential(t)^2 # second order time derivative operator
Dx  = Differential(x)   # first  order spatial derivative operator
Dxx = Differential(x)^2 # second order spatial derivative operator

# governing equations
eqwave  = Dtt(u(t, x)) ~ c(t, x)^2 * Dxx(u(t, x)) # wave equation with variable speed
eqspeed =  Dt(c(t, x)) ~ Dxx(c(t, x))             # diffusion equation for mass/speed

# boundary conditions
bcs = [
    # initial conditions
       u(0, x)  ~ exp(-x^2),   # initial spatial wave profile
    Dt(u(0, x)) ~ 0.,          # initial wave velocity
       c(0, x)  ~ exp(-x^2),   # initial spatial speed profile
    # boundary conditions (optional)
    Dx(u(t, -2pi)) ~ 0., # left side is rigid
       u(t,  2pi)  ~ 0.  # right side held at zero
]

# Space and time domains
domains = [
    t ∈ Interval(0.0, 20.0),      # arbitrary time domain
    x ∈ Interval(-2 * pi, 2 * pi) # arbitrary spatial domain
]

# PDE system
@named pdesys = PDESystem([equ, eqc], bcs, domains, [t, x], [u(t, x), c(t, x)])

# Method of lines discretization
dx = pi / 25. # doesn't work when denominator is a multiple of 25
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)

Versions

Packages

DomainSets has been held back from updating due to compatibility constraints.

[0c46a032] DifferentialEquations v7.11.0 ⌅ [5b8099bc] DomainSets v0.6.7 [7a1cc6ca] FFTW v1.7.1 [94925ecb] MethodOfLines v0.10.1 [961ee093] ModelingToolkit v8.72.2 [1dea7af3] OrdinaryDiffEq v6.59.0

Julia

1.9.3+0.x64.w64.mingw32

[xtalax]

I see the problem, you also need to supply boundary conditions for c, this needs a better error message though