RuleRewriteError: Chebyshev grid + UpwindScheme fails for parameterized advection terms

[ChrisRackauckas-Claude]

Description

When using a Chebyshev (non-uniform) grid with UpwindScheme(), discretize fails with a SymbolicUtils.RuleRewriteError for PDESystems that have parameterized advection terms like v * Dx(u).

Uniform grids with UpwindScheme() work fine for the same systems.

Error

SymbolicUtils.RuleRewriteError: Failed to apply rule
  ~(~a) * ($(Differential(x) ^ d))(u) * ~(~b) => upwind_difference(...)
on expression v*Differential(x, 1)(u(t, x))

MWE

using MethodOfLines, DomainSets, OrdinaryDiffEq, ModelingToolkit

@parameters x t v
@variables u(..)
Dx = Differential(x)
Dt = Differential(t)

eq = Dt(u(t, x)) ~ -v * Dx(u(t, x))
bcs = [u(0, x) ~ sin(2π * x), u(t, 0) ~ 0.0]
domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)], [v => 0.5])

# Uniform grid — works fine
disc_uniform = MOLFiniteDifference([x => 100], t, advection_scheme = UpwindScheme())
prob = discretize(pdesys, disc_uniform)  # OK

# Chebyshev grid — fails
gridnu = chebyspace(100, domains[2])
disc_cheby = MOLFiniteDifference([gridnu], t, advection_scheme = UpwindScheme())
prob = discretize(pdesys, disc_cheby)  # ERROR: RuleRewriteError

Parameter-free systems (e.g. -u(t,x) * Dx(u(t,x))) work with both grid types.

Versions

• MethodOfLines v0.11.11
• ModelingToolkit v11.17.0
• PDEBase v0.1.22
• SymbolicUtils v4.20.2

Impact

Affects all PDESystemLibrary systems with parameters when using Chebyshev discretizations in benchmarks.

[utkuyilmaz1903]

For non-uniform grids, the backward difference DerivativeOperator incorrectly gets D.offside == 0. In _upwind_difference evaluating weights = D.stencil_coefs[II[j] - D.offside] shifts the index by one. At the rightmost boundary N, this requests the Nth element from an N-1 length array and causes the crash.
I locally patched this by enforcing true_offside = D.stencil_length - 1. This eliminates the error and works perfectly for 1st-order upwinding. However, it is mathematically flawed for 2nd-order or higher on non-uniform grids because it fetches the coefficient for the wrong local evaluation point on the curve. I think the true solution is not to shift indexes locally. The upstream generation of D.offside and D.stencil_coefs needs to be fixed for non-uniform operators.

Thoughts on this @xtalax?

[utkuyilmaz1903]

@ChrisRackauckas Taking a closer look at this after my previous comment regarding the D.offside indexing, I've realized that the root cause here isn't actually an operator generation bug. The MWE is mathematically ill-posed.

Because v is defined as a parameter (@parameters v), its sign isn't known during the discretize phase (compile-time). To handle this dynamically at runtime, the codegen must generate an IfElse branch for both v > 0 (looking left) and v < 0 (looking right).

This means the system mathematically requires both boundary conditions to compile safely. However, the MWE only provides the left BC (u(t, 0) ~ 0.0). When the codegen attempts to build the v < 0 branch at the right boundary, it has no BC to fall back on. It blindly indexes out of bounds, which triggers a BoundsError that SymbolicUtils swallows and masks as a RuleRewriteError.

• Why Uniform grids "work": They blindly calculate x + dx to create a ghost node, bypassing the index crash (though producing physically questionable results if v turns negative).
• Why Chebyshev grids fail: They lack a scalar dx, so they hit the hard index limit.

Instead of patching the operator indexing, the true solution is that the user must provide the missing right BC (e.g., u(t, 1) ~ 0.0) so both parameter branches can compile.

I am currently working on PR #543 to catch these exact types of missing boundary conditions early in the flow. It will throw a clear, graceful ArgumentError before we ever reach the codegen ruleset, preventing these confusing mask errors entirely.