# Working around current MethodOfLines.jl known limitations

**URL:** <https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346>\
**Category:** Modelling & Simulations\
**Tags:** question, pde, modelingtoolkit\
**Created:** [April 11, 2022, 2:56pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346 "2022-04-11T14:56:29Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jenrei](https://avatars.discourse-cdn.com/v4/letter/j/e8c25b/32.png) [@jenrei](https://discourse.julialang.org/u/jenrei)\
**Post date:** [April 11, 2022, 2:56pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346/1 "2022-04-11T14:56:30Z")

</div>

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 Physics-Informed Learning · MethodOfLines.jl](http://methodoflines.sciml.ai/dev/)  
The original equation would look like:

```julia
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:

```julia
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:

```julia
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:

```julia
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

```julia
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](https://symbolics.juliasymbolics.org/v0.1/tutorials/symbolic_functions/#Registering-Functions-1) But really looking for some help or further docs on this issue. Thanks in advance!

> **heat equation**
>
> ```julia
> 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**
>
> ```julia
> 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**
>
> ```julia
> # 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [April 11, 2022, 3:32pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346/2 "2022-04-11T15:32:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [April 11, 2022, 3:44pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346/3 "2022-04-11T15:44:40Z")

</div>

I would also try the equation in the original form

```julia
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!

---

<div class="post-metadata">

**Author:** ![jenrei](https://avatars.discourse-cdn.com/v4/letter/j/e8c25b/32.png) [@jenrei](https://discourse.julialang.org/u/jenrei)\
**Post date:** [April 11, 2022, 9:09pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346/4 "2022-04-11T21:09:15Z")

</div>

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](https://github.com/SciML/MethodOfLines.jl/issues/78)

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [April 12, 2022, 4:25pm UTC](https://discourse.julialang.org/t/working-around-current-methodoflines-jl-known-limitations/79346/5 "2022-04-12T16:25:58Z")

</div>

Check out the issue for a fix 🙂