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

General Usage
1

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 :confused:.

Does anyone have any suggestions, ideas, tips?

Thank you,
Tiago

2

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

1 Like
3

Look at how DifferentialEquations.jl implements its solution handling interface: Solution Handling · DifferentialEquations.jl

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/solution_interface.jl

2 Likes
4

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

5

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

6

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
7

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