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().