Avoid auxiliary variables in the definition of boundary conditions?

[BernhardAhrens]

Consider the following simple PDE where I wrapped the convective Flux into an auxiliary variable:

# 1D PDE and boundary conditions
eq = [Dt(u(t, x)) ~ -Dx(Flux(t,x)),
Flux(t,x) ~ v / rho * u(t, x)
]

The system should describe the sedimentation and displacement of material with a sedimentation influx v and density rho, which translates v into an advection velocity (v / rho).

When I use the Flux(t,0) auxiliary variable in the boundary conditions, I get a concentration u higher than the theoretical maximum rho.

bcs = [u(0, x) ~ 0,
Flux(t,0) ~ v]

Using the exact definition of Flux(t,x) instead in the boundary conditions seems to give the correct solution.

bcs = [u(0, x) ~ 0,
v/rho * u(t, 0) ~ v]

My interpretation from this is that auxiliary variables cannot be used in the boundary conditions. It could be something for documentation if I haven't missed it.

Here is the complete simple example:

using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots
using OrdinaryDiffEq

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

@parameters v
@parameters rho
ps = [v => 0.5, rho => 145.0]

# 1D PDE and boundary conditions
eq = [Dt(u(t, x)) ~ -Dx(Flux(t,x)),
Flux(t,x) ~ v / rho * u(t, x)
]

# OPTION #1
# goes above rho
bcs = [u(0, x) ~ 0,
Flux(t,0.0) ~ v]

# OPTION #2
# does not go above rho, looks okay
#bcs = [u(0, x) ~ 0,
#v/rho * u(t, 0) ~ v]

# Space and time domains
domains = [t ∈ Interval(0.0, 200.0),
    x ∈ Interval(0.0, 2.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), Flux(t,x)], ps)

# Method of lines discretization
dx = 0.1
order = 1
discretization = MOLFiniteDifference([x => dx], t)
# explicitly specify upwind order

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

# Solve ODE problem
sol = solve(prob, Tsit5(), saveat=10.0)

p1 = plot(xlab="u(t,z)", ylab="depth", legendtitle="Time sol.t")
for i in 1:length(sol.t)
    p1 = plot!(collect(sol[u(t, x)][i, 1:end]), -collect(sol[x]), legend=:outertopright, label=round(sol.t[i]; sigdigits=2))
end
p1

[xtalax]

This is an unknown problem, thank you for bringing this to my attention. I'll add a note to the docs.

Can you try without aux variables? System transformation should handle this but we don't have too many test cases.

[BernhardAhrens]

Removing aux variables in eq and bcs produces the same results for

eq = [Dt(u(t, x)) ~ -Dx(Flux(t,x)), Flux(t,x) ~ v / rho * u(t, x)]

and

eq2 = [Dt(u(t, x)) ~ -Dx(v / rho * u(t, x))]

However, with the aux variable in bcs, you get the discrepancy that I mentioned in the first post: bcs = [u(0, x) ~ 0, Flux(t,0) ~ v] and bcs = [u(0, x) ~ 0, v/rho * u(t, 0) ~ v] for solving eq give different results. bcs = [u(0, x) ~ 0, v/rho * u(t, 0) ~ v] gives the correct results.

It seems like you can use aux variables in the equations but not the boundary conditions.

[Qfl3x]

If Flux(t,x) is in @named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), Flux(t,x)], ps), shouldn't MoL treat it as a dependent variable not any different from u(t,x)?

[xtalax]

Yes, to my knowledge this should be equivalent -- this is definitely unexpected. From experience the cause is likely upstream, but this is tough to verify. How does it do without structural_simplify i.e. using symbolic_discretize rather than discretize?

[Qfl3x]

The result is a Tuple{ODESystem, Tuple{Float64,Float64}} rather than a ODEProblem; solve crashes. Tried converting but it also crashes with some SymbolicUtils error.

[xtalax]

the tuple is (sys, tspan), some systems need simplification to be compatible though

[Qfl3x]

This system appears to be incompatible, tried transforming using https://github.com/SciML/MethodOfLines.jl/blob/feac4dc0851a23a8a8213a39605c823cd828f7e7/src/MOL_discretization.jl#L149 without simplification but it didn't work. When I used structural_simplify it worked but the same problem occured.

[xtalax]

Aux variables are treated identically by the system - I will double check the boundary parsing to make sure there isn't an extra assumption about the form of the bcs that is applying that I may have missed and refer this upstream if there's nothing