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

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)
``````

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

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.

``````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