Hello all,
I’ve been playing around with OpenCL and I’ve run into a problem that I can’t pin down. When running this n-body simulation, the amount of memory in use gradually increases with every time step (until I run out completely). I’ve tried using Bytecount[] on every variable that I can think of, and none of them seem to be the culprit.
Anyway, the OpenCL code in question is here:
__kernel void nbody_kern(
float timeStep,
float eps,
__global float4* position1,
__global float4* velocity1,
__global float4* acceleration1,
__global float4* position2,
__global float4* velocity2,
__global float4* acceleration2,
__local float4* localPosition
) {
const float4 dt = (float4)(timeStep,timeStep,timeStep,0.0f);
int idxGlobal = get_global_id(0);
int idxLocal = get_local_id(0);
int globalSize = get_global_size(0);
int localSize = get_local_size(0);
int blocks = globalSize / localSize;
float4 currentPosition = position1[idxGlobal];
float4 currentVelocity = velocity1[idxGlobal];
float4 currentAcceleration = acceleration1[idxGlobal];
float4 newPosition = (float4)(0.0f,0.0f,0.0f,0.0f);
float4 newVelocity = (float4)(0.0f,0.0f,0.0f,0.0f);
float4 newAcceleration = (float4)(0.0f,0.0f,0.0f,0.0f);
for(int currentBlock = 0; currentBlock < blocks; currentBlock++)
{
localPosition[idxLocal] = position1[currentBlock * localSize + idxLocal];
barrier(CLK_LOCAL_MEM_FENCE);
for(int idx = 0; idx < localSize; idx++)
{
float4 dir = localPosition[idx] - currentPosition;
float invRadius = rsqrt(dir.x * dir.x + dir.y * dir.y + dir.z * dir.z + eps);
float magnitude = localPosition[idx].w * invRadius * invRadius * invRadius;
newAcceleration += magnitude * dir;
}
barrier(CLK_LOCAL_MEM_FENCE);
}
// leapfrog integration
newPosition = currentPosition + currentVelocity * dt + 0.5f * currentAcceleration * dt * dt;
newVelocity = currentVelocity + 0.5f * (currentAcceleration + newAcceleration) * dt;
position2[idxGlobal] = newPosition;
velocity2[idxGlobal] = newVelocity;
acceleration2[idxGlobal] = newAcceleration;
}
and the Mathematica code is here:
initialVelocity[pt_] := Module[
{y, v0, \[Theta]},
y = Norm[pt[[1 ;; 3]]]/(2 rd);
v0 = 4*\[Pi]*\[CapitalSigma]0*rd*
y^2 (BesselI[0, y]*BesselK[0, y] - BesselI[1, y]*BesselK[1, y]);
\[Theta] = ArcTan[pt[[1]], pt[[2]]];
{-v0 Sin[\[Theta]], v0 Cos[\[Theta]], 0.0, 0.0} +
RandomVariate[NormalDistribution[0, .001], 4]
]
Needs["OpenCLLink`"];
source = {"D:\\Google Drive\\Mathematica\\nbody.cl"};
stepSimulation = OpenCLFunctionLoad[
source,
"nbody_kern",
{
"Float",
"Float",
{"Float", _, "Input"},
{"Float", _, "Input"},
{"Float", _, "Input"},
{"Float", _, "Output"},
{"Float", _, "Output"},
{"Float", _, "Output"},
{"Local", "Float"}
},
256, "ShellOutputFunction" -> Print
];
dt = .01;
eps = .0001;
rd = 1.0;
\[CapitalSigma]0 = 0.3;
len = 2^14;
particles = (RandomVariate[ExponentialDistribution[rd]]*{Cos[#],
Sin[#], 0.0, 0.0}) & /@ RandomReal[{0, 2 Pi}, len];
particles[[All, 4]] = 1.0/len;
velocity = initialVelocity /@ particles;
acceleration = Array[{0.0, 0.0, 0.0, 0.0} &, len];
position1 = OpenCLMemoryLoad[particles, "Float"];
velocity1 = OpenCLMemoryLoad[velocity, "Float"];
acceleration1 = OpenCLMemoryLoad[acceleration, "Float"];
position2 = OpenCLMemoryAllocate["Float", {len, 4}];
velocity2 = OpenCLMemoryAllocate["Float", {len, 4}];
acceleration2 = OpenCLMemoryAllocate["Float", {len, 4}];
pr = 8;
memList = {MemoryInUse[]};
Graphics3D[{
PointSize[Small],
Point[Dynamic[Refresh[
AppendTo[memList, MemoryInUse[]];
stepSimulation[dt, eps, position1, velocity1, acceleration1,
position2, velocity2, acceleration2, 256*4, len];
stepSimulation[dt, eps, position2, velocity2, acceleration2,
position1, velocity1, acceleration1, 256*4, len];
OpenCLMemoryGet[position1][[All, 1 ;; 3]],
UpdateInterval -> 0]]
]}, PlotRange -> {{-pr, pr}, {-pr, pr}, {-pr/4, pr/4}},
ImageSize -> 1000
]
Dynamic[ListLinePlot[memList]]
As far as I can tell, the code does what it’s supposed to and I’m generally pleased with the results. This is what I get after a few steps:
However, the memory problem limits how long I can run this thing for. This plot shows the amount of memory in use with respect to the number of steps completed in the simulation.
The only way I’ve been able to free up all that memory is by quitting the kernel, which is hardly ideal. I’m guessing that I’m just missing some important aspect of OpenCL programming. If anyone is more familiar with this stuff and has some ideas, I’d be thrilled to hear them.
