Issue MethodOfLines.jl Burger equation with periodic BCs

Specific Domains Modelling & Simulations
, , ,
1

Hello,

I would like to solve the 1D inviscid Burgers equation \frac{\partial u(x,t)}{\partial t} + \frac{1}{2} \frac{\partial u(x,t)^2}{\partial x} = 0 with periodic BCs on [0,1], i.e. u(0, t) = u(1,t), using MethodOfLines.jl.

I have started from the implementation of the inviscid Burgers equation on the interval [0,1] with Dirichlet BCs in the test folder of MethodOfLines.

However, I am getting an error message regarding the interface boundaries in the periodic case. Here is the code with the error message:

using DomainSets
using MethodOfLines
using ModelingToolkit

@parameters x t
@variables u(..)
Dx = Differential(x)
Dt = Differential(t)
x_min = 0.0
x_max = 1.0
t_min = 0.0
t_max = 6.0

eq = Dt(u(t, x)) ~ -0.5 * Dx(u(t, x)^2)

bcs = [u(0, x) ~ (1 + sin(2*π*x)),
          u(t, x_min) ~ u(t, x_max)]

domains = [t ∈ Interval(t_min, t_max),
    x ∈ Interval(x_min, x_max)]

dx = 0.05

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

disc = MOLFiniteDifference([x => dx], t, advection_scheme=WENOScheme())

prob = discretize(pdesys, disc)

sol = solve(prob, Tsit5())

x_disc = sol[x]
solu = sol[u(t, x)]

I get the following error message:

Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.
└ @ MethodOfLines ~/.julia/packages/MethodOfLines/JKuVA/src/system_parsing/pde_system_transformation.jl:39

The system of equations is:
Equation[0.5Differential(0.05)((u(t))[2]^2) + Differential(t)((u(t))[2]) ~ 0, 0.5Differential(0.1)((u(t))[3]^2) + Differential(t)((u(t))[3]) ~ 0, 0.5Differential(0.15)((u(t))[4]^2) + Differential(t)((u(t))[4]) ~ 0, Differential(t)((u(t))[5]) + 0.5Differential(0.2)((u(t))[5]^2) ~ 0, 0.5Differential(0.25)((u(t))[6]^2) + Differential(t)((u(t))[6]) ~ 0, Differential(t)((u(t))[7]) + 0.5Differential(0.3)((u(t))[7]^2) ~ 0, 0.5Differential(0.35)((u(t))[8]^2) + Differential(t)((u(t))[8]) ~ 0, 0.5Differential(0.4)((u(t))[9]^2) + Differential(t)((u(t))[9]) ~ 0, 0.5Differential(0.45)((u(t))[10]^2) + Differential(t)((u(t))[10]) ~ 0, Differential(t)((u(t))[11]) + 0.5Differential(0.5)((u(t))[11]^2) ~ 0, 0.5Differential(0.55)((u(t))[12]^2) + Differential(t)((u(t))[12]) ~ 0, 0.5Differential(0.6)((u(t))[13]^2) + Differential(t)((u(t))[13]) ~ 0, 0.5Differential(0.65)((u(t))[14]^2) + Differential(t)((u(t))[14]) ~ 0, 0.5Differential(0.7)((u(t))[15]^2) + Differential(t)((u(t))[15]) ~ 0, 0.5Differential(0.75)((u(t))[16]^2) + Differential(t)((u(t))[16]) ~ 0, Differential(t)((u(t))[17]) + 0.5Differential(0.8)((u(t))[17]^2) ~ 0, 0.5Differential(0.85)((u(t))[18]^2) + Differential(t)((u(t))[18]) ~ 0, Differential(t)((u(t))[19]) + 0.5Differential(0.9)((u(t))[19]^2) ~ 0, Differential(t)((u(t))[20]) + 0.5Differential(0.95)((u(t))[20]^2) ~ 0, Differential(t)((u(t))[21]) + 0.5Differential(1.0)((u(t))[21]^2) ~ 0, (u(t))[1] ~ (u(t))[21]]

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

ArgumentError: Differential w.r.t. multiple variables Any[0.95, 0.3, 0.45, 0.25, 0.35, 1.0, 0.7, 0.85, 0.15, 0.2, t, 0.9, 0.65, 0.8, 0.05, 0.5, 0.55, 0.1, 0.4, 0.75, 0.6] are not allowed.

Stacktrace:
 [1] check_equations(eqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/BsHty/src/utils.jl:178
 [2] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Any}, tspan::Nothing, 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, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/BsHty/src/systems/diffeqs/odesystem.jl:154
 [3] ODESystem
   @ ~/.julia/packages/ModelingToolkit/BsHty/src/systems/diffeqs/odesystem.jl:143 [inlined]
 [4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{SymbolicUtils.BasicSymbolic{Real}}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, gui_metadata::Nothing)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/BsHty/src/systems/diffeqs/odesystem.jl:218
 [5] generate_system(disc_state::PDEBase.EquationState, s::MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, u0::Vector{Pair{SymbolicUtils.BasicSymbolic{Real}, Float64}}, tspan::Tuple{Float64, Float64}, metadata::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, disc::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ PDEBase ~/.julia/packages/PDEBase/2OwkC/src/discretization_state.jl:41
 [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ PDEBase ~/.julia/packages/PDEBase/2OwkC/src/symbolic_discretize.jl:88
 [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ PDEBase ~/.julia/packages/PDEBase/2OwkC/src/discretization_state.jl:58
 [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ PDEBase ~/.julia/packages/PDEBase/2OwkC/src/discretization_state.jl:55
 [9] top-level scope
   @ In[13]:28
2

You need to define an auxiliary variable u_sq(t, x) ~ u(t,x)^2 and use that in your equation instead.

1 Like
3

Thank you for your help. I have tried to add a new variable u_sq and add the algebraic equation u_sq(t,x) ~ u(t,x)^2.
Can you clarify how to modify the code?

using DomainSets
using MethodOfLines
using ModelingToolkit

@parameters x t
@variables u(..), u_sq(..)
Dx = Differential(x)
Dt = Differential(t)
x_min = 0.0
x_max = 1.0
t_min = 0.0
t_max = 6.0

eq =  [u_sq(t,x) ~ u(t,x)^2, 
       Dt(u(t, x)) ~ -0.5 * Dx(u_sq(t,x))]

bcs = [u(0, x) ~ (1 + sin(2*π*x)),
          u(t, x_min) ~ u(t, x_max)]

domains = [t ∈ Interval(t_min, t_max),
    x ∈ Interval(x_min, x_max)]

dx = 0.05

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

disc = MOLFiniteDifference([x => dx], t, advection_scheme=WENOScheme())

prob = discretize(pdesys, disc)

sol = solve(prob, Tsit5())

x_disc = sol[x]
solu = sol[u(t, x)]
4

You need to add, u_sq(t, x) to the PDESystem depvars, then this will work. You may also need to give it a periodic condition but I think it should be fine.

5

Hello @xtalax,

I have added these modifications but I still get the same error message (see code below). Any chance that you could share with me a working code for this setting?

using DomainSets
using MethodOfLines
using ModelingToolkit
using OrdinaryDiffEq

@parameters x t
@variables u(..), u_sq(..)
Dx = Differential(x)
Dt = Differential(t)

x_min = 0.0
x_max = 1.0
t_min = 0.0
t_max = 1.0

eq =  [u_sq(t,x) ~ u(t,x)^2, 
       Dt(u(t, x)) ~ -0.5 * Dx(u_sq(t,x))]

bcs = [u(0, x) ~ (1 + sin(2*π*x)),
       u(t, x_min) ~ u(t, x_max),
       u_sq(t, x_min) ~ u_sq(t, x_max)]

domains = [t ∈ Interval(t_min, t_max),
    x ∈ Interval(x_min, x_max)]

dx = 0.001

@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x), u_sq(t, x)])

disc = MOLFiniteDifference([x => dx], t, advection_scheme=WENOScheme())

prob = discretize(pdesys, disc)

 Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.
└ @ MethodOfLines ~/.julia/packages/MethodOfLines/JKuVA/src/system_parsing/pde_system_transformation.jl:39

Code to solve after discretization of the PDE:

sol = solve(prob, SSPRK104(), dt = 1e-2)
plot(sol[u(t,x)][11,:])
6

Your code works for me, can you show me the output of ] st? Try updating packages