Unhelpful error message: "ERROR: MethodError: reducing over an empty collection is not allowed; consider supplying `init` to the reducer"

[gkantsidis]

The following is a small version of a larger problem:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables Nd(..)

Dt = Differential(t)

eq1 = Dt(Nd(t, x)) ~ 24.3 * (1.0 - Nd(t, x))

eq = [
    eq1,
]

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

boundary = [
    Nd(0, x) ~ 0.9,
]

@named pde_system = PDESystem(eq, boundary, domains, [t, x], [Nd(t, x)])

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

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

prob = discretize(pde_system, discretization)

The last call fails with the error: ERROR: MethodError: reducing over an empty collection is not allowed; consider supplying init to the reducer

I would assume that there is an error in setting up the system (apologies, I'm very new to this space), but the error message is not helpful in tracing down the problem.

(I'm using Julia v1.8.2, MethodLines v0.5.0)

[xtalax]

I have seen this error before but no idea why your setup might be causing it in this case. As a hunch, try expanding out your brackets as a workaround, I will take a closer look at this when I get time.

[gkantsidis]

Thanks for taking a look at it. I am not sure I understand what you mean by "brackets". If you mean to replace [x=>dx] by x=>dx, then I get ERROR: MethodError: no method matching keys(::Pair{Num, Float64}).

I have played around with the ranges of t and x, and with the value of the approximation order, but I still get the same error.

On the other hand, if I discretize time:

discretization = MOLFiniteDifference([x=>dx, t=>dt], nothing, approx_order=order, grid_align=center_align)

(instead of discretization = MOLFiniteDifference([x=>dx], t, approx_order=order, grid_align=center_align)) then it seems to work ok, i.e. the last statement above works, and then I can call the solve method.

[xtalax]

I mean try eq1 = Dt(Nd(t, x)) ~ 24.3 - 24.3*Nd(t, x)

[gkantsidis]

Same problem with eq1 = Dt(Nd(t, x)) ~ 24.3 * Nd(t, x) as before.

[BernhardAhrens]

I believe discretize expects not only an initial condition but also a boundary condition, although you do not use a space derivative. If you specify a zero-flux boundary at the top (I think that's what I have done), your code runs:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables Nd(..)

Dt = Differential(t)
Dx = Differential(x)

eq1 = Dt(Nd(t, x)) ~ 24.3 * (1.0 - Nd(t, x))

eq = [
    eq1
]

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

boundary = [
    Nd(0, x) ~ 0.9
    Nd(t,x_min) - Dx(Nd(t, x_min)) ~ 0.0 # zero flux boundary?
]

@named pde_system = PDESystem(eq, boundary, domains, [t, x], [Nd(t, x)])

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

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

prob = discretize(pde_system, discretization)