heat equation with variable material properties

[jenrei]

In accordance with this forum entry about solving the heat equation with variable material the problem as bug entry with full code. Two approaches are included

when trying the original equation

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

I get this:

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

ERROR: MethodError: no method matching collect_ivs_from_nested_operator!(::Set{Any}, ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, ::Type{Differential})
Closest candidates are:
  collect_ivs_from_nested_operator!(::Any, ::Term, ::Any) at ~/.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 ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:182
 [2] collect_ivs(eqs::Vector{Equation}, op::Type)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:155
 [3] collect_ivs
   @ ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:149 [inlined]
 [4] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:169
 [5] ODESystem(deqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}}, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Any}, 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, connections::Nothing, preface::Nothing, events::Vector{ModelingToolkit.SymbolicContinuousCallback}, tearing_state::Nothing, substitutions::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:120
 [6] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:174
 [7] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [9] top-level scope

With auxilary 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))]

the error is:

ERROR: AssertionError: Variable q(t, x) does not appear in any equation, therefore cannot be solved for
Stacktrace:
 [1] build_variable_mapping(m::Matrix{Int64}, vars::Vector{Any}, pdes::Vector{Equation})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/interiormap.jl:69
 [2] MethodOfLines.InteriorMap(pdes::Vector{Equation}, boundarymap::Dict{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}}}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/interiormap.jl:20
 [3] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:52
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [5] top-level scope
   @ ~/Julia/PCMfluiddynamics/heateq.jl:165

Full code:

using IfElse
using ModelingToolkit, MethodOfLines
using ModelingToolkit: Interval
using Symbolics

"""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

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

T₀ = 40+273.15
T₁ = 42+273.15
ξ(T) = smoothstep(T;edge1=T₀,edge2=T₁)
dξ(T) = dsmoothstep(T;edge1=T₀,edge2=T₁)
iξ(T) = ismoothstep(T;edge1=T₀,edge2=T₁)

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

"""derivitative of density"""
function dρ(T;ξ=ξ)
	ρₛ = 0.86
	ρₗ(T) =  -0.000636906389320033 * T + 0.975879174796585
	dρₗ(T) =  -0.000636906389320033
    return dξ(T) * ρₗ(T) + dρₗ(T) * ξ(T) - dξ(T) * ρₛ
end

"""thermal conductivity"""
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

"""derivative of thermal conductivity"""
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

"""specific 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

"""derivative specific enthalpie cp = const."""
function dh(T;ξ=ξ)
	cₗ = 2.0e3
	cₛ = 2.0e3
	Δh = 200e3
    return ξ(T)*cₗ + cₛ - ξ(T)*cₛ + dξ(T)*Δh
end

# variables, parameters, and derivatives
@parameters t x

# 
@variables T(..)
# when using of auxilary variables
@variables T(..) q(..) e(..)

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

@register ρ(T)
@register λ(T)
@register h(T)
@register dρ(T)
@register dλ(T)
@register dh(T)

Symbolics.derivative(::typeof(ρ), args::NTuple{1,Any}, ::Val{1}) = dρ(args[1])
Symbolics.derivative(::typeof(λ), args::NTuple{1,Any}, ::Val{1}) = dλ(args[1])
Symbolics.derivative(::typeof(h), args::NTuple{1,Any}, ::Val{1}) = dh(args[1])

# 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))]

Tᵢ = 273.15+30.0
Tⱼ = 273.15+50.0
# when using auxilary variables
eᵢ = ρ(Tᵢ)*h(Tᵢ)
qᵢ = 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.]
# when using auxilary variables
bcs = [T(t_min,x) ~ Tᵢ,
       T(t,x_min) ~ Tᵢ+(Tⱼ-Tᵢ)*H(t-10.0),
       Dx(T(t,x_max)) ~ 0.,
       e(t_min,x) ~ eᵢ,
       q(t_min,x) ~ 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)])
# when using auxilary variables
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [T(t,x),e(t,x),q(t,x)])

# method of lines discretization
discretization = MOLFiniteDifference([x=>dx], t; approx_order=order)
prob = ModelingToolkit.discretize(pdesys, discretization)
sol = solve(prob,Tsit5(),abstol=1e-8,saveat=1)

[xtalax]

Thank you, this has revealed 2 unknown issues, first with the nonlinear handling, second with the heuristic used to match equations to variables to be solved for - this will be important debugging fodder. I will try to get a fix in as soon as possible.

[xtalax]

@jenri I have a fix for this, check out the branch in the above pull request to see. Tests are failing as they need updating, so should already work. Hopefully this will be merged soon. This is for the auxiliary variable variant, as I realised that we have no special scheme for cases where a registered function wraps the inner derivative.

There are a few adjustments needed to your code to get this to work.
Firstly, please supply boundary conditions at either end of the domain for q.
I noticed that Dx(q(t, x_max)) ~ 0, so you just need to find one for the lower end.
Secondly, please use a solver that can use mass matricies such as Rosenbrock23().