Coupled PDE Systems

[vpuri3]

Problem

To solve the 2D Burgers equation with OrdinaryDiffEq.jl

https://github.com/vpuri3/PDEInterfaces.jl/blob/d1402040c35b80dc5815e508ea56daf6805b28ba/examples/fourier2d/coupled.jl#L2-L7

I am concatenating vx, vy into a ComponentArray

https://github.com/vpuri3/PDEInterfaces.jl/blob/d1402040c35b80dc5815e508ea56daf6805b28ba/examples/fourier2d/coupled.jl#L60-L66

and encapsulating the right-hand side into separate SciMLOperators for vx, and vy called Dtx, Dty respectively.

https://github.com/vpuri3/PDEInterfaces.jl/blob/d1402040c35b80dc5815e508ea56daf6805b28ba/examples/fourier2d/coupled.jl#L57-L58

My ODE function is now:

https://github.com/vpuri3/PDEInterfaces.jl/blob/d1402040c35b80dc5815e508ea56daf6805b28ba/examples/fourier2d/coupled.jl#L68-L73

The problem in this approach is that the update behaviour of Dtx depends on vy (and Dty depends on vx), but Dtx doesn't ever see vy (and Dty vx) in the ODE function.

Janky solution

The janky solution to the problem would be to modify the ODE function to manually update the SciMLOperators as follows:

https://github.com/vpuri3/PDEInterfaces.jl/blob/7b7d26df4aae491f8b2ea2148985d5f85c1ca7ca/examples/fourier2d/coupled_jank.jl#L63-L75

this is, of-course, not robust as different discretizations and function spaces would place vx, vy differently and this wouldn't extend to 3D easily either.

Quick fix

Even though Dtx, Dty act only on vx, vy respectively, they need access to the entire velocity vector during update_coefficients!. This can be done by supporting keyword arguments in update_coefficients!. The ODEFunction would then look like

https://github.com/vpuri3/PDEInterfaces.jl/blob/347bdcecafc16715c774a079d0a7246897ed1c07/examples/fourier2d/coupled.jl#L63-L71

update_coefficients! can be easily modified to use kwargs. This is both robust as the advection diagonal operator containing velocity can be prescribed advection behaviour as follows:

Vx = DiagonalOperator(zero(x); update_func= (diag,u,p,t;vel=vel) = copy!(diag, vel.vx))

Long term solution

The long term goal should be to just pass a SciMLOperator as your ODE function. So we'd have to write a SciMLOperator that can act on vector of vectors or ComponentArray types.

[vpuri3]

after personal communication with @ChrisRackauckas, best path forward seems to be to dump the entire velocity field in p.

Cx = advectionOp((zero(x), zero(x)), space, discr;
                 vel_update_funcs=(
                                   (v,u,p,t) -> copy!(v, p.vel.vx),
                                   (v,u,p,t) -> copy!(v, p.vel.vy),
                                  )
                )

Cy = advectionOp((zero(x), zero(x)), space, discr;
                 vel_update_funcs=(
                                   (v,u,p,t) -> copy!(v, p.vel.vx),
                                   (v,u,p,t) -> copy!(v, p.vel.vy),
                                  )
                )

function ddt(du, u, p, t)
    ps = ComponentArray(vel=u)
    Dtx(du.vx, u.vx, ps, t)
    Dty(du.vy, u.vy, ps, t)

    du
end