I’m working on a simulation that requires taking FFTs and IFFTs of pretty large systems (~1024x1024x2048 elements). Between transforms I need to advance all of the elements according to a differential equation. Because the arrays are so large I’m using PencilArrays.jl and PencilFFTs.jl to parallelize the transforms using mpi. The parallel/transform part of my code is running really well, but I’m stumbling on the ODE step.

I made a small demo code using normal arrays to test the ODE setup and it works fine.

```
#################### basic ODE test
A = Complex.([0.1,1.0])
dA = Complex.([0.0, 0.0])
p = Complex([1e-1])
function minimal_function(dA, A, p, t)
for i = 1:length(A)
dA[i] = p[1]*A[i]
end
return dA
end
# set up and solve ODE problem
tspan = (0.0,1.0)
f!(du, u, p, t) = minimal_function(du, u, p, t)
prob = ODEProblem{true}(f!, A, tspan, p)
sol = solve(prob)
```

However, when I tried adapting this to work with the pencil arrays I ran into trouble. The code runs just fine and gives a solution object `sol`

but the entries in `sol.u`

are all complex zeros. Any thoughts on what I’m doing wrong?

```
# Attempt using pencil arrays. One mpi node is sufficient to duplicate the problem
using MPI
MPI.Init()
using DifferentialEquations
using PencilArrays
using PencilFFTs
comm = MPI.COMM_WORLD
mpi_size = MPI.Comm_size(comm)
mpi_rank = MPI.Comm_rank(comm)
# set the box dimensions
Nx = 8
Ny = 8
Nz = 16
# divide the space up over the mpi workers
proc_dims = let pdims = zeros(Int, 2)
MPI.Dims_create!(mpi_size, pdims)
pdims[1], pdims[2]
end
# set the transform to in-place, complex-to-complex
transform = Transforms.FFT!()
plan = PencilFFTPlan((Nx,Ny,Nz), transform, proc_dims, comm)
# allocate memory for the pencil array that is compatible with the plan
A = allocate_input(plan)
dA = allocate_input(plan)
# local_R and local_dR are views into the pencil arrays and
# rngR annd rngdR are tuples of vectors of indices used to access the data.
local_R = global_view(first(A))
rngR = axes(local_R)
local_dR = global_view(first(dA))
rngdR = axes(local_dR)
# fill A with some random data
for i in rngR[1], j in rngR[2], k in rngR[3]
local_R[i,j,k] = rand()
end
function minimal_test(local_dS, local_S, p, t)
# apply a non-linearity
for i in rngR[1], j in rngR[2], k in rngR[3]
local_dR[i,j,k] = p[1]*local_R[i,j,k]^2
end
return local_dR
end
# set up and solve ODE problem
tspan = (0.0,1.0)
p = [0.25]
f!(du, u, p, t) = minimal_test(du, u, p, t)
prob = ODEProblem{true}(f!, local_R, tspan, p)
sol = solve(prob)
```