Summing the RPY kernel over blobs (class RPYVelocityEvaluator)

class RPYVelocityEvaluator.RPYVelocityEvaluator(a, mu, Npts)

The purpose of this class is to evaluate the RPY velocity on \(N\) blobs in Stokes flow due to \(N\) forces, where the sum is given by the RPY kernel

\[U_i = \sum_{j} M_{RPY}\left(X_i,X_j;\hat{a} \right) F_j.\]
The RPY kernel \(M_{RPY}\) describes the velocity on the surface of one sphere of radius \(\hat{a}\) due to force concentrated on the surface of another sphere of radius \(\hat{a}\); see this paper for details. Here we use the notation \(\hat{a}\) for spheres because we typically use notation \(a\) for the radius of the filaments (although this class does not know about fibers, only blobs).

In this abstract class, the sum is over free space and is computed using quadratic summation. Child classes implement periodic summation (and could be extended to use fast multipole methods in free space)

__init__(a, mu, Npts)

Constructor

Parameters:

  • a (double) – Hydrodynamic radius

  • mu (double) – System viscosity

  • Npts (int) – Number of points

calcBlobTotalVel(ptsxyz, forces, Dom, SpatialData, nThr)

Compute the velocities on the \(N\) points due to the forces there. This class calls a C++ function which does the calculation faster than python (it is still quadratic complexity and slow).

Parameters:

  • ptsxyz (array) – \(N \times 3\) array of point locations

  • forces (array) – \(N \times 3\) array of forces

  • Dom (Domain object) – The domain where the computation is done (not used in this abstract class)

  • SpatialData (SpatialDatabase object) – Used to sort points if we truncate the kernel (not used in this abstract class)

  • nThr (int) – Number of threads

Returns:

\(N \times 3\) array of velocities at the points ptsxyz

Return type:

array

calcMOneHalfW(ptsxyz, Dom, nThr=1)

This method is used to compute \(M[X]^{1/2} W\) when doing Brownian dynamics simulations. In free space (which is this abstract class), this method will actually compute the full dense RPY matrix \(M\), then do \(M^{1/2}\) densely using built-in python code, and multiply this result by a random Gaussian vector.

Parameters:

  • ptsxyz (array) – \(N \times 3\) array of point locations

  • Dom (Domain object) – The domain where the computation is done

  • nThr (integer) – Number of parallel threads

Returns:

\(3N\) array of Brownian velocities \(M[X]^{1/2} W\)

Return type:

array

class RPYVelocityEvaluator.GPUEwaldSplitter(a, mu, xi, PerDom, Npts, xiHalf=None)

This class implements the RPY summation on a triply periodic domain using Ewald splitting. The idea of Ewald splitting or Ewald summation is to split the RPY kernel into a “near field” part, which is short-ranged and can be truncated, and a “far field” part, which is long-ranged, but smooth, and can therefore be done by the NUFFT. In our case, the Ewald splitting has to be “positive,” so that both the near field and far field matrices that result are SPD, and we can compute

\[M^{1/2}W =^d M_{NF}^{1/2}W_1 + M_{FF}^{1/2}W_2\]
where \(W_1\) and \(W_2\) are uncorrelated Gaussian vectors.

More details on the positive split Ewald (PSE) method can be found in here, while details on the implementation (UAMMD by Raul Perez) can be found at the UAMMD_PSE_Python github page. Note that this class requires a GPU to run! While we have a CPU version for deterministic applications, that version is not maintained and does not support Brownian increments.

__init__(a, mu, xi, PerDom, Npts, xiHalf=None)

Constructor

Parameters:

  • a (double) – Hydrodynamic radius

  • mu (double) – System viscosity

  • xi (double) – Ewald parameter, determining how much of the computational cost goes into the near field and far field

  • PerDom (PeriodicDomain object) – The domain on which we do the calculation

  • Npts (int) – Number of points

  • xihalf (double, optional) – Supply if you want to use a separate Ewald parameter when computing \(M^{1/2}\). If not supplied, the same Ewald parameter will be used for \(M^{1/2}\) and \(M\).

calcBlobTotalVel(ptsxyz, forces, Dom, SpatialData, nThr=1)

Compute the velocities on the \(N\) points due to the forces there. This class calls the GPU implementation of PSE. Prior to doing so, it checks that the Ewald truncation distance (in the near field) is smaller than half the simulation box, so that there are not interactions with more than one periodic image (if there are, it increases the Ewald parameter so that the cutoff distance is smaller).

Parameters:

  • ptsxyz (array) – \(N \times 3\) array of point locations

  • forces (array) – \(N \times 3\) array of forces

  • Dom (Domain object) – The domain where the computation is done (to check the Ewald parameter and set the shear strain)

  • SpatialData (SpatialDatabase object) – Not used here, as the GPU does this for us.

  • nThr (int) – Number of threads; not used because the GPU sets this.

Returns:

\(N \times 3\) array of velocities at the points ptsxyz

Return type:

array

calcMOneHalfW(ptsxyz, Dom, nThr=1)

This method is used to call the GPU code to compute \(M[X]^{1/2} W\) when doing Brownian dynamics simulations. The GPU method splits the computation into far field (where the square root can be computed explicitly) and a near field (where the square root can be computed by Lanczos iteration).

Parameters:

  • ptsxyz (array) – \(N \times 3\) array of point locations

  • Dom (Domain object) – The domain where the computation is done (to check the Ewald parameter)

Returns:

\(3N\) array of Brownian velocities \(M[X]^{1/2} W\)

Return type:

array

updateFarFieldArrays()

Updates internal GPU variables if the Ewald parameter has to be enlarged.