# 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

# ρ = ρ(ξ(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