ArgumentError: Differential w.r.t. multiple variables Any[t, ...] are not allowed.

[Centauria]

Hi, I'm trying to solve a system of first order pde's of one variables (time + 1 other variables). I'm using the Method of Lines, but get an error message when trying to discretize the system.

Here is my equation:

eq = [
    Dt(v(x, t)) ~ ((T0 + k * L * (η0 / η(x, t) - 1)) * Dxx(u(x, t)) / (1 + Dx(u(x, t))^2)^2 - c * ρ_air * d * (v(x, t)^2 + w(x, t)^2)) / η0,
    Dt(u(x, t)) ~ v(x, t),
    Dx(Dt(η(x, t))) ~ -w(x, t) * Dx(η(x, t)),
    Dx(u(x, t))^2 ~ (η0 / η(x, t))^2 - 1
]

boundaries:

f(x; b = L/2) = h * min(x/b, (L-x)/(L-b))
bcs = [
    u(0, t) ~ 0.0,
    u(L, t) ~ 0.0,
    u(x, 0) ~ f(x; b=b),
    v(x, 0) ~ 0.0,
    η(x, 0) ~ m / (√(b^2+h^2) + √((L-b)^2+h^2) - L),
    Dt(η(x, 0)) ~ 0.0,
    w(0, t) ~ 0.0,
    w(L, t) ~ 0.0,
    w(x, 0) ~ 0.0,
]

constants:

inch = 2.54e-2
lbs = 0.4536
gravity = 9.80643738977
L = 34inch
T0 = 34lbs * gravity
d = 0.095inch
ρ = 7.9e3 # 7.85-8.05
λ = 0.4
η0 = ρ * π * (d/2)^2 * λ
m = η0 * L
k = 3.5e4
ρ_air = 1.29
c = 0.41
h = 0.01
b = 2L / 3

when discretizing using N = 6 and approx_order=2:

The system of equations is:
Equation[Differential(t)((v(t))[1]) - 69.20036612048055(((67.04165657466264(u(t))[1] - 167.6041414366566(u(t))[2] + 134.08331314932528(u(t))[3] - 33.52082828733132(u(t))[4])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[1]))) / ((1 + (-5.789717461787865(u(t))[1] + 5.789717461787865(u(t))[2])^2)^2) - 0.0012762357((v(t))[1]^2 + (w(t))[1]^2)) ~ 0, Differential(t)((v(t))[2]) - 69.20036612048055(((33.52082828733132(u(t))[1] - 67.04165657466264(u(t))[2] + 33.52082828733132(u(t))[3])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[2]))) / ((1 + (-5.789717461787865(u(t))[1] + 5.789717461787865(u(t))[2])^2)^2) - 0.0012762357((v(t))[2]^2 + (w(t))[2]^2)) ~ 0, Differential(t)((v(t))[3]) - 69.20036612048055(((33.52082828733132(u(t))[2] - 67.04165657466264(u(t))[3] + 33.52082828733132(u(t))[4])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[3]))) / ((1 + (-5.789717461787865(u(t))[2] + 5.789717461787865(u(t))[3])^2)^2) - 0.0012762357((v(t))[3]^2 + (w(t))[3]^2)) ~ 0, Differential(t)((v(t))[4]) - 69.20036612048055(((33.52082828733132(u(t))[3] - 67.04165657466264(u(t))[4] + 33.52082828733132(u(t))[5])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[4]))) / ((1 + (-5.789717461787865(u(t))[3] + 5.789717461787865(u(t))[4])^2)^2) - 0.0012762357((v(t))[4]^2 + (w(t))[4]^2)) ~ 0, Differential(t)((v(t))[5]) - 69.20036612048055(((33.52082828733132(u(t))[4] - 67.04165657466264(u(t))[5] + 33.52082828733132(u(t))[6])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[5]))) / ((1 + (-5.789717461787865(u(t))[4] + 5.789717461787865(u(t))[5])^2)^2) - 0.0012762357((v(t))[5]^2 + (w(t))[5]^2)) ~ 0, Differential(t)((v(t))[6]) - 69.20036612048055(((-33.52082828733132(u(t))[3] + 134.08331314932528(u(t))[4] - 167.6041414366566(u(t))[5] + 67.04165657466264(u(t))[6])*(151.23879999998886 + 30225.999999999996(-1 + 0.01445079059638154 / (η(t))[6]))) / ((1 + (-5.789717461787865(u(t))[5] + 5.789717461787865(u(t))[6])^2)^2) - 0.0012762357((v(t))[6]^2 + (w(t))[6]^2)) ~ 0, -(v(t))[2] + Differential(t)((u(t))[2]) ~ 0, -(v(t))[3] + Differential(t)((u(t))[3]) ~ 0, -(v(t))[4] + Differential(t)((u(t))[4]) ~ 0, -(v(t))[5] + Differential(t)((u(t))[5]) ~ 0, ifelse((w(t))[2] > 0, (w(t))[2]*(-5.789717461787865(η(t))[1] + 5.789717461787865(η(t))[2]), (w(t))[2]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3])) + Differential(0.17271999999999998)(Differential(t)((η(t))[2])) ~ 0, Differential(0.34543999999999997)(Differential(t)((η(t))[3])) + ifelse((w(t))[3] > 0, (w(t))[3]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3]), (w(t))[3]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4])) ~ 0, ifelse((w(t))[4] > 0, (w(t))[4]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4]), (w(t))[4]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5])) + Differential(0.51816)(Differential(t)((η(t))[4])) ~ 0, Differential(0.6908799999999999)(Differential(t)((η(t))[5])) + ifelse((w(t))[5] > 0, (w(t))[5]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5]), (w(t))[5]*(-5.789717461787865(η(t))[5] + 5.789717461787865(η(t))[6])) ~ 0, 1 + (-5.789717461787865(u(t))[1] + 5.789717461787865(u(t))[2])^2 - ((0.01445079059638154 / (η(t))[1])^2) ~ 0, 1 + (-5.789717461787865(u(t))[1] + 5.789717461787865(u(t))[2])^2 - ((0.01445079059638154 / (η(t))[2])^2) ~ 0, 1 + (-5.789717461787865(u(t))[2] + 5.789717461787865(u(t))[3])^2 - ((0.01445079059638154 / (η(t))[3])^2) ~ 0, 1 + (-5.789717461787865(u(t))[3] + 5.789717461787865(u(t))[4])^2 - ((0.01445079059638154 / (η(t))[4])^2) ~ 0, 1 + (-5.789717461787865(u(t))[4] + 5.789717461787865(u(t))[5])^2 - ((0.01445079059638154 / (η(t))[5])^2) ~ 0, 1 + (-5.789717461787865(u(t))[5] + 5.789717461787865(u(t))[6])^2 - ((0.01445079059638154 / (η(t))[6])^2) ~ 0, (u(t))[1] ~ 0.0, (u(t))[6] ~ 0.0, (w(t))[1] ~ 0.0, (w(t))[6] ~ 0.0]

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

in stderr:

ArgumentError: Differential w.r.t. multiple variables Any[t, 0.6908799999999999, 0.17271999999999998, 0.34543999999999997, 0.51816] are not allowed.

Stacktrace:
  [1] check_equations(eqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real})
    @ ModelingToolkit C:\Users\300300\.julia\packages\ModelingToolkit\BsHty\src\utils.jl:178
  [2] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Any}, tspan::Nothing, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 4, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing; checks::Bool)
    @ ModelingToolkit C:\Users\300300\.julia\packages\ModelingToolkit\BsHty\src\systems\diffeqs\odesystem.jl:154
  [3] ODESystem
    @ C:\Users\300300\.julia\packages\ModelingToolkit\BsHty\src\systems\diffeqs\odesystem.jl:143 [inlined]
  [4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, 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{SymbolicUtils.BasicSymbolic{Real}, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 4, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing)
    @ ModelingToolkit C:\Users\300300\.julia\packages\ModelingToolkit\BsHty\src\systems\diffeqs\odesystem.jl:218
  [5] generate_system(disc_state::PDEBase.EquationState, s::MethodOfLines.DiscreteSpace{1, 4, MethodOfLines.CenterAlignedGrid}, u0::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}, tspan::Tuple{Float64, Float64}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 4, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, disc::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase C:\Users\300300\.julia\packages\PDEBase\Aqj4G\src\discretization_state.jl:41
  [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase C:\Users\300300\.julia\packages\PDEBase\Aqj4G\src\symbolic_discretize.jl:89
  [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PDEBase C:\Users\300300\.julia\packages\PDEBase\Aqj4G\src\discretization_state.jl:58
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase C:\Users\300300\.julia\packages\PDEBase\Aqj4G\src\discretization_state.jl:55
  [9] top-level scope
    @ .\timing.jl:273 [inlined]
 [10] top-level scope
    @ .\In[85]:0

Thank you for your time and help in troubleshooting this model.

[cmhyett]

Hey @Centauria, I wrote your code into a minimal working example, and was able to reproduce your error. Please let me know if you see inconsistencies in this formulation.

using ModelingToolkit, MethodOfLines, DomainSets;

inch = 2.54e-2
lbs = 0.4536
gravity = 9.80643738977
L = 34inch
T0 = 34lbs * gravity
d = 0.095inch
ρ = 7.9e3 # 7.85-8.05
λ = 0.4
η0 = ρ * π * (d/2)^2 * λ
m = η0 * L
k = 3.5e4
ρ_air = 1.29
c = 0.41
h = 0.01
b = 2L / 3

@parameters t x
@variables u(..) v(..) w(..) η(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs = [ Dt(v(x, t)) ~ ((T0 + k * L * (η0 / η(x, t) - 1)) * Dxx(u(x, t)) / (1 + Dx(u(x, t))^2)^2 - c * ρ_air * d * (v(x, t)^2 + w(x, t)^2)) / η0,
        Dt(u(x, t)) ~ v(x, t),
        Dx(Dt(η(x, t))) ~ -w(x, t) * Dx(η(x, t)),
        Dx(u(x, t))^2 ~ (η0 / η(x, t))^2 - 1];

f(x; b = L/2) = h * min(x/b, (L-x)/(L-b))
bcs = [ u(0, t) ~ 0.0,
        u(L, t) ~ 0.0,
        u(x, 0) ~ f(x; b=b),
        v(x, 0) ~ 0.0,
        η(x, 0) ~ m / (√(b^2+h^2) + √((L-b)^2+h^2) - L),
        Dt(η(x, 0)) ~ 0.0,
        w(0, t) ~ 0.0,
        w(L, t) ~ 0.0,
        w(x, 0) ~ 0.0]

domains = [t in Interval(0.0, 1.0),
           x in Interval(0.0, L)]

@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [u(t,x), v(t,x), w(t,x), η(t,x)])

dx = 0.1;
order = 2;
discretization = MOLFiniteDifference([x=>6], t);
prob = discretize(pdesys, discretization)

[cmhyett]

As far as I can tell, this error is arising when attempting to discretize the mixed term Dx(Dt(η(x, t))). I show the output of a debug run of prob = ... at PDEBase/src/symbolic_discretize.jl:80

1|debug> 
In symbolic_discretize(pdesys, discretization) at /home/cmhyett/.julia/packages/PDEBase/Aqj4G/src/symbolic_discretize.jl:7
 63  
 64      ####
 65      # Loop over equations, Discretizing them and their dependent variables' boundary conditions
 66      ####
 67      for pde in pdeeqs
 68          # Read the dependent variables on both sides of the equation
 69          depvars_lhs = get_depvars(pde.lhs, v.depvar_ops)
 70          depvars_rhs = get_depvars(pde.rhs, v.depvar_ops)
 71          depvars = collect(depvars_lhs ∪ depvars_rhs)
 72          depvars = filter(u -> !any(map(x -> x isa Number, arguments(u))), depvars)
 73  
 74          eqvar = get_eqvar(vareqmap, pde)
 75  
 76          # * Assumes that all variables in the equation have same dimensionality except edgevals
 77          args = ivs(eqvar, v)
 78          indexmap = Dict([args[i] => i for i in 1:length(args)])
 79              # Generate the equations for the interior points
 80          discretize_equation!(disc_state, pde, vareqmap, eqvar, bcmap,
 81                               depvars, s, derivweights, indexmap, discretization)
>82      end

1|julia> disc_state.eqs[11:end]
4-element Vector{Equation}:
 ifelse((w(t))[2] > 0, (w(t))[2]*(-5.789717461787865(η(t))[1] + 5.789717461787865(η(t))[2]), (w(t))[2]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3])) + Differential(0.17271999999999998)(Differential(t)((η(t))[2])) ~ 0
 Differential(0.34543999999999997)(Differential(t)((η(t))[3])) + ifelse((w(t))[3] > 0, (w(t))[3]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3]), (w(t))[3]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4])) ~ 0
 ifelse((w(t))[4] > 0, (w(t))[4]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4]), (w(t))[4]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5])) + Differential(0.51816)(Differential(t)((η(t))[4])) ~ 0
 Differential(0.6908799999999999)(Differential(t)((η(t))[5])) + ifelse((w(t))[5] > 0, (w(t))[5]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5]), (w(t))[5]*(-5.789717461787865(η(t))[5] + 5.789717461787865(η(t))[6])) ~ 0

I'm unsure if MoL has support for these sort of terms... @xtalax have you tried this?