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

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.

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

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

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

1 Like

Look at how DifferentialEquations.jl implements its solution handling interface: 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:

2 Likes

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:

``````(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 in the official documentation:

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

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

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

``````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> (b::A)() = sum(b.a)

julia> a()
19
``````
2 Likes

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, see the remark about iteration in the readme (which also applies to getindex)

1 Like