Working around current MethodOfLines.jl known limitations

Hi, I’m trying to solve the energy/heat equation with temperature dependent material properties (density, enthalpy and thermal conductivity). Therefore I’m running in the issues that are mentioned as (current) known limitations for MethodOfLines.jl (I’m using version 0.2) MethodOfLines.jl: Automated Finite Difference for Phyiscs-Informed Learning · MethodOfLines.jl
The original equation would look like:

eqs = [
    Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x))*Dx(T(t,x)))
]

I’m looking for guidance and best practice how to apply the workarounds mentioned in the docs to my problem.
So far I tired using auxiliary variables:

eqs  = [Dt(e(t,x)) ~ Dx(q(t,x)),
        e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
        q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]

Here I get errors that number of equations and variables are not matching.
The other stategy was applying chain and product rules:

eqs = [ DT(h(T(t,x)))*Dt(T(t,x))*ρ(T(t,x)) + DT(ρ(T(t,x)))*Dt(T(t,x))*h(T(t,x)) ~ λ(T(t,x))*Dxx(T(t,x)) ]

Unfortunately errors still arise:

ERROR: MethodError: no method matching collect_ivs_from_nested_operator!(::Set{Any}, ::SymbolicUtils.Add{Real, Float64, Dict{Any, Number}, Nothing}, ::Type{Differential})        
Closest candidates are:
  collect_ivs_from_nested_operator!(::Any, ::Term, ::Any) at C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:178
Stacktrace:
 [1] collect_ivs_from_nested_operator!(ivs::Set{Any}, x::Term{Real, Nothing}, target_op::Type)
   @ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:182
 [2] collect_ivs(eqs::Vector{Equation}, op::Type)
   @ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:155
 [3] collect_ivs
   @ C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:149 [inlined]  
 [4] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit C:\Users\JentschR\.julia\packages\ModelingToolkit\f218S\src\utils.jl:169
[...]

I assume it’s because

DT(h(T(t,x)))

is still in there. And I have to use @register and predefine derivatives? Like in Symbolic Calculations and Building Fast Parallel Functions · Symbolics.jl But really looking for some help or further docs on this issue. Thanks in advance!

heat equation
using ModelingToolkit, MethodOfLines
using ModelingToolkit: Interval
using Plots

# Variables, parameters, and derivatives
@parameters t x

# Usage of auxilary variables
#@variables T(..) q(..) e(..)
@variables T(..)
Dt = Differential(t)
Dx = Differential(x)
DT = Differential(T)
Dxx = Differential(x)^2

# Loading parametrized material properties and auxilary functions
# ρ = ρ(ξ(T),T)
# h = h(ξ(T),T)
# λ = λ(ξ(T),T)
# ξ(T) = parametrized 7-th order smoothstep function https://en.wikipedia.org/wiki/Smoothstep
# H(x) = Stepfunction
include("matprops.jl");

# Domain edges
t_min, t_max = (0.,200.0)
x_min, x_max = (0.,5.e-3)
# Discretization parameters
dx=0.1e-3;
order=2;

# Original equation
#eqs = [
#    Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x)*Dx(T(t,x))))
#]
# Usage of auxiliary variables
#eqs  = [Dt(e(t,x)) ~ Dx(q(t,x)),
#        e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
#        q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]
#
# Usage of chain and product rules
eqs = [ DT(h(T(t,x)))*Dt(T(t,x))*ρ(T(t,x)) + DT(ρ(T(t,x)))*Dt(T(t,x))*h(T(t,x)) ~ λ(T(t,x))*Dxx(T(t,x)) ]


Tᵢ = 273.15+30.0
Tⱼ = 273.15+50.0
# Initial and boundary conditions
bcs = [T(t_min,x) ~ Tᵢ,
       T(t,x_min) ~ Tᵢ+(Tⱼ-Tᵢ)*H(t-10.0),
       Dx(T(t,x_max)) ~ 0.]

# Space and time domains
domains = [t ∈ Interval(t_min, t_max),
           x ∈ Interval(x_min, x_max)]
# PDE system
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [T(t,x)])
# Method of lines discretization
discretization = MOLFiniteDifference([x=>dx], t; approx_order=order)

prob = discretize(pdesys, discretization)
sol = solve(prob,Tsit5(),abstol=1e-8,saveat=1)
material properties
include("smoothstep.jl")

"""Gumbel-CDF-Funktion"""
function gumbel(T;μ,β)
    return exp(-exp(-(T-μ)/β))
end;

"""Ableitung der Gumbel-Funktion"""
function dgumbel(T;μ=μ₀,β=β₀)
    z = (T-μ)/β
    return (1/β)*exp(-(z+exp(-z)))
end;

"""Heaviside (Step) Funktion"""
function H(t)
       0.5 * (sign(t) + 1)
end

"""Rampen Funktion"""
function R(t;offset=0.0,slope=1.0)
	return slope*t*H(t)+offset
end

# choose between smoothstep or gumbel function
xi = "smooth"

if xi == "gumbel"
	println("Gumbel Funktion")
	μ₀=42.7+273.15
	β₀=6.360467e-1
	ξ(T) = gumbel(T;μ=μ₀,β=β₀)
	dξ(T) = dgumbel(T;μ=μ₀,β=β₀)
else
	println("Smoothstep Funktion")
	T₀ = 40+273.15 
	T₁ = 42+273.5
	ξ(T) = smoothstep(T;edge1=T₀,edge2=T₁)
	dξ(T) = dsmoothstep(T;edge1=T₀,edge2=T₁);
	iξ(T) = ismoothstep(T;edge1=T₀,edge2=T₁);
end

"""Temperaturabhängige Dichte"""
function ρ(T;ξ=ξ)
	ρₛ = 0.86
	ρₗ(T) =  -0.000636906389320033 * T + 0.975879174796585
	return ξ(T) * ρₗ(T) + (1-ξ(T)) * ρₛ
end;

"""Temperaturabhängige spezifische Wärmekapazität"""
function c(T;ξ=ξ,dξ=dξ)
	cₗ = 2.0e3
	cₛ = 2.0e3
	Δhᵥ = 200e3
	return ξ(T) * cₗ + (1-ξ(T)) * cₛ + dξ(T)*Δhᵥ
end

"""Temperaturabhängige Wärmeleitfähigkeit"""
function λ(T;ξ=ξ)
	λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
	λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
	return ξ(T) * λₗ(T) + (1-ξ(T)) * λₛ(T)
end;

"""Ableitung der Temperaturabhängigen Wärmeleitfähigkeit"""
function dλ(T;ξ=ξ)
	λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
	dλₗ(T) = 2*0.00000252159374600214 * T - 0.00178215201558874
	λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
	dλₛ(T) = -2*0.00000469854522539654 * T + 0.00207972452683292
    return dξ(T)*λₗ(T) + dλₗ(T)*ξ(T) - dξ(T)*λₛ(T) + (1-ξ(T))*dλₛ(T)
end;

function β(T;ρ₀)
	ρₛ = 0.86
	ρₗ(T) =  -0.000636906389320033 * T + 0.975879174796585
	dρₗ =  -0.000636906389320033
    return (-1.0/ρ₀)*(dξ(T)*ρₗ(T)+dρₗ*ξ(T)-dξ(T)*ρₛ)
end;

"""Temperaturabhängige spezifische Enthalpie cp = const."""
function h(T;ξ=ξ,href=0.,Tref=35+273.15)
	cₗ = 2.0e3
	cₛ = 2.0e3
	Δh = 200e3
	return href + iξ(T)*cₗ+((T-Tref)-iξ(T))*cₛ+ξ(T)*Δh
end

"""Temperaturleitfähigkeit"""
α(T) = λ(T)/(ρ(T)*c(T));
smoothstep
# Using IfElse package for ModelingToolkit compatability
using IfElse

# Scale, and clamp x to 0..1 range
"""clamp function"""
function clamp(x,edge1,edge2) 
    return (x - edge1) / ( edge2 - edge1)
end

"""derivative of clamp function"""
function dclamp(x,edge1,edge2)
    return 1.0 / ( edge2 - edge1)
end

"""7-th order smoothstep function"""
function smoothstep(x;edge1=0.,edge2=1.)
    xx = clamp(x,edge1,edge2)
    return IfElse.ifelse(xx < 0., 0.,
        IfElse.ifelse(xx > 1., 1.0, -20.0*xx^7.0+70.0*xx^6.0-84.0*xx^5.0+35.0*xx^4.0))
end

"""derivative of 7-th order smoothstep function"""
function dsmoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
dx = dclamp(x,edge1,edge2)
return IfElse.ifelse(xx < 0., 0.,
    IfElse.ifelse(xx > 1., 0., (-140.0*xx^6.0+420.0*xx^5.0-420.0*xx^4.0+140.0*xx^3.0)*dx))
end

"""integral of 7-th order smoothstep function"""
function ismoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
de = edge2 - edge1
iFct(c,j,de) = de*j^(c+1)/(c+1)
return IfElse.ifelse(xx < 0., 0.,
        IfElse.ifelse(xx > 1., (x-edge2)-20.0*de/8.0+70.0*de/7.0-84.0*de/6.0+35.0*de/5.0,
            -20.0*iFct(7.0,xx,de)+70.0*iFct(6.0,xx,de)-84.0*iFct(5.0,xx,de)+35.0*iFct(4.0,xx,de)))
end

The auxiliary variable approach seems the best here - did you also supply initial conditions for the auxillary variables? I would try this, and registering functions first. If that doesn’t work please post an issue with full code and the error as that would be a bug.

I would also try the equation in the original form

eqs = [
    Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x))*Dx(T(t,x)))
]

With registered λ, ρ and h as from inspection I’d expect this to be caught by the nonlinear laplacian handling.

I will add docs for this case as it’s one of the few where we have special handling.

Please let me know how you get on!

Thanks, for the fast reply! I tried with @register, still unsure if doing it right, but with no sucess. I created an issue at github Github Issue

Check out the issue for a fix :slight_smile: