Space varying parameter in PDE with ModelingToolkit

General Usage
,
1

Hi all,

I’m trying to solve a simple 1D heat equation. As the future objective is to identify some parameters of this equation, I use ModellingToolkit with MethodOfLines as showcased in the example of the MethodOfLines package.

The only difference is that I would like to solve the temperature profile in two materials, hence some parameters vary in space. I found in the examples a way to make a variable vary in time (Getting Started with ModelingToolkit.jl · ModelingToolkit.jl) but cannot make it work with space.

Here is my current take on the problem (adapted from the example Solving the Heat Equation · MethodOfLines.jl).

using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq


# Parameters, variables, and derivatives
@parameters t x 
@variables α(x) # the space varying parameter
@variables u(..) 
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

α_fun(x) = x >= 80 ? 10 : 20 # its value depends on x
@register_symbolic α_fun(x)

L = 100.
h = 10
Text = 10
tmax = 24*60*60

# 1D PDE and boundary conditions
eq  = [ α ~ α_fun(x),
        Dt(u(t, x)) ~ α*Dxx(u(t, x)) ] # Tried to reproduce the example with forcing function but not sure this is right...
bcs = [u(0,x) ~ 20.,
       Dt(u(t, 0)) ~ 500,
       Dt(u(t, L)) ~ -h*(u(t,L) - Text)] # Heat source on one side, convection on the other

# Space and time domains
domains = [t ∈ Interval(0.0, tmax),
        x ∈ Interval(0., L)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), α])

# Method of lines discretization
# Need a small dx here for accuracy
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization) # This has now been running for 3 hours straight...

I’m quite sure I’m not on the right path but could not find similar examples. Has someone achieved this before?
Thanks!

2

That’s a misnomer that cannot be mathematically true. It’s a shorthand we use in math, but not a truth. Something cannot both be a function and a scalar in the same interpretation. Just do:

α(x) = x >= 80 ? 10 : 20
eq  = [Dt(u(t, x)) ~ α(x)*Dxx(u(t, x)) ]
3

Thanks! I guess I was confused by the time variable forcing function example. This looks indeed more intuitive. Now it seems to go through, but it errors on my convection BC:
ArgumentError: Upper boundary condition Differential(t)(u(t, 100.0)) ~ 100 - 10u(t, 100.0) on time variable is not supported, please use a change of variables t => -τ to make this an initial condition.

I’m not quite sure what is meant by this change of variable, but I’ll investigate further.

4

@xtalax are time-dependent boundary conditions not handled? I don’t recall seeing a test on that, and the error looks like a missing feature.

5

Time dependent BCS do work, but I think this is missed in the tests. I think the condition for that error is erroneously including bcs in space that are on the upper end of the domain when they have time derivs, while it is only supposed to be thrown for final conditions in time. There are a lot of untested features, some refactors and changes added more features than I was intending to so I didn’t get around to testing them, especially with bcs.

Conditions like this one used to work before I changed the way initial conditions were parsed, and added this error message - the underlying discretization hasn’t changed though. This is a 10 line fix, I can get on this once I’m done debugging integrals