Particle systems

The main class. 

Subclass of kinetic particles. 

3D particle with a bodyorientation 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 radiusL[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 cellsi
andj
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 subdictionary
"parameters"
(eventually empty) contains the keyword arguments to be passed to the function found in the dictionarytarget_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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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 tox[i,:]
in the interaction kernel. Parameters
 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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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_ix_j) \left(\frac{(R_i+R_j)^2}{x_ix_j^2}  1\right) (x_ix_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 bylazy_overlapping_kernel()
. Parameters
who ((N,) BoolTensor or None, optional) – The rows of the interaction kernel is
pos[who,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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_ix_j) v_j}{\left\sum_j K(x_ix_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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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 MotschTadmor normalisation is:
\[\Omega_i = \frac{\sum_j K(x_ix_j) v_j}{\sum_j K(x_ix_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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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)
Nonnormalised average velocity.
Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the meanfield average velocity around particle \((x_i,v_i)\) (without normalisation) is:
\[J_i = \frac{1}{N}\sum_j K(x_ix_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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
if with_who is None). Default is None.**kwargs – Arbitrary keyword arguments.
 Returns
The meanfield 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 dimensiond
. Return type
Tensor
 max_kappa(kappa, kappa_max, who=None, with_who=None, **kwargs)
The meanfield target with a norm threshold.
Given a system of particles with positions \((x_i)_i\) and velocities \((v_i)_i\), the meanfield 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_ix_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 dimensiond
. 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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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 Qtensor:
\[Q_i = \sum_j K(x_ix_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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
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,:]
(orvel
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,:]
(orvel
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 bodyorientation in \(SO(3)\).
A bodyorientation 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(ytxz))^T\]\[e_2 = (2(yzxt), x^2  y^2 + z^2  t^2, 2(xy+zt))^T\]\[e_3 = (2(xz+yt), 2(ztxy), 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
Bodyorientation stored as a unit quaternion.
 Type
(N,4) Tensor
 property vel
The velocity is the first column vector of the bodyorientation 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 bodyorientation matrix.
e3 ((N,3) Tensor): Third column vector of the bodyorientation matrix.
 Return type
tuple
 property qtensor
Compute the normalised uniaxial Qtensors of the particles.
If \(q\) is a unit quaternion seen as a colum vector of dimension 4, the normalised uniaxial Qtensor 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
Qtensors 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 bodyorientation matrix.
 normalised(who=None, with_who=None, **kwargs)
Projection of the average bodyorientation on \(SO(3)\).
Given a system of particles with positions \((x_i)_i\) and bodyorientation matrices \((A_i)_i\), the meanfield average bodyorientation 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 Qtensor:
\[\overline{Q}_i = \sum_j K(x_i  x_j) Q_j\]where \(Q_j\) is the Qtensor 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,:]
(orpos
if who is None). Default is None.with_who ((N,) BoolTensor or None, optional) – The columns of the interaction kernel is
pos[with_who,:]
(orpos
if with_who is None). Default is None.**kwargs – Arbitrary keyword arguments.
 Returns
The unit quaternion associated to the projection of the meanfield average bodyorientation 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(N1)} \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 ofvel[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 bodyorientations of the particles which are outside the boundary.