discretize only generates N-1 variables for two-sided Dirichlet BC of Laplace equation

[rschiemann1]

The following code (hopes to) implement a discretization of the Laplace equation with Dirichlet BCs on the left and right end of the Interval into 5 grid points, but discretize constructs a problem that only represents four points. Is that by design, due to Dirichlet BC handling when grid_align = center_align?

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x t
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

∇²(u) = Dxx(u)

x_min = 0.0
x_max = 1.0

t_min = 0.0
t_max = 11.5

α = 10.

eq = [Dt(u(x,t)) ~  α*∇²(u(x,t))]

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


bcs = [u(x,0) ~ 20,
       Dx(u(0,t)) ~ 0,
       Dx(u(x_max,t)) ~ 0] 

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

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

# despite N = 5, prob only has 4 states.
prob = discretize(pdesys,discretization)

[xtalax]

That would be due to structural_simplify simplifying away points due to the specified Neumann BCs, please use sol[u] to retrieve all points including those that were simplified out.

[rschiemann1]

thanks for the hint. This does not seem to work:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x t
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

∇²(u) = Dxx(u)

x_min = 0.0
x_max = 1.0

t_min = 0.0
t_max = 11.5

α = 10.

eq = [Dt(u(x,t)) ~  α*∇²(u(x,t))]

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


bcs = [u(x,0) ~ 20,
       Dx(u(0,t)) ~ 0,
       Dx(u(x_max,t)) ~ 0] 

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

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

# despite N = 5, prob only has 4 states.
prob = discretize(pdesys,discretization)

sol = solve(prob, TRBDF2(), saveat=1)

# ERROR: ArgumentError: u⋆ is either an observed nor a state variable.
sol[u]

@variables u(x,t)

# ERROR: ArgumentError: u(x, t) is either an observed nor a state variable.
sol[u]

[xtalax]

@YingboMa I'm sure that this should work, we're doing it throughout the tests. All state/observed variables look like u[i](t). What is going wrong here?