PDE definition cannot handle simple multiplication by a constant?

[SimonEnsemble]

this gives an error, with D₀::Float64.

	diff_eq = ∂t(c(x, t)) ~ D₀ * ∂x(
		 ∂x(c(x, t)) / (1.0 + exp(α * (c(x, t) - χ)))
		)

my hack to get it to work:

	diff_eq = ∂t(c(x, t)) ~ ∂x(
		 ∂x(c(x, t)) / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀)
		)

the error is:

MethodError: no method matching collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Div{Real, SymbolicUtils.Term{Real, Nothing}, SymbolicUtils.Add{Real, Float64, Dict{Any, Number}, Nothing}, Nothing})

Closest candidates are:

collect_ivs_from_nested_differential!(::Any, !Matched::SymbolicUtils.Term)

[ChrisRackauckas]

Care you share a runnable example?

[SimonEnsemble]

sorry @ChrisRackauckas, will do for now on.

here is a minimal example. see the four different WAYs. all are the same mathematically, but two yield the error above.
v4.25.0 in DiffEqOperators.

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, DataFrames, CSV, Statistics, PlutoUI, DifferentialEquations


function toy_pde()
	@parameters x t

	@variables c(..)

	∂t  = Differential(t)
	∂x  = Differential(x)
	∂²x = Differential(x) ^ 2
	
	D₀ = 1.5
	α = 0.15
	χ = 1.2
	R = 0.1
	cₑ = 2.0
	ℓ = 1.0
	Δx = 0.01
	
	# WAY 1: MethodError
	# D = D₀ / (1.0 + exp(α * (c(x, t) - χ)))
	# diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
	
	# WAY 2: no error. it's the D₀ then?
	# D = 1.0 / (1.0 + exp(α * (c(x, t) - χ)))
	# diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))
	
	# WAY 3: no error.
	# diff_eq = ∂t(c(x, t)) ~ ∂x(1.0 / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀) * ∂x(c(x, t)))
	
	# WAY 4: error.
	diff_eq = ∂t(c(x, t)) ~ ∂x(D₀ / (1.0 + exp(α * (c(x, t) - χ))) * ∂x(c(x, t)))

	
	
	bcs = [
		# initial condition
		c(x, 0) ~ 0.0,
		# Robin BC
		∂x(c(0.0, t)) / (1 + exp(α * (c(0.0, t) - χ))) * R * D₀ + cₑ - c(0.0, t) ~ 0.0,		  
		# no flux BC
		∂x(c(ℓ, t)) ~ 0.0]
	

	# define space-time plane
	domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]

	# put it all together into a PDE system
	@named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]);

	# discretize space
	discretization = MOLFiniteDifference([x=>Δx], t)

	# convert the PDE into a system of ODEs via method of lines
	prob = discretize(pdesys, discretization)
	
	# compute number of spatial discretization points, boundaries included.
	Nₓ = length(prob.u0) + 1 # change to +2 for Dirichlet
	@assert (length(0:Δx:ℓ) == Nₓ)

	# solve system of ODEs in time.
	sol = solve(prob, saveat=0.01)
	
	# convert solution to data frame. solution at boundaries must be inferred.
	return DataFrame(sol)
end

toy_pde()

[ChrisRackauckas]

Thanks a ton. That'll be good debugging fodder.

[SimonEnsemble]

@ChrisRackauckas here is another example with a concentration-dependent diffusion coefficient that fails with an error. it works until discretize is called.

 using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, CairoMakie,ColorSchemes, DataFrames, CSV, Statistics, Optim, Polynomials, PlutoUI, DifferentialEquations, SpecialFunctions

function toy_pde2()
	@parameters x t

	@variables c(..)

	∂t  = Differential(t)
	∂x  = Differential(x)
	∂²x = Differential(x) ^ 2
	
	D₀ = 1.2
	c_max = 0.15
	c_e = 0.1
	η = 0.1
	cₑ = 2.0
	ℓ = 1.0
	Δx = 0.01

	# D = D(c)
	θ = c(x, t) / c_max
	β = sqrt(1 - 4 * θ * (1 - θ) * (1 - η))
	ϵ = (2 * θ + β - 1) / (2 * η * (1 - θ))
	D = D₀ * (1 + ϵ) ^ 3 / (1 - θ) / (1 + η * ϵ) ^ 4 * (1 + 2 * (1 - β) / β)

	diff_eq = ∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t)))

	bcs = [
		# initial condition
		c(x, 0) ~ 0.0,
		# Dirichelt BC
		c(0.0, t) ~ c_e,
		# no flux BC
		∂x(c(ℓ, t)) ~ 0.0]
	

	# define space-time plane
	domains = [x ∈ Interval(0.0, ℓ), t ∈ Interval(0.0, 5.0)]

	# put it all together into a PDE system
	@named pdesys = PDESystem(diff_eq, bcs, domains, [x, t], [c(x, t)]);

	# discretize space
	discretization = MOLFiniteDifference([x=>Δx], t)

	# convert the PDE into a system of ODEs via method of lines
	prob = discretize(pdesys, discretization)
end

toy_pde2()

the error is:

MethodError: no method matching collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Div{Real, SymbolicUtils.Mul{Real, Float64, Dict{Any, Number}, Nothing}, SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, Nothing})

Closest candidates are:

collect_ivs_from_nested_differential!(::Any, !Matched::SymbolicUtils.Term) at ~/.julia/packages/ModelingToolkit/r5ZaU/src/utils.jl:160

    collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Term{Real, Nothing})@utils.jl:164
    collect_differentials(::Vector{Symbolics.Equation})@utils.jl:138
    check_equations(::Vector{Symbolics.Equation}, ::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}})@utils.jl:151
    (::ModelingToolkit.var"#ODESystem#92#93")(::Bool, ::Type{ModelingToolkit.ODESystem}, ::Vector{Symbolics.Equation}, ::SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, ::Vector{SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}}}, ::Vector{Any}, ::Dict{Any, Any}, ::Vector{Any}, ::Vector{Symbolics.Equation}, ::Base.RefValue{Vector{Symbolics.Num}}, ::Base.RefValue{Any}, ::Base.RefValue{Any}, ::Base.RefValue{Matrix{Symbolics.Num}}, ::Base.RefValue{Matrix{Symbolics.Num}}, ::Symbol, ::Vector{ModelingToolkit.ODESystem}, ::Dict{Any, Any}, ::Nothing, ::Nothing, ::Nothing, ::Vector{ModelingToolkit.SymbolicContinuousCallback})@odesystem.jl:100
    var"#ODESystem#94"(::Vector{Symbolics.Num}, ::Vector{Symbolics.Equation}, ::Vector{ModelingToolkit.ODESystem}, ::Symbol, ::Dict{Any, Any}, ::Dict{Any, Any}, ::Dict{Symbolics.Num, Float64}, ::Nothing, ::Nothing, ::Nothing, ::Bool, ::Type{ModelingToolkit.ODESystem}, ::Vector{Symbolics.Equation}, ::Symbolics.Num, ::Vector{Symbolics.Num}, ::Vector{Symbolics.Num})@odesystem.jl:151
    symbolic_discretize(::ModelingToolkit.PDESystem, ::DiffEqOperators.MOLFiniteDifference{Vector{Pair{Symbolics.Num, Float64}}, Symbolics.Num})@MOL_discretization.jl:400
    discretize(::ModelingToolkit.PDESystem, ::DiffEqOperators.MOLFiniteDifference{Vector{Pair{Symbolics.Num, Float64}}, Symbolics.Num})@MOL_discretization.jl:406
    toy_pde2()@Other: 45
    top-level scope@Local: 1[inlined]

[SimonEnsemble]

we've isolated the problem as due to division:

	ϵ = (2 * θ + β - 1) # / (2 * η * (1 - θ))
	D = D₀ * (1 + ϵ) ^ 3#  / (1 - θ) / (1 + η * ϵ) ^ 4 * (1 + 2 * (1 - β) / β)

prevents the error.

[SimonEnsemble]

the key was to add expand_derivatives().

diff_eq = expand_derivatives(∂t(c(x, t)) ~ ∂x(D * ∂x(c(x, t))))

makes the error disappear. found this in:
SciML/ModelingToolkit.jl#1197

[xtalax]

This is due to limitations in @rule, in particular the lack of equivalance between *(a^-1, u(t, x), b and /(*(b, u(t, x), a) in the matching. a future PR will fix this, but this gets complicated very quickly as you have to match each possible set of nested functions individually

A workaround would be to expand / in to a flat equivalent * form before matching, but this doesn't avoid the issue of when some u is nested inside a function like exp

[xtalax]

@shashi What should I do about this?