\chapter{Algorithms}
\label{AppAlgos}
\section{Domain Decomposition Algorithms}

\subsubsection{Updating Particle Sector}

The \code{updateParticleSector(x,y)} function outlined in algorithm \ref{updateParticleSector_code} is called for each particle after their position is updated in the movement loop. Taking the particles position coordinates as input, the function calculates the sector coordinates (\code{sectX, sectY})for that position and the corresponding sector index using equations \eqref{pos-to-sectCoords_DD} and \eqref{pos-to-index_DD} respectively. A validation check is also included in the function where particles at invalid positions are flagged by assigning them the sector index $-1$. As particles are assigned to their respective sectors, we also keep track of sector populations, which are used in the subsequent \code{updateSpatialLookup()} function for sorting the particles. 


\subsubsection{Particle Sorting and Spatial Lookup}
The spatial lookup update algorithm outlined in \ref{updateSpatialLookup_code} shows the procedure for preparing the optimized spatial lookup. First, the populations of each sector are used to determine the size of the block needed for a given sector. Following the first sector which always starts at index 0, the start index of each subsequent sector block is calculated as the start index of the previous block plus its population. This first loop yields the \code{sectorStartIndices[M]} array where the $i^{th}$ entry is the index in the sorted particle list where the block for sector $i$ begins. Using the start indices array, particles can then be directly inserted into the sorted list based on their respective sector indices. As particles are placed into the list, the population of inserted sector members is tracked so that no two particles are placed at the same index.

While typical sorting methods like Quicksort and Merge Sort with complexity of $\O{N\log N}$ would offer a performance boost over the $\O{N^2}$ scaling of the original implementation, using the apriori knowledge about the sector populations allows us to produce the sorted particle list in only $\O{N}$ time. Thus, a significant boost to performance is expected from the implementation of this method. 

\subsubsection{Neighbourhood Search}
The final step of domain decomposition in the simulation loop is to make use of the sorted list and start indices produced from the spatial lookup. In algorithm \ref{updateNeighSearch_code}, the procedure for finding interacting neighbours is outlined. The outermost loop iterates over each sector $s$, first determining the start and end indices of $s$. If the start and end indices for a given sector are the same, it has a population of $0$ and is skipped by empty sector check. If the sector is populated, the first inner loop is entered which iterated over the particles of sector $s$. For each of these particles $p1$, the next For-loop goes through each of the (up to) nine surrounding sectors and identifies neighbouring particles $p2$ in those sectors. The separation distance between $p1$ and each $p2$ is then evaluated and the corresponding LJ force is applied to $p1$. Following the loop, each particle has its net deterministic force $F$, which is then used in the Velocity Verlet acceleration update. 

\begin{algorithm}[h]
\caption{updateParticleSector Function}
\label{updateParticleSector_code}
\begin{algorithmic}
    \State sectorPops[M] = [0,0,...,0] \Comment{Array to store sector populations}
    \Function{updateParticleSector}{$x$,$y$}
    \State sectX $\gets$ floor$((x + L/2)/(l_s))$
    \State sectY $\gets S_y -$ floor$((x + L/2)/(l_s))$
    \If{sectX $> S_x$ OR sectY $> S_y$} \Comment{Set out-of-bounds particles to null sector}
        \State sectorIndex $\gets -1$
    \Else 
        \State sectorIndex $\gets sectX + sectY*S_y$ \Comment{Assign sector index to particle}
        \State sectorPops[sectInd] $+= 1$ \Comment{Add one to sector population}
    \EndIf
    \EndFunction
\end{algorithmic}
\end{algorithm}


\begin{algorithm}[h]
\caption{updateSpatialLookup Function}
\label{updateSpatialLookup_code}
\begin{algorithmic}
    \State particles[N] \Comment{Incoming list of particles}
    \State sortedParticles[N] \Comment{An empty array to be filled with sorted particles}
    \State sectorPops[M] \Comment{Array of sector populations}
    \State sectorStartIndices[M] = [-1, -1, ... , -1] \Comment{Start index array initialized to -1}
    \Function{updateSpatialLookup}{particles[N]}
    \State Determine where each sector will start in the particle list:
    \For{$i = 1$ to $M$} \Comment{Loop over each sector}
        \If{$i==0$} \Comment{Special case for first sector: always starts at index 0}
            \State sectorStartIndices[i] $\gets 0$
        \Else
            \State sectorStartIndices[i] $\gets$ sectorStartIndices[i-1] + sectorPops[i-1]
            \State sectorPops[i-1] $\gets 0$ \Comment{Set sector pop to zero to recount during insertion step} 
        \EndIf
    \EndFor

    \State Insert particles into sorted list:
    \For{$p=1$ to $N$}
        \State $i \gets$ particle.sectorIndex
        \State insertIndex $\gets$ sectorStartIndices[i] + sectorPops[i]
        \State sectorPops[i] $+= 1$
    \EndFor

    \State particles[N] $\gets$ sortedParticles[N] \Comment{Copy sorted list into the original list}
    \EndFunction
\end{algorithmic}
\end{algorithm}



\begin{algorithm}[h]
\caption{updateNeighbourInteraction Function}
\label{updateNeighSearch_code}
\begin{algorithmic}
    \State 
    \Function{updateNeighbourInteraction}{particles[N], sectors[M]}
    \For{$s = 1$ to $M$} \Comment{Loop over sectors}
        \State sectorStart $\gets$ sectorStartIndices[s]
        \State sectorEnd $\gets$ sectorStart + sectorPops[s]
        \For{$p1 =$ sectorStart to sectorEnd} \Comment{Loop over particles in sector $s$}
            \For{$i =1$ to $9$} \Comment{Loop over neighbouring sectors}
                \State searchStart $\gets$ sectorStartIndices[i]
                \State searchEnd $\gets$ searchStart + sectorPops[i]
                \For{$p2 =$ searchStart to searchEnd} \Comment{Loop over particles in neighbouring sector $i$}
                    \If{$p1 == p1$}
                        continue \Comment{skip self-interaction}

                    \Else
                        \If{$||p1 - p2|| \leq r_{cut}$} \Comment{If separation distance is less than cutoff}
                            \State $F_{p1} +=$ applyLJ($p1,p2$) \Comment{Add LJ force to p1}
                        \EndIf
                    \EndIf
                \EndFor
            \EndFor
        \EndFor
    \EndFor  
    
    \EndFunction
\end{algorithmic}
\end{algorithm}