omlins/ParallelStencil.jl

dimensionality restrictions

koehlerson opened this issue · 3 comments

I'm considering to use this package for a PDE formulation where the grid is in a four dimensional space, therefore, I'd like to ask how hard is the dimensionality restriction (at least according to the readme) of 1,2 or 3D present in the code?

Sorry for the late reply. CellArrays can be used for additional dimensions. Currently you need to use the main branch if you want to use this functionality. You can find very basic example here: https://omlins.github.io/CellArrays.jl/dev/examples/memcopyCellArray3D_ParallelStencil/

Thanks for the hint! I am not sure if I understood it correctly. Do I understand the example correctly that it shows a three dimensional array whose values can be represented in a 4x4 way? If so, I was asking more towards if its possible to get another parallel index and another dimension of the underlying array so something like

T  = @zeros(nx, ny, nz, nu);
T2 = @zeros(nx, ny, nz, nu);
Ci = @zeros(nx, ny, nz, nu);

and

@parallel_indices (ix,iy,iz,iu) function copy3D_explicit!(T2::CellArray, T::CellArray, Ci::CellArray)
    T2[ix,iy,iz,iu] = T[ix,iy,iz,iu] + Ci[ix,iy,iz,iu];
    return
end

I'm not sure the example was fully clear to you. In the example it says:

nx, ny, nz = 128, 128, 1024
celldims   = (4, 4)
T  = @zeros(nx, ny, nz, celldims=celldims)

This will create a CellArray, which contains 128x128x1024 cells, where each cell consists of a matrix of 4x4 values. So, in the end it is a five dimensional array.
Operations on these cells follow the standard rules and syntax of array operations in julia.

if nevertheless you really need a 4th parallel index, then you can compute it yourself based on the one-dimensional, two-dimensional or three-dimensional indices. Just like in CUDA, only one-dimensional, two-dimensional or three-dimensional indices can be obtained directly. See for example: https://stackoverflow.com/questions/44368525/using-linear-index-to-map-to-4d-array