Combining ODE and PDE

[ctessum]

I would like to create an advection-reaction simulation that combines LaGrangian and Eulerian frames of reference, which means combining ODEs and PDEs together in one system:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
using Plots

@parameters x y t
@variables so4(..)
@variables xx(t) yy(t) so2_puff(t)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 7.0

N = 2

dx = (x_max-x_min)/N
dy = (y_max-y_min)/N

# Interaction between Eulerian and Lagrangian frames of reference
puff_conc(x,y,t) = IfElse(x-dx/2 < xx < x+dx/2 && y-dy/2 < yy < y+dy/2, so2_puff(t), 0.0)
@register puff_conc(x, y, t)

# Circular winds.
θ(x,y) = atan(y.-0.5, x.-0.5)
u(x,y) =  -sin(θ(x,y))
v(x,y) = cos(θ(x,y))

k=0.01 # k is reaction rate
eq = [
    # Lagrangian puff model.
    Dt(xx) ~ u(xx,yy),
    Dt(yy) ~ v(xx,yy),
    Dt(so2_puff) ~ -k*so2_puff,

    # Eulerian model.
    Dt(so4(x,y,t)) ~ -u(x,y)*Dx(so4(x,y,t)) - v(x,y)*Dy(so4(x,y,t)) + k*puff_conc(x,y,t),
]

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

# Periodic BCs
bcs = [so4(x,y,t_min) ~ 0.0,
       so4(x_min,y,t) ~ so4(x_max,y,t),
       so4(x,y_min,t) ~ so4(x,y_max,t),
]

@named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[so4(x,y,t), xx, yy, so2_puff])

discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=2, grid_align=center_align)
@time prob = discretize(pdesys,discretization)

I guess there might be a couple of problems with this, but the first one that shows up is this:

The system of equations is:
Equation[Differential(t)(xx(t)) ~ -sin(atan(yy(t) - 0.5, xx(t) - 0.5)), Differential(t)(yy(t)) ~ cos(atan(yy(t) - 0.5, xx(t) - 0.5)), Differential(t)(so2_puff(t)) ~ -0.01so2_puff(t), Differential(t)(so4[2, 2](t)) ~ 0.01puff_conc(0.5, 0.5, t) + 2.0so4[2, 3](t) - 2.0so4[2, 2](t), Differential(t)(so4[3, 2](t)) ~ 0.01puff_conc(1.0, 0.5, t) + 2.0so4[3, 3](t) - 2.0so4[3, 2](t), Differential(t)(so4[2, 3](t)) ~ 0.01puff_conc(0.5, 1.0, t) + 1.2246467991473532e-16so4[2, 2](t) + 2.0so4[3, 3](t) - 2.0so4[2, 3](t), Differential(t)(so4[3, 3](t)) ~ 0.01puff_conc(1.0, 1.0, t) + 1.414213562373095so4[2, 3](t) + 1.4142135623730951so4[3, 2](t) - 2.82842712474619so4[3, 3](t), so4[2, 1](t) ~ so4[2, 3](t), so4[3, 1](t) ~ so4[3, 3](t), so4[1, 2](t) ~ so4[3, 2](t), so4[1, 3](t) ~ so4[3, 3](t), so4[1, 1](t) ~ 0][Base.OneTo(12)]

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

Term{Real, Base.ImmutableDict{DataType, Any}}[so2_puff(t), so4[1, 1](t), so4[2, 1](t), so4[3, 1](t), so4[1, 2](t), so4[2, 2](t), so4[3, 2](t), so4[1, 3](t), so4[2, 3](t), so4[3, 3](t), xx(t), yy(t)][Base.OneTo(12)]
ERROR: ArgumentError: Term{Real, Base.ImmutableDict{DataType, Any}}[xx(t), yy(t), so2_puff(t)] are missing from the variable map.
Stacktrace:
...

Thanks in advance for any information you might have regarding whether this type of simulation should be expected to work, or any suggestions for going about fixing it. Thanks!

[xtalax]

I've not seen anyone get this working, but a few things to try:

1. Don't register puff_conc, should be fine as is.
2. Give xx, yy and so2_puff initial condiitions (This may be the fix).
3. Don't broadcast in θ.

Let me know if that changes anything.

[valentinsulzer]

There is a test that does a simpler version of this:

MethodOfLines.jl/test/pde_systems/MOL_1D_Linear_Diffusion.jl

Lines 742 to 789 in 17f4cbf

	@testset "Test 13: one linear diffusion with mixed BCs, one ODE" begin
	# Method of Manufactured Solutions
	u_exact = (x,t) -> exp.(-t) * sin.(x)
	v_exact = (t) -> exp.(-t)
	@parameters t x
	@variables u(..) v(..)
	Dt = Differential(t)
	Dx = Differential(x)
	Dxx = Dx^2
	# 1D PDE and boundary conditions
	eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
	Dt(v(t)) ~ -v(t)]
	bcs = [u(0,x) ~ sin(x),
	v(0) ~ 1,
	u(t,0) ~ 0,
	Dx(u(t,1)) ~ exp(-t) * cos(1)]
	# Space and time domains
	domains = [t ∈ Interval(0.0,1.0),
	x ∈ Interval(0.0,1.0)]
	# PDE system
	@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t)])
	# Method of lines discretization
	l = 100
	dx = range(0.0,1.0,length=l)
	dx_ = dx[2]-dx[1]
	order = 2
	discretization = MOLFiniteDifference([x=>dx_],t)
	# Convert the PDE problem into an ODE problem
	prob = discretize(pdesys,discretization)
	# Solve ODE problem
	sol = solve(prob,Tsit5(),saveat=0.1)
	x_sol = dx[2:end-1]
	t_sol = sol.t
	# Test against exact solution
	for i in 1:length(sol)
	@test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:length(x_sol)], atol=0.01))
	@test v_exact(t_sol[i]) ≈ sol.u[i][end] atol=0.01
	end
	end

So it is definitely possible

[ctessum]

Thanks! The problem turned out to be that @variables xx(t) yy(t) so2_puff(t) needs to be @variables xx(..) yy(..) so2_puff(..) and

    Dt(xx) ~ u(xx,yy),
    Dt(yy) ~ v(xx,yy),
    Dt(so2_puff) ~ -k*so2_puff,

needs to be

    Dt(xx(t)) ~ u(xx(t),yy(t)),
    Dt(yy(t)) ~ v(xx(t),yy(t)),
    Dt(so2_puff(t)) ~ -k*so2_puff(t),

and then corresponding boundary conditions need to be added.

I'm still having a problem getting the interaction between the two frames of reference to work, but I'll create a separate issue if I'm not able to figure it out. Thanks again!