Arrays with periodic boundaries

New to Julia
1

Hi! I am a student in the STEM field. I’m not totally new to Julia, but I really have little experience as a programmer.

A problem which I encounter very often is dealing with arrays with periodic/cyclic boundaries (es. working on periodic lattices). I supposed this was not uncommon, so I was surprised to notice that there are no packages dedicated to this purpose (I am aware of FFTViews, but it seemed to me that it was conceived with a particular application in mind).

Another alternative was to implement by myself a module with, say, a PeriodicArray type and overload the relevant functions, but this was really out of my league.

So my question is: what is the proper/julian way to deal with periodic arrays? Is out there any project I am not aware of that can fulfill my needs?

1 Like
2

You might find the extrapolation methods from https://github.com/JuliaMath/Interpolations.jl useful.

3

The proper way is to do a PeriodicArray type like you said. If you can’t do it yourself yet then you just have to hope someone else does it.

In theory it should be quite simple, you just have to overload size, getindex and setindex! (c.f. https://docs.julialang.org/en/release-0.5/manual/interfaces/#indexing), but in practice it can get a bit complicated - like in the FFTViews case - when you do the general case and you want optimal performance.

getindex is the function that get called when you do x[1] (x[1] is the same as getindex(x,1)). setindex is the same thing for x[1] = pi (it calls setindex!(x,pi,1)). So all you need to do is to transform the index with a modulo to make it periodic.

I made a minimal implementation for 1D arrays, it seems to work, but it’s probably not very good:

module tmp
    
    immutable CircularArray{T} <: AbstractArray{T,1}
        data::AbstractArray{T,1}
    end
    
    circindex(i::Int,N::Int) = 1 + mod(i-1,N)
    circindex(I,N::Int) = [circindex(i,N)::Int for i in I]

    Base.size(A::CircularArray) = size(A.data)
    
    Base.getindex(A::CircularArray, i::Int) = getindex(A.data, circindex(i,length(A.data)))
    Base.getindex(A::CircularArray, I) = getindex(A.data, circindex(I,length(A.data)))
    
    Base.setindex!(A::CircularArray, v, i::Int) = setindex!(A.data, v, circindex(i,length(A.data)))
    Base.setindex!(A::CircularArray, v, I) = setindex!(A.data, v, circindex(I,length(A.data)))    
end

c = tmp.CircularArray(rand(10))

@assert c[11] == c[1] 
@assert c[-9:0] == c[1:10]
2 Likes
4

If you’re working with e.g. finite-difference approximations and similar lattice-based computations, often it is more flexible to use e.g. “ghost elements” (also called “ghost points”, “ghost cells”, “overlap regions”, “ghost zones”, …): you have an extra array element at each edge of the array that you update with whatever boundary condition you want as needed (e.g. at each time step). The motivation here is:

  • For typical lattice-based simulations, you only need to access elements adjacent to the boundaries. You don’t want to pay the price of a specialized getindex(a, i) = a.a[i % n] method for every array access just to be able to access one or two extra elements from each end.

  • It is much more flexible. It becomes trivial to implement boundary conditions such as: periodic, Dirichlet, quasi/Bloch periodic (periodic × phase), mirror reflection or other symmetries, domain decomposition (boundary conditions that couple different arrays, e.g. on different processors), and so on.

(See e.g. the section on “overlap regions” in chapter 5 of Principles of Parallel Programming by Lin and Snyder.)

12 Likes
5

@jonathanBieler Indeed I ended up doing something very similar. However I had a problem with broadcasting, for which I don’t know how the implementation works.

@stevengj Thanks for your suggestion. Lately I’ve written a small script for the solution of the traveling salesman problem with simulated annealing, as proposed in Computational Physics by Mark Newman, for fun and practicing with Julia. In this particular case I needed to swap elements of an array with periodic conditions on the boundary, however I think that this can be handled nicely using a swap function and the “ghost elements” technique. This approach works well for Ising simulations too (implementing a flip function that accounts for the boundary conditions).

6

See padarray in JuliaImages