How do I construct the initial condition in FEM for ODEProblem once, not twice?

General Usage
, , ,
1

I am looking to obtain a solution approximation u = u(x, t) for a partial differential equation via FEM. I am using DifferentialEquations with the prospect utilizing its parameter estimation functionality.

The parameters determine the spatial discretization, which in turn determines the mass and stiffness matrices involved. As a MWE,

function pde!(du, u, p, t)
    # Spatial Grid
    x = spatial_grid(p)

    # FEM Matrices
    M, S = fem_matrices(x)

    # Solve
    du[:] = M\(S*u)
end

for Julia functions spatial_grid and fem_matrices self-explanatory.

In order to define the ODEProblem I need an initial condition u0, which depends on the grid defined inside the pde! function.

So my question is: Is there a more elegant way to obtain the spatial grid x than having to run x = spatial_grid(p) twice?

It would be nice to somehow fetch the grid or define the initial conditions inside pde!, since the method for producing the grid inside pde! is much more complicated than is shown in the MWE above. Another option is declaring a global which I would love to avoid.

The current implementation idea would be as follows.

function spatial_grid(p::AbstractVector{T}) where T <: Real
    ...
end
function fem_matrices(x::AbstractVector{T}) where T <: Real
    ...
end
function pde!(du, u, p, t)
    x = spatial_grid(p)
    M, S = fem_matrices(x)
    du[:] = M\(S*u)
end
function f0(x::Real)
    ...
end
p = [...]
tspan = (0.0, 30.0)
x = spatial_grid(p)
u0 = f0.(x)
prob = ODEProblem(pde!, u0, tspan, p)
sol = solve(prob)
2

Use a closure to enclose it into the function

pde!(du,u,p,t) = pde!(du,u,p,t,x)

and define x outside.

2 Likes
3

I’m confused. Which variable(s) are you talking about becoming closed? I couldn’t figure out what you meant.

Like this?

function solve_pde(p)
    function pde!(du, u, p, t, x)
        M, S = fem_matrices(x)
        du[:] = M\(S*u)
    end
    x = spatial_grid(p)
    pde!(du, u, p, t) = pde!(du, u, p, t, x)
    u0 = f0.(x)
    tspan = (0.0, 30.0)
    prob = ODEProblem(pde!, u0, tspan, p)
    sol = solve(prob)
end

where x is defined outside? Though this means that when solve(...) passes p to pde!(du, u, p, t), it doesn’t use p at all inside it, which doesn’t make the function ready for parameter estimation.

Or like this?

function pde!(du, u, p, t, x)
    M, S = fem_matrices(x)
    du[:] = M\(S*u)
end

function pde!(du, u, p, t)
    x = spatial_grid(p)
    pde!(du, u, p, t, x)
end

which I think isn’t an implementation of the concept of closure, and I don’t see how this allows me to define initial conditions.

4 
x = spatial_grid(p)
function pde!(du, u, p, t)
    M, S = fem_matrices(x)
    du[:] = M\(S*u)
end

and such. You can then use a callable type to enclose it in a type-stable way.

1 Like