Particle systems

Particles(pos, interaction_radius, box_size)

The main class.

KineticParticles(pos, vel, ...[, ...])

Subclass of kinetic particles.

BOParticles(pos, bo, interaction_radius, ...)

3D particle with a body-orientation in \(SO(3)\).
class sisyphe.particles.Particles(pos, interaction_radius, box_size, vision_angle=6.283185307179586, axis=None, boundary_conditions='periodic', block_sparse_reduction=False, number_of_cells=1024)

The main class.

N

Number of particles.

Type

int

d

Dimension.

Type

int

pos

Positions of the particles.

Type

(N,d) Tensor

R

Radius of the particles.

Type

float or (N,) Tensor

L

Box size.

Type

(d,) Tensor

angle

Vision angle.

Type

float

axis_is_none

True when an axis is specified.

Type

bool

axis

Axis of the particles.

Type

(N,d) Tensor or None

bc

Boundary conditions. Can be one of the following:

  • list of size \(d\) containing 0 (periodic) or 1 (wall with reflecting boundary conditions) for each dimension.

  • "open" : no boundary conditions.

  • "periodic" : periodic boundary conditions.

  • "spherical" : reflecting boundary conditions on the sphere of radius L[0]/2 and center L/2.

Type

list or str

target_method

Dictionary of target methods (see the function compute_target()).

Type

dict

target_option_method

Dictionary of options for the targets (see the function compute_target()).

Type

dict

blocksparse

Use block sparse reduction or not.

Type

bool

iteration

Iteration.

Type

int

centroids

Centroids.

Type

Tensor

eps

Size of the cells.

Type

Tensor

keep

keep[i,j] is True when the cells i and j are contiguous.

Type

BoolTensor

compute_blocksparse_parameters(number_of_cells)

Compute the centroids, the size of the cells and the keep matrix.

Parameters

number_of_cells (int) – The maximal number of cells. It will be rounded to the nearest lower power \(d\).

Returns

(nb_cells, d) Tensor, (nb_cells,) Tensor, (nb_cells, nb_cells) BoolTensor:

The centroids the size of the cells and the keep matrix.

best_blocksparse_parameters(min_nb_cells, max_nb_cells, step=1, nb_calls=100, loop=1, set_best=True)

Returns the optimal block sparse parameters.

This function measures the time to execute nb_calls calls to the __next__() for different values of the block sparse parameters given by the number of cells to the power \(1/d\).

Parameters

  • min_nb_cells (int) – Minimal number of cells to the power \(1/d\).

  • max_nb_cells (int) – Maximal number of cells to the power \(1/d\).

  • step (int, optional) – Step between two numbers of cells. Default is 1.

  • nb_calls (int, optional) – Number of calls to the function __next__() at each test. Default is 100.

  • loop (int, optional) – Number of loops. Default is 1.

  • set_best (bool, optional) – Set the block sparse parameters to the optimal ones. Default is True.

Returns

  • fastest (int): The number of cells which gave the

    fastest result.

  • nb_cells (list): The number of cells tested.

  • average_simu_time (list): The average simulation

    times.

  • simulation_time (numpy array): The total simulation

    times.

Return type

tuple

add_target_method(name, method)

Add a method to the dictionary target_method.

Parameters

  • name (str) – The name of the method.

  • method (func) – A method.

add_target_option_method(name, method)

Add an option to the dictionary target_option_method

Parameters

  • name (str) – The name of the option.

  • method (func) – An option.

compute_target(which_target, which_options, **kwargs)

Compute a target and apply some options.

Parameters

  • which_target (dict) –

    Dictionary of the form:

    {"name"       : "method_name",
 "parameters" : {"param1" : param1,
                 "param2" : param2,
                  ...
                }
}

    The sub-dictionary "parameters" (eventually empty) contains the keyword arguments to be passed to the function found in the dictionary target_method with the keyword “method_name”.

  • which_options (dict) –

    Dictionary of the form:

    {"name_of_option1" : parameters_of_option1,
 "name_of_option2" : parameters_of_option2,
...
}

    The dictionary `parameters_of_option1` contains the keyword arguments to be passed to the method found in the dictionary target_option_method with the keyword “name_of_option1”.

  • **kwargs – Keywords arguments to be passed to the functions which compute the target and the options.

Returns

Compute a target using which_target and then apply the options in which_options successively.

Return type

Tensor

linear_local_average(*to_average, who=None, with_who=None, isaverage=True, kernel=<function lazy_interaction_kernel>)

Fast computation of a linear local average.

A linear local average is a matrix vector product between an interaction kernel (LazyTensor of size (M,N)) and a Tensor of size (N,dimU).

Parameters

  • *to_average (Tensor) – The quantities to average. The first dimension must be equal to N.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • isaverage (bool, optional) – Divide the result by N if isaverage is True. Default is True.

  • kernel (func, optional) – The function which returns the interaction kernel as a LazyTensor. Default is lazy_interaction_kernel().

Returns

The linear local averages.

Return type

tuple of Tensors

nonlinear_local_average(binary_formula, arg1, arg2, who=None, with_who=None, isaverage=True, kernel=<function lazy_interaction_kernel>)

Fast computation of a nonlinear local average.

A nonlinear local average is a reduction of the form

\[\sum_j K_{ij} U_{ij}\]

where \(K_{ij}\) is an interaction kernel ( LazyTensor of size (M,N)) and \(U_{ij}\) is a LazyTensor of size (M,N,dim) computed using a binary formula.

Parameters

  • binary_formula (func) – Takes two arguments arg1 and arg2 and returns a LazyTensor of size (M,N,dim).

  • arg1 ((M,D1) Tensor) – First argument of the binary formula.

  • arg2 ((N,D2) Tensor) – Second argument of the binary formula.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • isaverage (bool, optional) – Divide the result by N if isaverage is True. Default is True.

  • kernel (func, optional) – The function which returns the interaction kernel as a LazyTensor. Default is lazy_interaction_kernel().

Returns

The nonlinear local average.

Return type

(M,dim) Tensor

number_of_neighbours(who=None, with_who=None)

Compute the number of neighbours.

The number of neighbours of x[i,:] is the number of ones in the row corresponding to x[i,:] in the interaction kernel.

Parameters

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

Returns

The number of neighbours.

Return type

(M,) Tensor

morse_target(Ca, la, Cr, lr, p=2, mass=1.0, who=None, with_who=None, local=False, isaverage=True, kernel=<function lazy_interaction_kernel>, **kwargs)

Compute the force exerted on each particle due to the Morse potential.

Parameters

  • Ca (float) – Attraction coefficient.

  • la (float) – Attraction length.

  • Cr (float) – Repulsion coefficient.

  • lr (float) – Repulsion length.

  • p (int, optional) – Exponent. Default is 2.

  • mass ((N,) Tensor or float, optional) – Mass of the particles. Default is 1.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • local (bool, optional) – If local is true then the sum only considers the y particles within the interaction kernel. Default is False.

  • isaverage (bool, optional) – If isaverage is True then divide the total force by N where N is the number of True in :Tensor”with_who. Default is True.

  • kernel (func, optional) – The interaction kernel to use. Default is lazy_interaction_kernel()

  • **kwargs – Arbitrary keyword arguments.

Returns

The force exerted on each particle located at pos[who,:].

Return type

Tensor

quadratic_potential_target(who=None, with_who=None, kernel=<function lazy_interaction_kernel>, isaverage=True, **kwargs)

Compute the force exerted on each particle due to the quadratic potential.

Parameters

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • isaverage (bool, optional) – If isaverage is True then divide the total force by N where N is the number of True in :Tensor”with_who. Default is True.

  • kernel (func, optional) – The interaction kernel to use. Default is lazy_interaction_kernel()

  • **kwargs – Arbitrary keyword arguments.

Returns

The force exerted on each particle located at pos[who,:].

Return type

Tensor

overlapping_repulsion_target(who=None, with_who=None, **kwargs)

Compute the repulsion force exerted on each particle due to the overlapping.

The force exerted on a particle located at \(x_i\) derives from a logarithmic potential and is given by:

\[F = \sum_j K(x_i-x_j) \left(\frac{(R_i+R_j)^2}{|x_i-x_j|^2} - 1\right) (x_i-x_j)\]

where \(x_j\) is the position of the other particles; \(R_i\) and \(R_j\) are the radii of the particles at \(x_i\) and \(x_j\) respectively and \(K\) is the overlapping kernel.

Note

in the present case, it is simpler to write explicitly the reduction formula but the same function could be implemented using a binary formula in nonlinear_local_average() with the kernel given by lazy_overlapping_kernel().

Parameters

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The force exerted on each particle located at pos[who,:].

Return type

Tensor

class sisyphe.particles.KineticParticles(pos, vel, interaction_radius, box_size, vision_angle=6.283185307179586, axis=None, boundary_conditions='periodic', block_sparse_reduction=False, number_of_cells=1024)

Subclass of kinetic particles.

vel

Velocities.

Type

(N,d) Tensor

property axis

If a None axis is given, then the default axis is the normalised velocity.

Returns

(N,d) Tensor

property order_parameter

Compute the norm of the average velocity.

Returns

The order parameter.

Return type

float

check_boundary()

Update the positions and velocities of the particles which are outside the boundary.

normalised(kappa, who=None, with_who=None, **kwargs)

The normalised average velocity.

Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the normalised average velocity around particle \((x_i,v_i)\) is

\[\Omega_i = \frac{\sum_j K(x_i-x_j) v_j}{\left|\sum_j K(x_i-x_j) v_j\right|},\]

where \(K\) is the interaction kernel.

Parameters

  • kappa ((1,) Tensor or (M,) Tensor) – The concentration parameter. If the size of kappa is bigger than 1 then its size must be equal to the number of True values in who.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The normalised averaged velocity multiplied by the concentration parameter.

Return type

Tensor

motsch_tadmor(kappa, who=None, with_who=None, **kwargs)

The average velocity of the neighbours only.

Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the average velocity around particle \((x_i,v_i)\) with the Motsch-Tadmor normalisation is:

\[\Omega_i = \frac{\sum_j K(x_i-x_j) v_j}{\sum_j K(x_i-x_j)},\]

where \(K\) is the interaction kernel.

Parameters

  • kappa ((1,) Tensor or (M,) Tensor) – The concentration parameter. If the size of kappa is bigger than 1 then its size must be equal to the number of True values in who.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The average velocity of the neighbours multiplied by the concentration parameter.

Return type

Tensor

mean_field(kappa, who=None, with_who=None, **kwargs)

Non-normalised average velocity.

Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the mean-field average velocity around particle \((x_i,v_i)\) (without normalisation) is:

\[J_i = \frac{1}{N}\sum_j K(x_i-x_j) v_j,\]

where \(K\) is the interaction kernel.

Parameters

  • kappa ((1,) Tensor or (M,) Tensor) – The concentration parameter. If the size of kappa is bigger than 1 then its size must be equal to the number of True values in who.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The mean-field average velocity multiplied by a scaling parameter and by the concentration parameter. The scaling parameter is equal to the volume of the box divided by the volume of the ball of radius R in dimension d.

Return type

Tensor

max_kappa(kappa, kappa_max, who=None, with_who=None, **kwargs)

The mean-field target with a norm threshold.

Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the mean-field average velocity around particle \((x_i,v_i)\) with a cutoff is:

\[\Omega_i = \frac{1}{\frac{1}{\kappa}+\frac{|J_i|}{\kappa_\mathrm{max}}} J_i\]

where

\[J_i = \frac{1}{N}\sum_j K(x_i-x_j) v_j,\]

and \(K\) is the interaction kernel multiplied by the scaling parameter. The scaling parameter is equal to the volume of the box divided by the volume of the ball of radius R in dimension d. The norm of \(\Omega_i\) is thus bounded by \(\kappa_\mathrm{max}\).

Parameters

  • kappa ((1,) Tensor or (M,) Tensor) – The concentration parameter. If the size of kappa is bigger than 1 then its size must be equal to the number of True values in who.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The cutoff average velocity.

Return type

Tensor

nematic(kappa, who=None, with_who=None, **kwargs)

Normalised average of the nematic velocities.

Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the nematic average velocity around particle \((x_i,v_i)\) is:

\[\Omega_i = (v_i\cdot \overline{\Omega}_i)\overline{\Omega}_i,\]

where \(\overline{\Omega}_i\) is any unit eigenvector associated to the maximal eigenvalue of the average Q-tensor:

\[Q_i = \sum_j K(x_i-x_j) \left(v_j\otimes v_j - \frac{1}{d} I_d\right),\]

where \(K\) is the interaction kernel.

Parameters

  • kappa ((1,) Tensor or (M,) Tensor) – The concentration parameter. If the size of kappa is bigger than 1 then its size must be equal to the number of True values in who.

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The nematic average velocity multiplied by the concentration parameter.

Return type

Tensor

bounded_angular_velocity(target, angvel_max, dt, who=None, **kwargs)

Bounded angle between the velocity and the target.

The maximum angle between the velocity and the new target cannot exceed the value angvel_max * dt. The output tensor new_target is of the form:

new_target = a * vel + b * target

where \(a,b>0\) are such that the cosine between new_target and vel is equal to:

max_cos := max(cos(angvel_max * dt), cos_target)

where cos_target is the cosine between target and vel. The solution is:

b = (1 - max_cos ** 2) / (1 - cos_target ** 2)
a = |target| / |vel| * (max_cos - b * cos_target)
Parameters

  • target ((M,d) Tensor) – Target.

  • angvel_max (float) – Maximum angular velocity.

  • dt (float) – Interval of time.

  • who ((N,) BoolTensor, optional) – The velocities of the particles associated to the target are vel[who,:] (or vel if who is None). Default is None.

  • **kwargs – Arbitrary keywords arguments.

Returns

The new target.

Return type

(M,d) Tensor

pseudo_nematic(target, who=None, **kwargs)

Choose the closest to the velocity between the target and its opposite.

Parameters

  • target ((M,d) Tensor) – Target.

  • who ((N,) BoolTensor, optional) – The velocities of the particles associated to the target are vel[who,:] (or vel if who is None). Default is None.

  • **kwargs – Arbitrary keywords arguments.

Returns

Take the opposite of target if the angle with the velocity is larger than \(\pi/2\).

Return type

(M,d) Tensor

class sisyphe.particles.BOParticles(pos, bo, interaction_radius, box_size, vision_angle=6.283185307179586, axis=None, boundary_conditions='periodic', block_sparse_reduction=False, number_of_cells=1024)

3D particle with a body-orientation in \(SO(3)\).

A body-orientation is a rotation matrix is \(SO(3)\) or equivalently a unit quaternion. If \(q = (x, y, z, t)\) is a unit quaternion, the associated rotation matrix is \(A = [e_1, e_2, e_3]\) where the column vectors are defined by:

\[e_1 = (x^2 + y^2 - z^2 + t^2, 2(xt+yz), 2(yt-xz))^T\]
\[e_2 = (2(yz-xt), x^2 - y^2 + z^2 - t^2, 2(xy+zt))^T\]
\[e_3 = (2(xz+yt), 2(zt-xy), x^2 - y^2 - z^2 + t^2)^T\]

Conversely, if \(A\) is the rotation of angle \(\theta\) around the axis \(u = (u_1, u_2, u_3)^T\) then the associated unit quaternion is:

\[q = \cos(\theta/2) + \sin(\theta/2)(iu_1 + ju_2 + ku_3).\]

or its opposite.

bo

Body-orientation stored as a unit quaternion.

Type

(N,4) Tensor

property vel

The velocity is the first column vector of the body-orientation matrix.

Returns

The velocity.

Return type

(N,3) Tensor

property axis

If a None axis is given, then the default axis is the (normalised) velocity.

Return type

(N,3) Tensor

property orth_basis

An othogonal basis of the orthogonal plane of the velocity.

Returns

  • e2 ((N,3) Tensor): Second column vector of the body-orientation matrix.

  • e3 ((N,3) Tensor): Third column vector of the body-orientation matrix.

Return type

tuple

property qtensor

Compute the normalised uniaxial Q-tensors of the particles.

If \(q\) is a unit quaternion seen as a colum vector of dimension 4, the normalised uniaxial Q-tensor is defined by:

\[Q = q \otimes q - \frac{1}{4}I_4\]

where \(\otimes\) is the tensor product and \(I_4\) the identity matrix of size 4.

Returns

Q-tensors of the N particles.

Return type

(N,4,4) Tensor

property omega_basis

Local frame associated to the direction of motion.

The unit vector \(\Omega\) is the direction of motion (it is equal to vel) and \([\Omega, e_\theta, e_\varphi]\) is the local orthomormal frame associated to the spherical cooordinates.

Returns

  • p ((N,3) Tensor): \(e_\varphi\).

  • q ((N,3) Tensor): \(-e_\theta\).

Return type

tuple

property u_in_omega_basis

Coordinates of \(u\) in the frame \([\Omega, p, q]\) given by omega_basis where \(u\) is the second column vector of the body-orientation matrix.

normalised(who=None, with_who=None, **kwargs)

Projection of the average body-orientation on \(SO(3)\).

Given a system of particles with positions \((x_i)_i\) and body-orientation matrices \((A_i)_i\), the mean-field average body-orientation around the particle \((x_i,A_i)\) is is the arithmetic average:

\[J_i = \sum_j K(x_i - x_j) A_j\]

where \(K\) is the interaction kernel. The projection on \(SO(3)\) is

\[P_{SO(3)}(J_i) = argmax_{A\in SO(3)}( A \cdot J_i )\]

for the dot product

\[A \cdot B = \frac{1}{2} Tr(A^T B).\]

The unit quaternion associated to \(P_{SO(3)}(J_i)\) is any eigenvector associated to the maximal eigenvalue of the average Q-tensor:

\[\overline{Q}_i = \sum_j K(x_i - x_j) Q_j\]

where \(Q_j\) is the Q-tensor of particle \((x_j,A_j)\) given by qtensor.

Parameters

  • who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is pos[who,:] (or pos if who is None). Default is None.

  • with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is pos[with_who,:] (or pos if with_who is None). Default is None.

  • **kwargs – Arbitrary keyword arguments.

Returns

The unit quaternion associated to the projection of the mean-field average body-orientation on \(SO(3)\).

Return type

(M,4) Tensor

property theta

Polar angle of the direction of motion

property phi

Azimuthal angle of the direction of motion

property global_order

Returns

\[\frac{1}{ N(N-1)} \sum_{i,j} (q_i \cdot q_j)^2.\]
build_reflection_quat(who, axis)

Returns the quaternion associated to the rotation of angle \(2\alpha\) and axis \(n\) where \(\alpha\) is the angle between vel[who,:] and the plane orthogonal to axis and \(n\) is the element of this plane obtained by a rotation of \(\pi/2\) around axis of the normalised projection of vel[who,:].

Parameters

  • who ((N,) BoolTensor) –

  • axis (int) – The axis number either 0 (\(x\) axis), 1 (\(y\) axis) or 2 (\(z\) axis).

Returns

Rotation (unit quaternion).

Return type

(M,4) Tensor

check_boundary()

Update the positions and body-orientations of the particles which are outside the boundary.