MethodOfLines.jl MethodOfLines complains about my model...

Here is my model:

This 1D model describes the flow of ideal gas in a pipe. I assume that inlet pressure $p_0$ is a known function of time, and that the outlet mass flow rate $\dot{m}_L$ is a known function of time.

I get a few warnings after the statement:

@named pipe = PDESystem(eqs,bcs,doms,[t,x],[p(t,x), ṁ(t,x)], 
                    [RdM=>RdM0, T=>T0, A_p=>Ap0, g=>g0, H=>H0, L=>L0, ℘=>℘0, f=>f0])

but the code still runs. However, when I then do:

prob = discretize(pipe,disc);

I get more warnings, and eventually the running crashes.

Question: Is there a problem with the term:

?? Is only the differentiation of a simple function allowed?

If so, how do I get around it?

Differentiating out, or
Adding an algebraic definition of the argument of the differentiation?

If you guys need it, I can provide the entire code.

Mar 08 '23 13:03 B-LIE

In case you are interested in my code [note! I have successfully solved at least 3 PDEs using MethodsOfLine.jl, so I have a little experience with it]:

# PACKAGES
using ModelingToolkit
using DifferentialEquations
using Plots, LaTeXStrings
using MethodOfLines, DomainSets

# MODEL
parameters RdM T A_p g H L ℘ f 
@variables t x p(..) ṁ(..) # p_0(t) ṁ_L(t)
Dt = Differential(t)
Dx = Differential(x)
#
eqs = [Dt(p(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)),
        Dt(ṁ(t,x)) ~ -(RdM*T)/(A_p)*Dx((ṁ(t,x))^2/p(t,x)) - A_p*Dx(p(t,x)) +
                (g*H)/L*A_p/(RdM*T)*p(t,x) - RdM*T*℘/A_p^2/2*f*(ṁ(t,x))^2/p(t,x)]

#=  In case term needs to be expanded out:
eqs = [Dt(p(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)),
        Dt(ṁ(t,x)) ~ -(RdM*T)/(A_p)*( 2*ṁ(t,x)/p(t,x)*Dx(ṁ(t,x)) - (ṁ(t,x))^2/p(t,x)*Dx(p(t,x)) ) - A_p*Dx(p(t,x)) +
                (g*H)/L*A_p/(RdM*T)*p(t,x) - RdM*T*℘/A_p^2/2*f*(ṁ(t,x))^2/p(t,x)]
=#

# NUMERIC VALUES
# From A. Carlson and B. Lie, 1992
D0 = 1.422 # m 
Ap0 =pi*D0^2/4
℘0 = pi*D0
L0 = 122e3 # m 
RdM0 = 493 # J/(kg.K)
T0 = 12 + 273.15 # K
f0 = 0.28e-3 # friction factor
g0 = 9.81
H0 = 0
p0 = 8.4e6 # Pa
ṁ0 = 1.5e6/3600 # kg/s

# BOUNDARY FUNCTIONS
@register_symbolic p_0(t)
@register_symbolic ṁ_L(t)

p_0(t) = p0
ṁ_L(t) = ṁ0

# DOMAINS
L = L0
doms = [t ∈ Interval(0.0, 30*3600), x ∈ Interval(0, L)]

# (INITIAL &) BOUNDARY CONDITIONS
bcs = [ p(0,x) ~ p0,
        ṁ(0,x) ~ ṁ0,
        p(t,0) ~ p_0(t),
        ṁ(t,L) ~ ṁ_L(t)]

# SYMBOLIC MODEL
@named pipe = PDESystem(eqs,bcs,doms,[t,x],[p(t,x), ṁ(t,x)], 
                    [RdM=>RdM0, T=>T0, A_p=>Ap0, g=>g0, H=>H0, L=>L0, ℘=>℘0, f=>f0])

# The above leads to warnings, but nothing that stops the execution

# DISCRETIZATION INTO NUMERIC MODEL
Nmol = 10
order = 2 
disc = MOLFiniteDifference([x=>Nmol], t, approx_order=order)

prob = discretize(pipe,disc)

# The above `discretize` function runs for 1 minute ++ on my 12th gen intel laptop, then crashes

# SOLVING MODEL AND EXTRACTING RESULTS
sol = solve(prob)

solp = sol[p(t,x)]
solṁ = sol[ṁ(t,x)]
solx = sol[x]
solt = sol[t]

# Code crashes before this headline

Mar 08 '23 14:03 B-LIE

Hm... I also tried with the following re-formulation (the model is the simplified Euler equations, i.e., without energy balance):

This model avoids division by some variables, at least in the PDEs themselves (but shifting those operations to the algebraic constraints). And also makes sure that the differentiation operators are used on simple functions (instead of expressions). In this update, I have corrected a tiny non-physicality of the previous model to ensure that friction opposes flow also when the flow direction changes.

Anyways, I get the same error message with these equations.

Mar 09 '23 09:03 B-LIE

I think that the Dx((ṁ(t,x))^2 is causing problems with the rule application, what error are you getting? Please post full output including warnings and stacktrace.

Question: Is there a problem with the term:

Can you write that as Dx(m(t,x)^2/p(t,x))? This system should be compatible with automatic transformation, so you could try writing it like that. The registered BCs may cause an issue there though, can you try without registering, replacing any control flow with ifelse. Also, make sure you are running the latest version.

Mar 09 '23 18:03 xtalax

Can you write that as Dx(m(t,x)^2/p(t,x))?

Yes, I tried that prior to writing Dx((m(t,x))^2/p(t,x)).

I have tried to avoid registered functions in BCs... Here is the latest attempt:

@parameters RdM T A_p g H L ℘ f 
@variables t x p(..) ṁ(..) # p_0(t) ṁ_L(t)
Dt = Differential(t)
Dx = Differential(x)
#
eqs = [Dt(p(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)),
        Dt(ṁ(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)^2/p(t,x)) - A_p*Dx(p(t,x)) +
                (g*H)/L*A_p/(RdM*T)*p(t,x) - RdM*T*℘/A_p^2/2*f*(ṁ(t,x))^2/p(t,x)]

# From A. Carlson and B. Lie, 1992
D0 = 1.422 # m 
Ap0 =pi*D0^2/4
℘0 = pi*D0
L0 = 122e3 # m 
RdM0 = 493 # J/(kg.K)
T0 = 12 + 273.15 # K
f0 = 0.28e-3 # friction factor
g0 = 9.81
H0 = 0
p0 = 8.4e6 # Pa
ṁ0 = 1.5e6/3600 # kg/s

L = L0
doms = [t ∈ Interval(0.0, 30*3600), x ∈ Interval(0, L)]

# Boundary conditions, including Initial conditions

bcs = [ p(0,x) ~ p0,
        ṁ(0,x) ~ ṁ0,
        p(t,0) ~ p0,
        ṁ(t,L) ~ ṁ0]

@named pipe = PDESystem(eqs,bcs,doms,[t,x],[p(t,x), ṁ(t,x)], 
                    [RdM=>RdM0, T=>T0, A_p=>Ap0, g=>g0, H=>H0, L=>L0, ℘=>℘0, f=>f0])

Nmol = 10

order = 2 # This may be increased to improve accuracy of some schemes

# Integers for x and y are interpreted as number of points. Use a Float to directtly specify stepsizes dx and dy.
disc = MOLFiniteDifference([x=>Nmol], t, approx_order=order)

prob = discretize(pipe,disc)

The latest line (discretize) cause a crash, even when I have removed registered functions from the BCs.

Here is the error message:

┌ Warning: Incompatible term found. Adding auxiliary equation var"⟦(ṁ(t, x)^2) / p(t, x)⟧"(t, x) ~ (ṁ(t, x)^2) / p(t, x) to the system.
└ @ MethodOfLines [C:\Users\Bernt\.julia\packages\MethodOfLines\44Kiz\src\system_parsing\pde_system_transformation.jl:214](file:///C:/Users/Bernt/.julia/packages/MethodOfLines/44Kiz/src/system_parsing/pde_system_transformation.jl:214)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
┌ Warning: : no method matching get_unit for arguments (Pair{Num, Float64},).
└ @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\validation.jl:154](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/validation.jl:154)
The system of equations is:
Equation[Differential(t)((p(t))[2]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[2] - 7.377049180327869e-5(ṁ(t))[1])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[3] - 7.377049180327869e-5(ṁ(t))[2])) / A_p) ~ 0, Differential(t)((p(t))[3]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[3] - 7.377049180327869e-5(ṁ(t))[2])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[4] - 7.377049180327869e-5(ṁ(t))[3])) / A_p) ~ 0, Differential(t)((p(t))[4]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[4] - 7.377049180327869e-5(ṁ(t))[3])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[5] - 7.377049180327869e-5(ṁ(t))[4])) / A_p) ~ 0, Differential(t)((p(t))[5]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[5] - 7.377049180327869e-5(ṁ(t))[4])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[6] - 7.377049180327869e-5(ṁ(t))[5])) / A_p) ~ 0, Differential(t)((p(t))[6]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[6] - 7.377049180327869e-5(ṁ(t))[5])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[7] - 7.377049180327869e-5(ṁ(t))[6])) / A_p) ~ 0, Differential(t)((p(t))[7]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[7] - 7.377049180327869e-5(ṁ(t))[6])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[8] - 7.377049180327869e-5(ṁ(t))[7])) / A_p) ~ 0, Differential(t)((p(t))[8]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[8] - 7.377049180327869e-5(ṁ(t))[7])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[9] - 7.377049180327869e-5(ṁ(t))[8])) / A_p) ~ 0, Differential(t)((p(t))[9]) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(ṁ(t))[9] - 7.377049180327869e-5(ṁ(t))[8])) / A_p, (RdM*T*(7.377049180327869e-5(ṁ(t))[10] - 7.377049180327869e-5(ṁ(t))[9])) / A_p) ~ 0, (-A_p*H*g*(p(t))[2]) / (L*RdM*T) + (RdM*T*f*℘*((ṁ(t))[2]^2)) / (2(A_p^2)*(p(t))[2]) + Differential(t)((ṁ(t))[2]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[2] - 7.377049180327869e-5(p(t))[1]), A_p*(7.377049180327869e-5(p(t))[3] - 7.377049180327869e-5(p(t))[2])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[2] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[1])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[3] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[2])) / A_p) ~ 0, (RdM*T*f*℘*((ṁ(t))[3]^2)) / (2(A_p^2)*(p(t))[3]) + (-A_p*H*g*(p(t))[3]) / (L*RdM*T) + Differential(t)((ṁ(t))[3]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[3] - 7.377049180327869e-5(p(t))[2]), A_p*(7.377049180327869e-5(p(t))[4] - 7.377049180327869e-5(p(t))[3])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[3] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[2])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[4] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[3])) / A_p) ~ 0, (-A_p*H*g*(p(t))[4]) / (L*RdM*T) + (RdM*T*f*℘*((ṁ(t))[4]^2)) / (2(A_p^2)*(p(t))[4]) + Differential(t)((ṁ(t))[4]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[4] - 7.377049180327869e-5(p(t))[3]), A_p*(7.377049180327869e-5(p(t))[5] - 7.377049180327869e-5(p(t))[4])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[4] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[3])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[5] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[4])) / A_p) ~ 0, (-A_p*H*g*(p(t))[5]) / (L*RdM*T) + (RdM*T*f*℘*((ṁ(t))[5]^2)) / (2(A_p^2)*(p(t))[5]) + Differential(t)((ṁ(t))[5]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[5] - 7.377049180327869e-5(p(t))[4]), A_p*(7.377049180327869e-5(p(t))[6] - 7.377049180327869e-5(p(t))[5])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[5] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[4])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[6] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[5])) / A_p) ~ 0, (RdM*T*f*℘*((ṁ(t))[6]^2)) / (2(A_p^2)*(p(t))[6]) + (-A_p*H*g*(p(t))[6]) / (L*RdM*T) + Differential(t)((ṁ(t))[6]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[6] - 7.377049180327869e-5(p(t))[5]), A_p*(7.377049180327869e-5(p(t))[7] - 7.377049180327869e-5(p(t))[6])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[6] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[5])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[7] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[6])) / A_p) ~ 0, (RdM*T*f*℘*((ṁ(t))[7]^2)) / (2(A_p^2)*(p(t))[7]) + (-A_p*H*g*(p(t))[7]) / (L*RdM*T) + Differential(t)((ṁ(t))[7]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[7] - 7.377049180327869e-5(p(t))[6]), A_p*(7.377049180327869e-5(p(t))[8] - 7.377049180327869e-5(p(t))[7])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[7] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[6])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[8] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[7])) / A_p) ~ 0, (-A_p*H*g*(p(t))[8]) / (L*RdM*T) + (RdM*T*f*℘*((ṁ(t))[8]^2)) / (2(A_p^2)*(p(t))[8]) + Differential(t)((ṁ(t))[8]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[8] - 7.377049180327869e-5(p(t))[7]), A_p*(7.377049180327869e-5(p(t))[9] - 7.377049180327869e-5(p(t))[8])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[8] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[7])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[9] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[8])) / A_p) ~ 0, (RdM*T*f*℘*((ṁ(t))[9]^2)) / (2(A_p^2)*(p(t))[9]) + (-A_p*H*g*(p(t))[9]) / (L*RdM*T) + Differential(t)((ṁ(t))[9]) + ifelse(A_p > 0, A_p*(7.377049180327869e-5(p(t))[9] - 7.377049180327869e-5(p(t))[8]), A_p*(7.377049180327869e-5(p(t))[10] - 7.377049180327869e-5(p(t))[9])) + ifelse(((RdM*T) / A_p) > 0, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[9] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[8])) / A_p, (RdM*T*(7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[10] - 7.377049180327869e-5(var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[9])) / A_p) ~ 0, (-((ṁ(t))[2]^2)) / (p(t))[2] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[2] ~ 0, (-((ṁ(t))[3]^2)) / (p(t))[3] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[3] ~ 0, (-((ṁ(t))[4]^2)) / (p(t))[4] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[4] ~ 0, (-((ṁ(t))[5]^2)) / (p(t))[5] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[5] ~ 0, (-((ṁ(t))[6]^2)) / (p(t))[6] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[6] ~ 0, (-((ṁ(t))[7]^2)) / (p(t))[7] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[7] ~ 0, (-((ṁ(t))[8]^2)) / (p(t))[8] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[8] ~ 0, (-((ṁ(t))[9]^2)) / (p(t))[9] + (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[9] ~ 0, (p(t))[1] ~ 8.4e6, (p(t))[10] ~ 5.0(p(t))[9] + 10.0(p(t))[7] + (p(t))[5] - 5.0(p(t))[6] - 10.0(p(t))[8], (ṁ(t))[10] ~ 416.6666666666667, (ṁ(t))[1] ~ 5.0(ṁ(t))[2] + 10.0(ṁ(t))[4] + (ṁ(t))[6] - 5.0(ṁ(t))[5] - 10.0(ṁ(t))[3], (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[10] ~ 173611.11111111112 / (p(t))[10], (var"var\"⟦(ṁ(t, x)^2) / p(t, x)⟧\""(t))[1] ~ 1.1904761904761904e-7((ṁ(t))[1]^2)]

Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

MethodError: no method matching hasmetadata(::Float64, ::Type{Symbolics.GetindexParent})

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Closest candidates are:
  hasmetadata(!Matched::SymbolicUtils.Symbolic, ::Any) at [C:\Users\Bernt\.julia\packages\SymbolicUtils\Vzo2s\src\types.jl:600](file:///C:/Users/Bernt/.julia/packages/SymbolicUtils/Vzo2s/src/types.jl:600)
  hasmetadata(!Matched::Complex{Num}, ::Any) at [C:\Users\Bernt\.julia\packages\Symbolics\VIBnK\src\Symbolics.jl:162](file:///C:/Users/Bernt/.julia/packages/Symbolics/VIBnK/src/Symbolics.jl:162)
  hasmetadata(!Matched::Num, ::Any) at [C:\Users\Bernt\.julia\packages\Symbolics\VIBnK\src\Symbolics.jl:162](file:///C:/Users/Bernt/.julia/packages/Symbolics/VIBnK/src/Symbolics.jl:162)

Stacktrace:
 [1] collect_var_to_name!(vars::Dict{Any, Any}, xs::Vector{Any})
   @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\utils.jl:267](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/utils.jl:267)
 [2] process_variables!(var_to_name::Dict{Any, Any}, defs::Dict{Any, Any}, vars::Vector{Any})
   @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\utils.jl:252](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/utils.jl:252)
 [3] ODESystem(deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing)
   @ ModelingToolkit [C:\Users\Bernt\.julia\packages\ModelingToolkit\3ePwL\src\systems\diffeqs\odesystem.jl:203](file:///C:/Users/Bernt/.julia/packages/ModelingToolkit/3ePwL/src/systems/diffeqs/odesystem.jl:203)
 [4] generate_system(alleqs::Vector{Any}, bceqs::Vector{Any}, ics::Vector{Any}, discvars::Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Num}}, defaults::Dict{Num, Float64}, ps::Vector{Num}, tspan::Tuple{Float64, Float64}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization})
   @ MethodOfLines [C:\Users\Bernt\.julia\packages\MethodOfLines\44Kiz\src\scalar_discretization.jl:61](file:///C:/Users/Bernt/.julia/packages/MethodOfLines/44Kiz/src/scalar_discretization.jl:61)
 [5] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ MethodOfLines [C:\Users\Bernt\.julia\packages\MethodOfLines\44Kiz\src\MOL_discretization.jl:128](file:///C:/Users/Bernt/.julia/packages/MethodOfLines/44Kiz/src/MOL_discretization.jl:128)
 [6] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ MethodOfLines [C:\Users\Bernt\.julia\packages\MethodOfLines\44Kiz\src\MOL_discretization.jl:132](file:///C:/Users/Bernt/.julia/packages/MethodOfLines/44Kiz/src/MOL_discretization.jl:132)
 [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ MethodOfLines [C:\Users\Bernt\.julia\packages\MethodOfLines\44Kiz\src\MOL_discretization.jl:131](file:///C:/Users/Bernt/.julia/packages/MethodOfLines/44Kiz/src/MOL_discretization.jl:131)
 [8] top-level scope
   @ [c:\Users\Bernt\OneDrive\Documents\booksBLSOL\ModDynSyst\Cases_Julia\MoB_Case-3_Gas_Pipeline.ipynb:1](file:///C:/Users/Bernt/OneDrive/Documents/booksBLSOL/ModDynSyst/Cases_Julia/MoB_Case-3_Gas_Pipeline.ipynb:1)

Mar 10 '23 08:03 B-LIE

Ok this is a bug, sorry about this. I will investigate closer soon.

Mar 10 '23 11:03 xtalax

Bugs is a part of life; let me know if you figure it out. In the mean time, I'm racing on with trying to learn the intricacies of ModelingToolkit.

Mar 10 '23 11:03 B-LIE

I am also fighting with a similar bug using MethodOfLines 0.11.1 and 0.11.2. Would it be helpful to contribute another example where it occurs?

ERROR: MethodError: no method matching hasmetadata(::Pair{Num, Float64}, ::Type{Symbolics.VariableDefaultValue})
...
Stacktrace:
 [1] hasdefault(v::Pair{Num, Float64})
   @ ModelingToolkit ~/scratch/twutz/julia_cluster_depots/packages/ModelingToolkit/GJiqn/src/utils.jl:237
 [2] collect_defaults!(defs::Any, vars::Any)
   @ ModelingToolkit ~/scratch/twutz/julia_cluster_depots/packages/ModelingToolkit/GJiqn/src/utils.jl:255 [inlined]
 [3] process_variables!(var_to_name::Dict{Any, Any}, defs::Dict{Any, Any}, vars::Vector{Pair{Num, Float64}})
   @ ModelingToolkit ~/scratch/twutz/julia_cluster_depots/packages/ModelingToolkit/GJiqn/src/utils.jl:248
...

Aug 13 '24 15:08 bgctw

0.11? That must be a typo?

Aug 13 '24 21:08 ChrisRackauckas

Indeed, corrected to MethodOfLines 0.11.1 and 0.11.2. I cannot try 0.11.0 because it does not resolve due other (Symbolics?) constraints.

Aug 14 '24 05:08 bgctw

@bgctw I have this same issue you're describing. I've tried a variety of different package versions (within the compat constraints) and I continue to get the same issue.

I don't have anything helpful to offer besides, "I understand what you're going through."

Aug 20 '24 18:08 bolognam

@bgctw it seems like my issue has been solved by pinning PDEBase for now. Here's what is in my Project.toml compat:

PDEBase = "=0.1.12"

Aug 21 '24 15:08 bolognam

@sathvikbhagavan was there an issue due to your PDEBase release?

Sep 02 '24 03:09 ChrisRackauckas

@bolognam, can I get an MWE for it? I think I know where it might be going wrong. Instead of

@named pdesys = PDESystem(.., [p1 => v1, p2 => v2])

try

@named pdesys = PDESystem(.., [p1, p2]; defaults = Dict(p1 => v1, p2 => v2))

I think MOL also needs a release, try using master.

Sep 03 '24 08:09 sathvikbhagavan

Are the tutorials and such updated to that? That should probably be listed as a breaking change for MethodOfLines.jl, though it's not for PDEBase since this was a non-standard use of the interface in MethodOfLines.jl that is being corrected.

Sep 03 '24 09:09 ChrisRackauckas

No, I will update it.

Sep 03 '24 09:09 sathvikbhagavan