Differential Calculation Functions

This is, essentially, the core module of FINDS. This is where the calculations that govern the mechanics of the particle simulation occur, and this is written in accordance with the prior work of Mohamed Niged Mabrouk and Dr. Daniel Floryan’s past work with the far-field model of swimmers [3, 4].

The core goal of this module is the computation of the time-derivative of the system matrix \(\mathbf{X}\). Using this calculation, numerical integration schemes are implemented under the ‘simulation’ module to allow for time-stepping.

The following is a simple definition of the far-field model as stated in these papers alongside descriptions of how this model was implemented in the program [1].

[1]

Please note that no in-depth justification for this model will be given in this documentation! A simple mathematical justification is that the leading-order term of the multipole expansion for a body of constant volume in fluid is a dipole, but many more justifying reasons are stated in the papers of Mabrouk and Floryan [3, 4].

The Far-Field Model of Fish

Each fish-particle consists of two fundamental variables, a center-of-mass position \(\mathbf{x}_c\) and an orientation unit vector \(\hat{n}\), and the fish can be represented in full as an array or tuple that stores both of these variables, \(\mathbf{s} = (\mathbf{x}_c, \hat{n})\), and the system of swimmers \(\mathbf{X}\) can be defined as a matrix consisting of each of these fish vectors as rows. Thus, for a system of size \(N\), the system matrix is \(N \times 6\) or of shape \((N,6)\).

The fundamental approximation of the far-field model is treating fish as self-propelled dipoles moving in free space, with the head of the fish being the source and the tail of the fish being the sink. The head and the tail for each fish are separated \(\ell\) apart, or the length of the fish.

Mathematically, these two ends can differ only up to a sign (so the head and tail also form positive and negative ends respectively) due to the conservation of mass of the fluid (the two ends must cancel each other out). Due to this duality, the sources and sinks can be considered more generally as “features” of a fish.

We store all of the fish in the system in a matrix \(\mathbf{X}\) where each fish vector comprises a row, and this is called the system (or state) matrix.

Computing the Derivative of the System Matrix

For a given fish \(\mathbf{X}_i\), the features (front/source and back/sink) can be computed as

(1)\[\begin{split}\mathbf{x}_{f,i} &= \mathbf{x}_{c,i} + \vec{\delta}_i \\ \mathbf{x}_b &= \mathbf{x}_{c,i} - \vec{\delta}_i\end{split}\]

where

(2)\[\vec{\delta}_i = \frac{1}{2} \ell \hat{n}_i\]

and this is returned as the following matrix

(3)\[\begin{split}\mathbf{F} = \begin{bmatrix} x_{f1} & y_{f1} & z_{f1} & x_{b1} & y_{b1} & z_{b1} \\ x_{f2} & y_{f2} & z_{f2} & x_{b2} & y_{b2} & z_{b2} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ x_{fN} & y_{fN} & z_{fN} & x_{bN} & y_{bN} & z_{bN} \end{bmatrix}\end{split}\]

representing the feature positions.

The core equation governing the motion of bodies within the simulation is the induced velocity for a particle a displacement \(\mathbf{r}\) away from a feature (some source or sink of a fish),

(4)\[\mathbf{u} = \pm \frac{\sigma}{4\pi} \frac{\mathbf{r}}{|\mathbf{r}|^3}\]

where \(\sigma\) is the volumetric flow rate of the source or sink. Sources produce outward velocities as they repel bodies outward and sinks produce inward velocities as they attract bodies toward them.

In essence, this equation is used to develop the differential time-derivative of the full state, as the velocity for a single feature is defined as the total velocity contribution from all other features in the system, and the feature velocity and orientation of a fish define the velocity of its orientation.

This idea leads to the feature velocity equation

(5)\[\mathbf{v}_{\alpha,i} = U \hat{n}_i + \frac{\sigma}{4\pi} \mathbf{h}_{\alpha,i}\]

(see finds.calculations.calculate_feature_velocities()) that computes the velocity of some feature \(\alpha\) (which would be either the front or back) for a fish. This is comprised of two terms, an internal contribution to the velocity resulting from the interaction between this feature and its polar opposite on the same fish (also known as the self-propelled velocity) defined as

(6)\[U \hat{n}_i = \frac{\sigma}{4 \pi \ell^2} \hat{n}_i,\]

and the external contribution term that encodes the total velocity contribution from all features (sources and sinks) external to this fish, where the vectors \(\mathbf{h}\) and \(\mathbf{c}\) are defined to represent the value of the original term \(\mathbf{r}/r^3\), which I call the “interaction” between features.

In this context, \(\mathbf{h}_{\alpha,i}\) refers to the total interaction between feature \(\alpha\) of fish \(i\) with all of the other features (barring its polar opposite on fish \(i\)) within the greater system \(\mathbf{X}\),

(7)\[\frac{\sigma}{4\pi} \mathbf{h}_{\alpha,i} = \sum_{j \neq i}^N \frac{\sigma}{4\pi} \mathbf{h}_{\alpha,ij}.\]

The sum \(\sum_{j \neq i}^N \dots\) is the core bottleneck of this program, as this computation for each fish produces an \(\mathcal{O}(N^2)\) runtime when fully pairwise (see finds.calculations.compute_interaction_pairwise()). Thus, this is where the system diverges into using a tree-based computation method, the Barnes-Hut algorithm (see finds.calculations.compute_interaction_barnes_hut() and Barnes-Hut Clustering Approximation).

Furthermore, \(\mathbf{h}_{\alpha,ij}\) refers to the total interaction between feature \(\alpha\) of fish \(i\) and both features of fish \(j\) (both the front and the back), and this is always defined as the individual interaction to the front of fish \(j\) minus the individual interaction to the back of fish \(j\) due to the aforementioned sign change between sources and sinks,

(8)\[\mathbf{h}_{\alpha,ij} = \mathbf{c}_{\alpha f} - \mathbf{c}_{\alpha b}.\]

(see finds.calculations.calculate_fish_interaction()) Lastly, \(\mathbf{c}_{\alpha \beta}\) is defined as the individual interaction between feature \(\alpha\) of fish \(i\) and feature \(\beta\) of fish \(j\) within \(\mathbf{X}\), and this is defined simply as

(9)\[\mathbf{c}_{\alpha\beta} = \frac{\mathbf{r}_{\alpha\beta}}{|\mathbf{r}_{\alpha\beta}|^3}\]

(see finds.calculations.calculate_feature_interaction()) where \(\mathbf{r}\) is defined as the displacement between the two features

(10)\[\mathbf{r}_{\alpha\beta} = \mathbf{x}_{\alpha,i} - \mathbf{x}_{\beta,j}.\]

With this, the derivative for one fish can be defined as the derivative of the center-of-mass position and the orientation

(11)\[\dot{\mathbf{X}}_i = \left(\dot{\mathbf{x}}_{ci}, \frac{d \hat{n_i}}{dt}\right),\]

(see finds.calculations.calculate_system_derivative()) and equivalently, \(\dot{\mathbf{x}}_c\) is the translational velocity and \(\frac{d \hat{n}}{dt}\) is the rotational velocity of the fish. The translational velocity for a fish is just the average of the front and back velocities,

(12)\[\dot{\mathbf{x}}_c = \frac{\mathbf{v}_f + \mathbf{v}_b}{2}\]

and the rotational velocity is defined as

(13)\[\frac{d \hat{n}}{dt} = \frac{\Delta \mathbf{v} + 2 \lambda \hat{n}}{\ell}\]

where \(\Delta \mathbf{v} = \mathbf{v}_f - \mathbf{v}_b\) is the difference in velocity between the front and back of the fish and

(14)\[\lambda = \frac{-\Delta \mathbf{v} \cdot \hat{n}}{2}\]

is a Lagrange multiplier used to ensure that the length of the fish \(\ell\) is kept constant.

Barnes-Hut Clustering Approximation

The Barnes-Hut algorithm simplifies calculating the interaction vectors by clustering fish that are sufficiently “far away” from a given fish (as determined by the Barnes-Hut ratio \(\theta\)). A complete, interactive description of the Barnes-Hut algorithm can be found online [2].

This begins by constructing an Octree that partitions the three-dimensional space around the origin. The fish-particles are inserted in list order, and for each additional point, the Octree expands by further subdividing the three- dimensional space. Then, once the Octree is fully built, each node of the tree (essentially representing every possible division of the three-dimensional space) is “clustered,” meaning that several of the fish-particles are computed into a fish-particle representing the average of all of them.

For example, if we have fish \(\mathbf{X}_i\) and fish \(\mathbf{X}_j\) in \(\mathbf{X}\), their clustered form would simply be the average of the two, \((\mathbf{X}_i + \mathbf{X}_j)/2\). Then, we traverse this tree for each fish. If the node is a leaf, then we automatically compute the interaction. For each branch node in the octree, we then compute a size-to- distance ratio \(\phi\). This is computed as the fraction of the “size” or the side length of the cube comprising the subdivision over the distance from the center-of-mass of the subdivision (or the position of the cluster).

For example, if we have a subdivision storing fish 1 through \(k\) clustered into a cluster-fish with position \(\mathbf{x}_{c}\) with a side length of size \(l\), then the Barnes-Hut ratio with respect to a fish at position \(\mathbf{x}_{c0}\) would be

(15)\[\phi = \frac{l}{||\mathbf{x}_c - \mathbf{x}_{c0}||}.\]

Once we’ve calculated \(\phi\) for the given node, if \(\phi \ge \theta\), then we continue travering the tree to the children of the node. If \(\phi < \theta\), then we compute the interaction.

Functions

finds.calculations.calculate_feature_positions(system)

Computes the front and back positions for each fish in the system per equation (1) in format (3).

Parameters:

system (NDArray) – The system.

Returns:

The matrix storing the head and tail positions for each fish.

Return type:

NDArray

class finds.calculations.OctreeNode(center, side_length)

Bases: object

Represents the node of an Octree, representing the octant partition and the clustered data for this node.

Variables:

• center (NDArray) – The center position of this octant.

• side_length (float) – The side length of this octant.

• children (list[OctreeNode]) – The list of child octants. The size should always be kept constant at 8.

• cluster (NDArray) – The sum of all of the fish stored at or underneath this node.

• cluster_size (int) – The number of fish stored at this node.

• average (NDArray) – The average of all of the fish clustered at this node, defined as the cluster sum divided by the cluster size.

• front_pos (NDArray) – The position of the front of the clustered fish for this octant.

• back_pos (NDArray) – The position of the back of the clustered fish for this octant.

• is_leaf (bool) – Whether or not this node has children.

children

center

side_length

is_leaf

cluster_size

cluster

data

calculate_octant_index(position)

Computes the index of the corresponding octant that a particle located at position would be found in relative to the center of this octant.

Parameters:

position (NDArray) – The position of the object we’re comparing against the center of this octant.

Returns:

The integer index corresponding to the octant the position belongs in, between 0 to 7.

Return type:

int

The logic behind this is simple, and follows binary. Given the center position \(\mathbf{r} = \langle r_x, r_y, r_z \rangle\) of this octant and the comparison point \(\mathbf{p} = \langle p_x, p_y p_z \rangle\), we loop through each of the dimensions \(i\). For each dimension, we evaluate whether or not \(p_i > r_i\), and if so, we set the binary digit at that index to 1. Otherwise, that digit is set to 0. For brevity, rather than computing directly in binary, this same effect can be produced by adding \(2^i\) if \(p_i > r_i\) is true.

Each of the octants map in the following manner:

Octant	Binary	Index
-x, -y, -z	000	0
+x, -y, -z	001	1
-x, +y, -z	010	2
+x, +y, -z	011	3
-x, -y, +z	100	4
+x, -y, +z	101	5
-x, +y, +z	110	6
+x, +y, +z	111	7

For example, given a center point \(\mathbf{r} = \langle 0, 0, 0 \rangle\) and a test point \(\mathbf{p} = \langle 0, -1, 1 \rangle\), the predicate evaluated on each index produces a vector \(\langle \text{False}, \text{False}, \text{True} \rangle\), which corresponds to the binary string 100 (as binary is built in reverse) where 1 is True and 0 is False. From there, the base 10 representation of 100 is 4, so therefore, \(\mathbf{p}\) would be said to lie in octant 4, which corresponds to the –+ octant.

calculate_child_center(child_octant_index)

Computes the center position of a corresponding child octant relative to the center of this parent octant.

Parameters:

child_octant_index (int) – The octant index of the child octant.

Returns:

The center position of the child octant.

Return type:

NDArray

First, the offset from the parent center position is calculated as the side length for this octant divided by four [2]. Then, the child’s octant index can be reinterpreted as binary and looped through for each dimension. For each dimension, if the predicate bit corresponding to that dimension is true [3], then it sets the offset vector for that dimension to have the positive offset. Otherwise, it sets the offset vector at that dimension to have the negative offset.

For example, if the side length of this octant is \(\ell\), and we want to compute the center of an octant with index

[2]

This octant is a cube and the child octant is an octant cube contained within the parent, thus all child octant center positions will be offset from the parent center position by exactly one-fourth of the side length on each dimension.

[3]

This can be directly calculated with base-10 using the bitwise AND operator. Rather than doing these calculations directly in binary in this code, if we have an octant index \(o\) on dimension \(i\), the bit being True or False can be checked by performing \(o\) & \(2^i\).

insert_into_children(fish)

Inserts a new fish into the child nodes of this octant.

Parameters:

fish (NDArray) – The new fish to insert.

Given a new fish \(\mathbf{s} = (\mathbf{x}_c, \mathbf{n})\), we first compute the child octant index \(o\) for the fish based on its position relative to the center of this octant. If the child octant has not been initialized yet (i.e., there is no data presently there), then we initialize a new octree node for that octant with half the side length of this node and an offset center. Then, we insert the fish into that octant.

insert_data(fish)

Inserts a new fish into this octant.

Parameters:

fish (NDArray) – The fish to be inserted.

This follows the Barnes-Hut hierarchical tree generation algorithm as written from Ventimiglia and Wayne [5].

Constructing the Barnes-Hut tree

To construct the Barnes-Hut tree, insert the bodies one after another. To insert a body b into the tree rooted at node x, use the following recursive procedure:

1. If node x does not contain a body, put the new body b here.

2. If node x is an internal node, update the center-of-mass and total mass of x. Recursively insert the body b in the appropriate quadrant.

3. If node x is an external node, say containing a body named c, then there are two bodies b and c in the same region. Subdivide the region further by creating four children. Then, recursively insert both b and c into the appropriate quadrant(s). Since and c may still end up in the same quadrant, there may be several subdivisions during a single insertion. Finally, update the center-of-mass and total mass of x.

calculate_feature_positions()

Recursively computes the positions of the front and back for each cluster average within the greater Octree via finds.calculations.calculate_feature_positions().

compute_interaction(fish_pos, fish_front, fish_back, maximum_ratio)

Recursively computes the interaction of the fish at fish_pos with all of the fish (and clustered fish) within this tree.

Parameters:

• fish_pos (NDArray) – The position of the reference fish we’re computing the interaction with.

• fish_front (NDArray) – The position of the front of the reference fish.

• fish_back (NDArray) – The position of the back of the reference fish.

• minimum_ratio (float) – The Barnes-Hut Ratio \(\theta\), or the minimum ratio \(s/d < \theta\) of side length \(s\) to the distance from the reference fish to the center of the octant \(d\) to compute the interaction at this node rather than summing the interaction at the child nodes via recursion.

Returns:

A 1-dimensional array of 6 values containing the front and back interaction for the reference fish.

Return type:

NDArray

average

front_pos

back_pos

finds.calculations.build_octree(system)

Builds the Barnes-Hut Octree.

Parameters:

system (NDArray) – The system matrix.

Returns:

The built Octree.

Return type:

OctreeNode

This function builds an Octree out of a system matrix in three steps. First, it determines the longest distance along any given axis between two fish within the system as the maximum between the differences of the maximums and minimums on each axis.

For example, if we had the points

\[(1,1,3), (2,0,8), (3, 5, 4), (6, 2, 4), (11, 8, 3)\]

the minimums for each dimension would be \((1,0,3)\) and the maximums would be \((11,8,8)\), making the differences \((10,8,5)\) and the longest distance would be 10.

Then, this value (increased slightly by \(0.1\%\)) is used as the side length for the root cube of the Octree, and the position of the center is defined as half of the difference vector (which would be \((5,4,2.5)\) in the previous example).

Then, each of the fish in the system are inserted sequentially into the Octree using OctreeNode.insert_data(), and once all of the fish are inserted into the tree (and all of the nodes are clustered), all of the feature positions are calculated for the tree with OctreeNode.calculate_feature_positions().

finds.calculations.calculate_feature_interaction(feature_a_pos, feature_b_pos)

Computes the individual interaction vector between feature \(\alpha\) of fish \(i\) and feature \(\beta\) of fish \(j\), defined as the displacement between their positions divided by the cube of its norm:

(16)\[\mathbf{c}_{\alpha\beta} = \frac{\mathbf{r}_{\alpha\beta}} {r_{\alpha\beta}^3}\]

where

(17)\[\mathbf{r}_{\alpha\beta} = \mathbf{x}_{\alpha,i} - \mathbf{x}_{\beta,j}.\]

Parameters:

• feature_a_pos (NDArray) – The position of the first feature, \(\mathbf{x}_{\alpha,i}\).

• feature_b_pos (NDArray) – The position of the second feature, \(\mathbf{x}_{\beta,j}\).

Returns:

The interaction vector between features \(\alpha\) and \(\beta\), \(\mathbf{c}_{\alpha\beta}\).

Return type:

NDArray

finds.calculations.calculate_fish_interaction(fish_front, fish_back, other_front, other_back)

Returns the front and back interaction vectors between two fish per (8).

Return type:

NDArray

finds.calculations.compute_interaction_barnes_hut(system, bh_ratio, show_progress=False)

Computes the interaction vectors using the Barnes-Hut approximation as described in Barnes-Hut Clustering Approximation.

Parameters:

• system (NDArray) – The system.

• bh_ratio (float) – The maximum ratio \(\theta\) of partition size to particle distance for which to compute on clustered nodes.

• show_progress (bool) – Whether or not to show the calculation progress on each time-step.

Returns:

The array of interaction vectors.

Return type:

NDArray

finds.calculations.compute_interaction_pairwise(system)

Computes the sum of the pairwise interactions for each fish for both head and tail interactions.

Returns:

A matrix with shape \((N,3)\) consisting of three- dimensional rows, each with units \(m^{-2}\) encoding the sum of all of the interaction vectors.

Return type:

NDArray

Parameters:

system (NDArray) – The system.

finds.calculations.calculate_feature_velocities(system, use_barnes_hut, bh_ratio, show_progress=False)

Calculates the velocities for the head and tail of all fish in a system per equation (5).

Returns:

The first-derivative of the feature_positions matrix, or a matrix containing the velocities of all of the features in the system.

Return type:

NDArray

Parameters:

• system (NDArray) – The system.

• use_barnes_hut (bool) – Whether or not to simplify the calculations using the Barnes-Hut algorithm.

• bh_ratio (float) – The Barnes-Hut ratio.

• show_progress (bool) – Whether or not to show the calculation progress of this time-step

finds.calculations.calculate_system_derivative(system, use_barnes_hut, bh_ratio, show_progress=False)

Computes the derivative of the system matrix per equation :eq:fish_derivative:.

Returns:

The derivative of the system matrix.

Return type:

NDArray

Parameters:

• system (NDArray) – The system matrix.

• use_barnes_hut (bool) – Whether or not to use Barnes-Hut approximation to simplify the calculation.

• bh_ratio (float) – The Barnes-Hut Ratio

• show_progress (bool) – Whether or not to display the progress of this time-step