Kernels

Kernels describe to a looping method what computation should be performed. A ParticleLoop will execute a kernel on each particle in the dictonary of ParticleDat instances that are passed to the constructor. In a ParticleLoop a kernel may only read and write data to the particle currently under consideration or a ScalarArray. A PairLoop kernel follows the same rules as a ParticleLoop kernel except that a PairLoop kernel may also read data from another particle. Kernels must take into account that the order of execution over the set of particles or pairs of particles is not guaranteed.

Constants

We can substitute constant variables for “hardcoded” values by creating a new kernel.Constant. For example to make dht take the value 0.0005 we would write:

dht_c = kernel.Constant('dht', 0.0005)

ParticleLoop and PairLoop constructors expect to be passed a list of constants if such substitutions are required.

ParticleLoop Example

This example is modified from the second step of a Velocity-Verlet integrator. The kernel presented below reads values from ParticleDat data structures labeled F, M and updates a ParticleDat labeled by V. Finally a value is reduced into a ScalarArray labeled K. To access data in a ParticleDat within the kernel we take our chosen name, for example M, and append .i to indicate that this operation concerns the first particle under consideration, to access the second particle under consideration, for example in a PairLoop, we append a .j. In a ParticleLoop we may only consider one particle at a time hence in this example all our data access concerns only one particle by using M.i, V.i, F.i. Components of ParticleDat data structures are indexed through square brackets.

kernel_code = '''
const double M_tmp = 1.0 / M.i[0];

// Update the values in a ParticleDat
V.i[0] += dht * F.i[0] * M_tmp;
V.i[1] += dht * F.i[1] * M_tmp;
V.i[2] += dht * F.i[2] * M_tmp;

To access the ScalarArray data we do not append any .i or .j suffix and instead access the data as if it were a C array.

// Here we reduce into a ScalarArray (kinetic energy example)
K[0] += V.i[0]*V.i[0] + V.i[1]*V.i[1] + V.i[2]*V.i[2]
'''