State Containers

The state object couples multiple ParticleDats with a simulation domain that notionally contains the particles. The framework uses a spatial domain decomposition approach for course grain parallelism across MPI ranks. Spatial decomposition requires that when a particle leaves a subdomain owned by a process the data associated with the particle is moved to the new process.

By adding ParticleDat instances to a state object boundary conditions and the movement of data between MPI ranks can be handled automatically by the framework. Instances of state objects follow the convention that npart returns the total number of particles in the system. npart_local returns the number of particles currently stored on the calling MPI rank.

Base Case

The base case state object is designed with NVE ensembles in mind. In the example below we create a state object called A.

from ppmd import *
State = state.State
PBC = domain.BoundaryTypePeriodic

A = State()

Here we add boundary conditions to the state called A.

A.domain = domain.BaseDomainHalo(extent=(10., 10., 10.))
A.domain.boundary_condition = PBC()

The second requirement for a state object is that a particular type of ParticleDat is added called a PositionDat. With a domain and a PositionDat the state object can handle the movement of particle data between processes.

PositionDat = data.PositionDat
A.pos = PositionDat(ncomp=3)

When added ParticleDats to a state class the npart argument should not be passed. As the framework uses spatial domain decomposition the ParticleDat holds valid particle data in the first A.npart_local rows.

print A.pos[:A.npart_local:, ::]

Particle Data Operations

When we add ParticleDats to a state we define which properties particles hold. For example the following code defines properties for positions, velocities, forces and global ids.

A.pos = PositionDat(ncomp=3)
A.vel = ParticleDat(ncomp=3)
A.force = ParticleDat(ncomp=3)
A.gid = ParticleDat(ncomp=1, dtype=ctypes.c_int64)

Particles are added through the State.modify() interface. This interface allows particles to be added and removed. The use of this interface is collective over all MPI ranks. A “window” where particles can be added and removed is created, on all MPI ranks, with:

with A.modify() as m: # this line must be called on all MPI ranks.
    # Add/remove here

Adding Particles

After State.modify() is called particle data is added with a call to add. In the following example with add N=1000 particles with uniform random positions, normal distributed forces and global ids. Particle properties where no values are passed are initialised to zero.

N = 1000
with A.modify() as m:
    if A.domain.comm.rank == 0:
        m.add({
            A.pos: np.random.uniform(-5, 5, size=(N, 3)),
            A.vel: np.random.normal(0, 1, size=(N, 3)),
            A.gid: np.arange(N).reshape((N, 1))
        })

print(A.npart) # should print 1000
print(A.npart_local) # should print the number of particles on this MPI rank

The if A.domain.comm.rank == 0 statement ensures that only rank 0 adds particles. There are no restrictions on which MPI ranks may add particles. Furthermore a MPI rank can add particles anywhere in the simulation domain. Particle data is added on completion of the State.modify() context. Particles are automatically communicated to the MPI rank that owns the subdomain for those particular particles.

Removing Particles

Particles are removed by passing local indices to the remove call of (in this case) m. For example, to remove the first two particles on rank 1:

with A.modify() as m:
    if A.domain.comm.rank == 1:
        m.remove((0, 1))

It is the users responsibility to determine the local index of the particles they wish to remove. Particles are removed on completion of the State.modify() context. Removing particles will force a reordering of particles.

Modifying Particle Data

The data held in a ParticleDat is modified in a similar manner to modifying a State. First a collective call is made to ParticleDat.modify_view() on all MPI ranks. This call returns a NumPy view of size A.npart_local x ParticleDat.ncomp which can be modified safely. For example to draw new random velocities we perform:

with A.vel.modify_view() as m:
    m[:, :] = np.random.normal(0, 1, size=(A.npart_local, 3))

If particle positions are modified then the particles will be communicated between MPI ranks on completion of the ParticleDat.modify_view context. This may well cause particles to change MPI ranks and hence cause a reordering of particle data.

Initialise Data Scatter Example

Scattering data is a legacy method to broadcast particle data across the MPI ranks. Using State.modify() is the replacement.

In this example we will create initial data on rank 0 then scatter that data across available MPI ranks. When scattering data from a rank the total number of particles State.npart should be set on the state object prior to scattering.

import numpy as np
import ctypes
from ppmd import *

# aliases start

State = state.State
ParticleDat = data.ParticleDat
PositionDat = data.PositionDat
PBC = domain.BoundaryTypePeriodic

# aliases end


N = 1000
A = State()

# Total number of particles is set.
A.npart = N

A.domain = domain.BaseDomainHalo(extent=(10., 10., 10.))
A.domain.boundary_condition = PBC()

A.p = PositionDat(ncomp=3)
A.v = ParticleDat(ncomp=3)
A.f = ParticleDat(ncomp=3)
A.gid = ParticleDat(ncomp=1, dtype=ctypes.c_int)


if mpi.MPI_HANDLE.rank == 0:
    A.p[:] = utility.lattice.cubic_lattice((10, 10, 10), (10., 10., 10.))
    A.v[:] = np.random.normal(0.0, 1.0, size=(N, 3))
    A.f[:] = 0.0
    A.gid[:, 0] = np.arange(N)

A.scatter_data_from(0)

The state will use the positions in the PositionDat to filter which data belongs on which subdomain.