New report and code for final submission.HEADmaster

1 files changed, 168 insertions, 0 deletions

diff --git a/report/Chapters/simulation.tex b/report/Chapters/simulation.tex
new file mode 100644
index 0000000..ce530fc
--- /dev/null
+++ b/report/Chapters/simulation.tex
@@ -0,0 +1,168 @@
+% "Methods" section - describe basic structure of LD algo without optimizations
+\chapter{Simulation}
+
+\section{Simulation Setup}
+
+The physics were first developed as pseudo-code before implementation in FORTRAN. Particle mass is normalized to $m = 1$ throughout, initial velocities are sampled using the Box--Muller transform to produce Gaussian-distributed random numbers, and positions and velocities are integrated using the velocity-Verlet scheme.
+
+\subsection{Particle Initialization}
+Before a simulation is conducted, particles must be prescribed initial positions and velocities. Positions can be chosen arbitrarily with the exception being that they should not overlap. If two particles are placed sufficiently close, a strong interaction can occur yielding significant forces and an addition of energy to the system. It is also nonphysical to have particles overlap, as in the real world they have a finite volume preventing this. A scheme was introduced to place particles at randomly chosen discrete lattice points in lieu of a random continuous position and alleviate the relevant issues. A visual overview of the placement is shown in Figure \ref{fig:box-mesh}:
+
+\begin{figure}[H]
+  \centering
+  \resizebox{0.45\textwidth}{!}{\input{Figures/mesh.pdf_tex}}
+  \caption{A visual example of stochastic particle initialization using a fine mesh structure. A lattice is defined using the outer box geometry, and lattice points are picked at random for lattice placement. Each particle is placed on the discrete grid points such that particle-particle overlap is mitigated. A border region (10\%) is excluded around the domain boundary.}
+  \label{fig:box-mesh}
+\end{figure}
+
+While this technique removes a certain degree of `randomness,' it does not constitute playing an unfair game. The positions are still randomly sampled, they just truncated to discrete points using a ceiling function. The ceiling function was chosen specifically to avoid zero-indexing mesh positions in Fortran. The algorithm for placing points is presented in Listing \ref{alg:mesh-pts}.
+
+\begin{algorithm}[H]
+\caption{Initial Particle Positions Algorithm}
+\label{alg:mesh-pts}
+\begin{algorithmic}[1]
+\State $n_{mesh} = L/r_c$ \Comment{Define the number of mesh points in each direction}
+\State $d=0.9\times L/n_{mesh}$ \Comment{And the spacing between mesh points}
+\While{$k < N$} \Comment{Loop over all particles in the system}
+    \State $l_{x,y} \gets \mathrm{PNR} \in (0,1) $ \Comment{Generate 2 pseudorandom numbers}
+    \State $i_{x,y} = \mathrm{ceiling}(l_{x,y}\times n_{mesh})$ \Comment{Transform the random numbers to a lattice position}
+    \If{$\mathrm{empty}(i_x,i_y)$} \Comment{If the position does not have a particle, place one there}
+        \State $p_{x,y} = (i_{x,y}\times d) - L/2$ \Comment{Update the $k$-th particle to have the position}
+        \State $k = k + 1$ \Comment{Move on to the next particle}
+    \EndIf
+\EndWhile
+\end{algorithmic}
+\end{algorithm}
+
+Particle velocities were sampled from a Gaussian (normal) distribution. The Box-Muller transform was applied to generate pairs of normally-distributed numbers using uniformly-distributed numbers, which were then applied for prescribing an initial $x$- and $y$-velocity to each particle.
+%\note{algorithm 2 for initial particles positions for the particle initialization is showing in the 2.1.2 particle motion section. Fix | REPLY: Add a [h] after the begin\{\} statements to force objects like figs, tables, and algos to be placed '\textbf{h}ere.}
+
+\subsection{Particle Motion}
+Each particle experiences three forces at every time step: viscous drag, a stochastic thermal kick, and WCA repulsion from neighbouring particles. With $m = 1$ throughout, force and acceleration are interchangeable.
+
+The drag force is evaluated at the half-step velocity, consistent with the Verlet update order. Initial velocities are drawn from the Maxwell distribution using the Box--Muller transform, giving each component a Gaussian distribution with zero mean and variance proportional to the system temperature, consistent with thermal equilibrium. Positions and velocities are integrated using the velocity-Verlet scheme:
+
+\begin{algorithm}[H]
+\caption{Velocity Verlet Integration Loop}
+\begin{algorithmic}[1]
+\State Initialise particle positions, velocities, accelerations
+\While{$t < t_{\max}$}
+    \For{$i = 1$ to $N$}
+        \State $v_{hx} \gets v_x + \frac{1}{2} a_x \, dt$ \Comment{half-step velocity}
+        \State $x \gets x + v_{hx} \, dt$                  \Comment{update position}
+        \State Impose boundary conditions
+    \EndFor
+    \For{$i = 1$ to $N$}
+        \State Compute forces and accelerations at new position
+        \State $v_x \gets v_{hx} + \frac{1}{2} a_x \, dt$ \Comment{full-step velocity}
+        \State Record data
+    \EndFor
+    \State Compute ensemble averages
+    \State $t \gets t + dt$
+\EndWhile
+\end{algorithmic}
+\end{algorithm}
+
+The order of operations is important: positions must be updated before forces are computed, and half-step velocities must be used consistently throughout each step.
+
+\subsection{Boundary Conditions}
+
+The simulation domain is a square box $[-L/2,\,L/2]$ with reflective walls, particles bounce off boundaries rather than wrapping around (which would be periodic boundary conditions). The reflection rule for, e.g., the top wall is:
+
+\begin{algorithm}[H]
+\caption{Reflective boundary condition with consistency check}
+\label{alg:bc}
+\begin{algorithmic}[1]
+\If{$y(i) > L/2$}
+    \State $y(i) \gets L - y(i)$ \Comment{reflect position}
+    \State $vhy(i) \gets -vhy(i)$ \Comment{reverse normal velocity}
+\EndIf
+\end{algorithmic}
+\end{algorithm}
+
+The same logic applies to all four walls. An \texttt{outside} flag is set if a particle cannot be reflected back inside (e.g.\ after an extreme overshoot), in which case it is excluded from further updates. After every call to \texttt{impose\_BC(i)}, the consistency check verifies the particle is either inside the box or flagged as outside; if this never triggers, the BC logic is internally consistent.
+
+Three tests were used to validate the boundary condition implementation.
+
+\paragraph{Test 1: Logical consistency check.}
+The consistency check was embedded in Algorithm \ref{alg:bc} was run throughout all simulations. No BC errors were triggered in any test case.
+
+\paragraph{Test 2: Individual trajectories.}
+Trajectories of particles 1, 2, and 3 were logged to separate files and plotted to visually confirm smooth, physically reasonable motion with correct reflections at the walls.
+
+\begin{figure}[H]
+    \centering
+    \subcaptionbox{Working principle of a specular reflective boundary 
+    condition. Red: incident particle penetrates boundary. Magenta: 
+    position corrected as if bounced. Blue: particle continues in next 
+    timestep.\label{fig:reflective-BC}}[0.50\textwidth]
+        {\resizebox{0.45\textwidth}{!}{\input{Figures/reflectiveBC.pdf_tex}}}
+    \hfill
+    \subcaptionbox{Trajectories of particles 1 (in red), 2 (in green), and 3 (in blue) showing correct 
+    reflections at the boundary walls.\label{fig:trajectories}}[0.45\textwidth]
+        {\def\svgwith{0.40\textwidth}\input{Figures/trajectories.tex}}
+    \caption{Boundary condition behaviour.}
+    \label{fig:BC-combined}
+\end{figure}
+
+\paragraph{Test 3: Stress test with large time step.}
+Setting $dt = 1.0$ forces large overshoots to test whether \texttt{impose\_BC} handles extreme cases. No failure case was triggered, which is consistent with correct behaviour but not conclusive on its own — a particle can overshoot by more than $L$ in a single step and a single reflection may not be sufficient to bring it back inside.
+
+\subsection{Particle Interactions}
+
+Adding pairwise interactions increases the computational complexity from $\mathcal{O}(N)$ to $\mathcal{O}(N^2)$, since each particle must be compared with every other particle at each timestep. This pairwise force evaluation constitutes the dominant computational cost of the simulation.
+
+To ensure particles actually interact during testing, the particle density was increased so that particles are close enough to fall within $r_c$ of each other. Without this, the interaction code would never be exercised.
+
+The interaction loop follows this structure:
+
+\begin{algorithm}[H]
+\caption{Force calculation with Langevin dynamics}
+\begin{algorithmic}[1]
+\For{$i = 1$ to $N$}
+    \State Apply drag and random (Langevin) forces to particle $i$
+    \ForAll{$j \neq i$}
+        \State $d_{ij} \gets \lVert x_i - x_j \rVert$
+        \If{$d_{ij} < r_c$}
+            \State $F \gets -\dfrac{d}{dr}V(r)$
+            \State Apply $F$ to particle $i$ along $r_{ij}$
+        \EndIf
+    \EndFor
+\EndFor
+\end{algorithmic}
+\end{algorithm}
+
+%\note{I wonder for this section if we could just have the algorithm or the equations 2.1 and 2.2. they essentially say the same thing | REPLY: The rubric says that WCA should be mentioned in the intro. Maybe put the equations there and keep the algo here?}
+% done - addy
+
+An output file was added to log whenever an interaction was detected, recording the positions of the interacting pair. This allowed direct verification that interactions were occurring at the right locations and not spuriously triggered.
+
+\paragraph{Energy conservation check.}
+Kinetic energy vs.\ time was plotted to verify the thermostat maintains a stable average energy. An initial version showed energy not being conserved; this was diagnosed and fixed.
+\begin{figure}[H]
+  \centering
+  %\includegraphics[width=0.62\linewidth]{Figures/walltime_vs_N.pdf}
+  \input{Figures/kinetic_energy.tex}
+  \caption{Sum of System Kinetic Energy vs. Time. Simulated with $N=1000$ particles and in one case no stochastic/drag forces, normalized to start at equal energies. Note the dips in kinetic energy when interactions are taking place and kinetic energy is being transferred to potential energy in the form of particle-particle repulsion, which later return to nominal after interactions conclude. This is not seen in the run with stochastic forces enabled.}
+  \label{fig:kinetic-energy}
+\end{figure}
+%\note{does anyone have this figure? I know Lennaert made it while in class, but it might be cool to include if we have time - Addy}
+
+\section{Performance}
+\subsection{Simulation Complexity}
+
+\begin{figure}[H]
+  \centering
+  %\includegraphics[width=0.90\linewidth]{Figures/walltime_vs_N.pdf}
+  \input{Figures/timing.tex}
+  \caption{Tracked wall times vs number of particles}
+  \label{fig:Newwalltime_vs_N}
+\end{figure}
+
+Wall time grows with increasing $N$ number of particles, with underlying quadratic behaviour, consistent with the expected $\mathcal{O}(N^2)$ complexity. 
+
+\subsection{Module Performance}
+
+To identify computational bottlenecks, wall time was profiled by subroutine. The dominant cost is \textbf{pairwise interactions}: the double loop over all $i$--$j$ pairs evaluates $N(N-1)/2$ distances per time step, with each iteration requiring a square root and, when within $r_c$, several power-law evaluations for the WCA force. \textbf{Disk I/O} is a significant secondary cost, as positions and 
+velocities are written for every particle at every time step; reducing the output frequency is a straightforward way to recover performance. \textbf{Random number generation} accumulates measurable cost over many time steps due to two length-$N$ array draws per step. By contrast, \textbf{boundary condition} checks and \textbf{velocity--Verlet} integration are both $\mathcal{O}(N)$ and contribute 
+minimally to total runtime. These bottlenecks motivate the optimizations described in Chapter~\ref{ch:Optimization}.
\ No newline at end of file