How to solve a combination of ODEs and PDEs and organize the solution

General Usage
,
1

Is it possible to use a struct to hold the dependent variables when using the solvers in DifferentialEquations.jl, or is it needed to create a single vector u of all dependent variables?

For example, I have a system with 3 ODEs

  • dXt/dt = …
  • dSt/dt = …
  • dLf/dt = …

and 2 PDES that depend on one spatial dimension

  • dPb(z)/dt = …
  • dSb(z)/dt = …

I’m solving these by discretizing the spatial derivatives in the PDEs by creating a grid with N grid cells and using finite differences to create N ODEs representing each PDE. This gives a total of 3 + 2*N ODEs that are organized into a single vector solved with DifferentialEquations.jl. Everything works, but the solution is difficult to work with since it requires splitting the variables into Xt, St, Lf, Pb, and Sb. I should note that it’s actually a little more complicated, but this is the general idea.

It would be much easier to have a struct that organizes the solution, e.g.,

struct sol 
    Xt :: Float
    St :: Float
    Lf :: Float
    Pb :: Vector{Float}
    Sb :: Vector{Float}
end
2

I have used ComponentArrays to do similar things, I only used an ODE solver but I suspect it works with the PDE solvers as well.

1 Like
3

I usually just use a bigger vector, like is shown in the showcase example:

https://docs.sciml.ai/Overview/stable/showcase/gpu_spde/

However, a great library for writing this kind of code is ComponentArrays.jl which then allows the pieces to be named.

Writing your PDE+ODE with component arrays is pretty natural and gives an efficient form.

1 Like
4

Thanks @contradict and @ChrisRackauckas for the suggestion to use ComponentArrays. I have implemented it in my code and it works well.

I did find it non-trival to access the parts of the DiffEq solution using the component labels. Ultimately, I found the label2index() function which is part of ComponentArrays can be used to find the appropriate indices which can be used with the provided solution analysis tools or plot recipes in DifferentialEquations.jl.

If there is a way to access the components without using label2index(), please share!

Below is a MWE that highlights how this can work.


using ComponentArrays
using OrdinaryDiffEq
using Parameters: @unpack
using Plots

# Define function to get index of 
# each labeled component in solution
idx(sol,var) = label2index(sol.u[1],var)

function odetest()
    # Define and solve 
    # ODE dXo/dt = Xo  and 
    # PDE d(Xp=A,B,C)/dt = Xp
    # ---------------------

    # Define right-hand-side of ODEs
    function dsol_dt!(dsol, sol, p, t)
        @unpack Xo,Xp = sol
        dsol.Xo  = Xo
        dsol.Xp  = Xp
        return nothing
    end

    # Initial conditions
    Xo₀ = 1.0  # 
    Xp₀ =  ComponentArray(A = 1.0, B = 2.0, C = 3.0)
    sol₀ = ComponentArray(Xo=Xo₀,Xp=Xp₀)

    # Setup problem and solve
    tspan = (0.,1.0)
    prob = ODEProblem(dsol_dt!,sol₀,tspan)
    sol = solve(prob,Tsit5(),)


    # Access final solution 
    # --------------------- 
    
    # Access components
    Xo = sol[idx(sol,"Xo"),:]
    Xp = sol[idx(sol,"Xp"),:]

    # Access components with interpolation
    t = 0.0:0.1:1.0
    Xo_dense = sol(t,idxs = idx(sol,"Xo"))
    Xp_dense = sol(t,idxs = idx(sol,"Xp"))

    # Plot solution and components
    p1 = plot(sol                 ,title="Full solution")
    p2 = plot(sol,idxs = idx(sol,"Xo"),title="Only Xo")
    p3 = plot(sol,idxs = idx(sol,"Xp"),title="Only Xp")
    p  = plot(p1,p2,p3,layout=(1,3))
    display(p)
end
odetest()

odetest

1 Like
5

sol[end].Xo?

6

That does work to access the solution at the final time, but access more than one time, e.g. at all the times with Xp = sol[1:end].Xp, does not work. This gives the error ERROR: type ODESolution has no field Xp

7

It works on each i, so just use a compeehension

1 Like