MethodOfLines.jl

MethodOfLines.jl copied to clipboard

Hi, I was wondering, is it possible to connect the ODESystem created by this package with other ODESystems, and if so, how?

I'd like to model a wave that on it's left boundary is set to the simple ODE from the MTK tutorial. This works if I add that equation directly to the PDESystem boundary condition equations, but ideally I'd like to specify that boundary later by creating a separate bc_sys::ODESystem and connnecting and composing the two systems.

The example below is basically a modified heat tutorial. It uses the latest release of this package.

After alias elimination both pdesys₊u[1](t) and bc_sys₊u_bc(t) remain in the 11 states, and there are 12 equations, including:

0 ~ pdesys₊u[1](t) - bc_sys₊u_bc(t)
10.0pdesys₊u[2](t) + 2.094(pdesys₊u[1](t)^-0.4)*Differential(t)(pdesys₊u[1](t)) - 10.0pdesys₊u[1](t) ~ 0
10.0pdesys₊u[2](t) + 2.094(pdesys₊u[2](t)^-0.4)*Differential(t)(pdesys₊u[2](t)) - 10.0pdesys₊u[1](t) ~ 0
Differential(t)(bc_sys₊u_bc(t)) ~ 0.3333333333333333 - 0.3333333333333333bc_sys₊u_bc(t)

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables u(..) u_bc(t) = 0.0
Dt = Differential(t)
Dx = Differential(x)

# kinematic wave with hardcoded parameters
eq = Dx(u(t, x)) + 3.49 * 0.6 * u(t, x)^(0.6 - 1) * Dt(u(t, x)) ~ 0

# the left BC is left out here, since we want that to be defined by a separate ODESystem
# which should work similar to adding this equation to `bcs`:
# Dt(u(t, 0)) ~ (1 - u(t, 0)) / 3.0
bcs = [
    u(0, x) ~ 0.5 - abs(x - 0.5)  # triangle, zero at edges
    # no left BC here
    u(t, 1) ~ 0
]

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

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE system into an ODE system
molsys, tspan = symbolic_discretize(pdesys, discretization)

# Create a separate ODESystem to control the boundary
@named bc_sys = ODESystem(Dt(u_bc) ~ (1 - u_bc) / 3.0)

# Connect the two systems and compose
@named _sys = ODESystem([bc_sys.u_bc ~ molsys.u[1]], t, [], [])
@named sys = compose(_sys, molsys, bc_sys)

sim = structural_simplify(sys)
# More equations than variables, here are the potential extra equation(s):
#  0 ~ pdesys₊u[1](t) - bc_sys₊u_bc(t)

Hi, thanks for your interest.

Unfortunately this isn't possible at the moment, and likely won't be for some time until we have better integration with ModelingToolkit, since equations aren't marked by whether they come from a boundary or interior. So unfortunately you will have to create the PDESystem after your ODE BC, and re discretize every time you want to change it.

It may be possible by creating a marker bc equation type that is replaced when compose is called, I will have a think about this.

So unfortunately you will have to create the PDESystem after your ODE BC, and re discretize every time you want to change it.

This sounds like a good solution, but I haven't been able to adapt the example from above to this approach. How would you do it? For this example I can directly add Dt(u(t, 0)) ~ (1 - u(t, 0)) / 3.0 as a BC, but how would I do it if bc_sys is a larger system of which I just want to use one state u_bc as a boundary condition?

You could create a variable u_bc and use it in a bc equation, as long as a boundary i.e. u(t, 0) is referenced in the equation it should work. Try adding a pde equation tying the state to the boundary i.e. u_bc(t) ~ u(t, 0).

Perhaps I don't fully understand, but are you suggesting to add the equation bc_sys.u_bc ~ u(t, 0) both to the PDE equations as well as the boundary conditions? I tried this and got this error at symbolic_discretize:

There must be the same number of equations and unknowns, got 2 equations and 1 unknowns

Adding u_bc to the PDESystem depvars argument [u(t, x), u_bc] doesn't help, since it gets filtered out here; u_bc becomes u(t, 0) in alldepvars, which gets filtered out since it has a number argument, the 0. Here is the full example and error as I tried it:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables u(..) u_bc(t) = 0.0
Dt = Differential(t)
Dx = Differential(x)

# Create a separate ODESystem to control the boundary
@named bc_sys = ODESystem(Dt(u_bc) ~ (1 - u_bc) / 3.0)

# kinematic wave with hardcoded parameters
eqs = [
    Dx(u(t, x)) + 3.49 * 0.6 * u(t, x)^(0.6 - 1) * Dt(u(t, x)) ~ 0
    bc_sys.u_bc ~ u(t, 0)
]

bcs = [
    u(0, x) ~ 0.5 - abs(x - 0.5)  # triangle, zero at edges
    u(t, 0) ~ bc_sys.u_bc
    u(t, 1) ~ 0
]

domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), u_bc])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE system into an ODE system
molsys, tspan = symbolic_discretize(pdesys, discretization)

# Connect the two systems and compose
@named _sys = ODESystem([bc_sys.u_bc ~ molsys.u[1]], t, [], [])
@named sys = compose(_sys, molsys, bc_sys)

sim = structural_simplify(sys)

ERROR: AssertionError: There must be the same number of equations and unknowns, got 2 equations and 1 unknowns
Stacktrace:
 [1] 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, 1, MethodOfLines.CenterAlignedGrid})        
   @ MethodOfLines d:\visser_mn\.julia\dev\MethodOfLines\src\system_parsing\interior_map.jl:18
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines d:\visser_mn\.julia\dev\MethodOfLines\src\MOL_discretization.jl:55

The following is equivalent, you could build up a larger ODE system by adding your ODEs to eqs. Hopefully you will be able to adapt this to your needs

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables u(..) u_bc(..) = 0.0
Dt = Differential(t)
Dx = Differential(x)

# kinematic wave with hardcoded parameters
eqs = [
    Dx(u(t, x)) + 3.49 * 0.6 * u(t, x)^(0.6 - 1) * Dt(u(t, x)) ~ 0,
    Dt(u_bc(t)) ~ (1 - u_bc(t)) / 3.0]

bcs = [
    u(0, x) ~ 0.5 - abs(x - 0.5)  # triangle, zero at edges
    u(t, 0) ~ u_bc(t)
    u(t, 1) ~ 0
]

domains = [t ∈ Interval(0.0, 1.0), x ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem(eqs, bcs, domains, [t, x], [u(t, x), u_bc(t)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t)

# Convert the PDE system into an ODE system
molsys, tspan = symbolic_discretize(pdesys, discretization)

sim = structural_simplify(molsys)

Please let me know if you have further questions.