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

**URL:** <https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688>\
**Category:** General Usage\
**Tags:** question, package, diffeq, fem\
**Created:** [October 20, 2020, 2:16pm UTC](https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688 "2020-10-20T14:16:32Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 20, 2020, 2:16pm UTC](https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688/1 "2020-10-20T14:16:32Z")

</div>

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,

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 20, 2020, 3:18pm UTC](https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688/2 "2020-10-20T15:18:58Z")

</div>

Use a closure to enclose it into the function

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

and define `x` outside.

---

<div class="post-metadata">

**Author:** ![kapple](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kapple/32/218915_2.png) [@kapple](https://discourse.julialang.org/u/kapple)\
**Post date:** [October 21, 2020, 9:36am UTC](https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688/3 "2020-10-21T09:36:50Z")

</div>

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

Like this?

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

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 21, 2020, 4:20pm UTC](https://discourse.julialang.org/t/how-do-i-construct-the-initial-condition-in-fem-for-odeproblem-once-not-twice/48688/4 "2020-10-21T16:20:16Z")

</div>

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