# Suggestions to handle, smartly, a solution of a stochastic PDE

**URL:** <https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [November 5, 2020, 11:49pm UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650 "2020-11-05T23:49:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 5, 2020, 11:49pm UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/1 "2020-11-05T23:49:16Z")

</div>

Hello,

I’m currently building a stochastic PDE solver based on a specific algorithm, under the theme of my master’s thesis.  
The main idea behind my code is simple, given some inputs (time span, physical domain, number of simulations, etc) provided by a user, the code must return a numerical solution computed in the time instants previously specified. Since I’m working with 2D solutions, my output will be one matrix for each instant for each simulation.

Thus, I want my output to be in such way that the time instants are linked with the respective matrix.

```julia
sol = solveSPDE(inputs)
sol[t1] = V1 # V1 is the solution correspondent to t1

```

I don’t know how to do this in a clever and optimized way.  
First, I thought that `sol` could be an object of some structure that I would define, but in that case I would lose that “friendly” syntax, mentioned above, and use the `dot` to refer to the fields of the object,

```julia
sol.t[1] = t1
sol.V[1] = V1

```

Then I came up with the idea of using a dictionary to store my solution and time instants with the respective link between them.  
But as I said, I’m dealing with matrices (approx 500x500), stored in 4/5 time steps for each simulation. This is a lot of data and probably using dictionaries in this situation would be a mess in terms of performance.

Probably, I could try to define a `struct` that somehow have similar properties to dictionaries. But I’m not getting there 😕.

Does anyone have any suggestions, ideas, tips?

Thank you,  
Tiago

---

<div class="post-metadata">

**Author:** ![martenlienen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martenlienen/32/18572_2.png) [@martenlienen](https://discourse.julialang.org/u/martenlienen)\
**Post date:** [November 6, 2020, 12:02am UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/2 "2020-11-06T00:02:04Z")

</div>

You can implement `getindex` for your struct to get the `sol[t1] # => v1` syntax for accessing the solutions.

---

<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:** [November 6, 2020, 1:02am UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/3 "2020-11-06T01:02:23Z")

</div>

Look at how DifferentialEquations.jl implements its solution handling interface: [Solution Handling · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/basics/solution/)

You can use RecursiveArrayTools.jl’s DiffEqArray by just matching the interface (defining the portions of your struct as `.u` and `.t`) and then all of the indexing exists, add a call function so it’s a continuous interpolation (linear for SPDEs), add retcodes, automated plot recipes, etc. This might help you get started:

[https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/ode\_solutions.jl#L1-L19](https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/ode_solutions.jl#L1-L19)  
[https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/solution\_interface.jl](https://github.com/SciML/DiffEqBase.jl/blob/master/src/solutions/solution_interface.jl)

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 6, 2020, 11:20am UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/4 "2020-11-06T11:20:35Z")

</div>

Hi Chris!

Thanks for your response. Actually, after having seen one of your videos about the `DifferentialEquations.jl` package, inspired me to create my own “user friendly” solver for the Stochastic Neural-Field Delayed Equations that I’m studying.

I’ve to admit, before I posted this topic here, I had already seen the code that you’ve posted here, and didn’t understand as well as I wished.  
For example, i don’t even understand these lines:

```julia
(sol::ODESolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::ODESolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

```

Probably, you’re treating the object `sol` as a function too in order to accept some arguments to give an output. I searched in the documentation to find some information about it and found the [Case study: Rational](https://docs.julialang.org/en/v1/manual/constructors/index.html#Case-Study:-Rational) in the official documentation:

```julia
julia> struct OurRational{T<:Integer} <: Real
           num::T
           den::T
           function OurRational{T}(num::T, den::T) where T<:Integer
               if num == 0 && den == 0
                    error("invalid rational: 0//0")
               end
               g = gcd(den, num)
               num = div(num, g)
               den = div(den, g)
               new(num, den)
           end
       end

OurRational(n::T, d::T) where {T<:Integer} = OurRational{T}(n,d)

OurRational(n::Integer, d::Integer) = OurRational(promote(n,d)...)
...

```

I’m not sure if is this concept that you used after defining the structure `ODESolution` or not.

I’ll look to `RecursiveArrayTools.jl` as you said, looks really promising after reading about it.

Thank you so much for your time.  
Cheers,  
Tiago

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [November 6, 2020, 11:31am UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/5 "2020-11-06T11:31:15Z")

</div>

Hello,  
Can you be a little more specific? Before I wrote this post, I had already seen the function `getindex` but, unfortunately I could not figured out how to use it in a structure to do what I’m trying to achieve.

It seems that returns only the indexes of the specified elements, as can be seen here

```julia
  julia> A = [1 2; 3 4]
  2×2 Array{Int64,2}:
   1 2
   3 4
  
  julia> getindex(A, 1)
  1
  
  julia> getindex(A, [2, 1])
  2-element Array{Int64,1}:
   3
   1
  
  julia> getindex(A, 2:4)
  3-element Array{Int64,1}:
   3
   2
   4

```

Using this function I don’t how to assign the time instant to the respective indexes of the matrix.

Thank you,  
Tiago

---

<div class="post-metadata">

**Author:** ![martenlienen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martenlienen/32/18572_2.png) [@martenlienen](https://discourse.julialang.org/u/martenlienen)\
**Post date:** [November 6, 2020, 4:27pm UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/6 "2020-11-06T16:27:26Z")

</div>

You can extend the function `Base.getindex` with a method for your own struct to achieve the interface you want.

```julia
julia> struct A
           a::Vector{Int}
       end

julia> a = A([5, 10, 3, 1])
A([5, 10, 3, 1])

julia> Base.getindex(a::A, i) = a.a[i]

julia> a[4]
1

```

Also that syntax in your other reply that you don’t understand defines the function call syntax for objects.

```julia
julia> (b::A)() = sum(b.a)

julia> a()
19

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [November 6, 2020, 4:43pm UTC](https://discourse.julialang.org/t/suggestions-to-handle-smartly-a-solution-of-a-stochastic-pde/49650/7 "2020-11-06T16:43:59Z")

</div>

My choice for an interface for trajectories `b` would be to allow call overload for interpolation/evaluation at time `b(t)` and then providing `getindex` not directly on trajectories `b` but on `values(b)` or `pairs(b)`, e.g. `values(b)[end]`. Like in the lightweight [https://github.com/mschauer/Trajectories.jl](https://github.com/mschauer/Trajectories.jl), see the remark about iteration in the readme (which also applies to getindex)