diff --git a/src/FluidSimulation.cpp b/src/FluidSimulation.cpp
index e1fec92dee890635bca69c597ac3b9f4aa1657c3..09752c8cd1a5ef031a92f3ae7c35285450115cc4 100644
--- a/src/FluidSimulation.cpp
+++ b/src/FluidSimulation.cpp
@@ -17,6 +17,7 @@ namespace FluidSimulation {
         m_viscosity = h_VISCOSITY;
 
         m_simulationStepsPerFrame = 2;
+        m_mode = SimulationMode::OneLargeCube;
     }
 
     FluidSimulation::~FluidSimulation() {
@@ -32,7 +33,7 @@ namespace FluidSimulation {
         std::vector<float3> h_positions(m_numParticles);
         std::vector<float3> h_velocities(m_numParticles);
         float gridSpacing = 0.5 * m_smoothingRadius;
-        initParticles(h_positions, h_velocities, m_numParticles, gridSpacing);
+        initParticles(h_positions, h_velocities, m_numParticles, gridSpacing, m_mode);
     
         FluidSimulationGPU::resetSimulation(
             h_positions.data(), h_velocities.data(), m_numParticles);
@@ -47,6 +48,12 @@ namespace FluidSimulation {
         updateConstants();
     }
 
+    void FluidSimulation::setSimulationMode(SimulationMode newMode)
+    {
+        m_mode = newMode;
+        Reset();
+    }
+
     void FluidSimulation::updateConstants() {
     // Send updated parameters to CUDA
         FluidSimulationGPU::updateConstantsOnDevice(
@@ -61,7 +68,7 @@ namespace FluidSimulation {
         std::vector<float3> h_positions(m_numParticles);
         std::vector<float3> h_velocities(m_numParticles);
         float gridSpacing = 0.5 * m_smoothingRadius;
-        initParticles(h_positions, h_velocities, m_numParticles, gridSpacing);
+        initParticles(h_positions, h_velocities, m_numParticles, gridSpacing, m_mode);
 
         // Initialize particle device arrays
         FluidSimulationGPU::initializeFluidDeviceArrays(h_positions.data(), h_velocities.data(), m_numParticles);
diff --git a/src/FluidSimulationApp.cpp b/src/FluidSimulationApp.cpp
index 62c9fb8e63418fa43cc4b2e83f1ccaec85e40fae..cb6f9fab76e3469e05a819d45e6672da9f0cd1f3 100644
--- a/src/FluidSimulationApp.cpp
+++ b/src/FluidSimulationApp.cpp
@@ -3,6 +3,7 @@
 #include <iostream>
 #include "imgui.h"
 #include "settings/constants.h"
+#include "include/SimulationMode.h"
 
 namespace FluidSimulation {
 
@@ -60,7 +61,7 @@ namespace FluidSimulation {
 
             // Particle Mass
             ImGui::Text("Particle Mass");
-            ImGui::SliderFloat("##ParticleMass", &particleMass, 0.01f, 0.1f, "%.3f");
+            ImGui::SliderFloat("##ParticleMass", &particleMass, 0.01f, 1.4f, "%.3f");
 
             // Rest Density
             static int restDensityInt = static_cast<int>(restDensity);
@@ -141,12 +142,25 @@ namespace FluidSimulation {
             m_fluidSimulation->Reset();
         }
 
+        ImGui::PushItemWidth(-1); // Make sliders take up the full available width
+        
+        ImGui::Separator();
+
+        static int selectedOption = static_cast<int>(m_fluidSimulation->getSimulationMode()); // Sync with current mode
+
+        const char* options[] = { "One Large Cube", "Two Small Cubes" };
+
+        if (ImGui::Combo("##DropDownMenu", &selectedOption, options, IM_ARRAYSIZE(options)))
+        {
+            m_fluidSimulation->setSimulationMode(static_cast<SimulationMode>(selectedOption));
+        }
+
+
         ImGui::Separator();
 
         static int stepsPerFrame = m_fluidSimulation->getSimulationStepsPerFrame();
 
         // Full-width sliders
-        ImGui::PushItemWidth(-1); // Make sliders take up the full available width
 
         // Particle Mass
         ImGui::Text("Simulation steps/frame");
diff --git a/src/include/FluidSimulation.h b/src/include/FluidSimulation.h
index a63ba6a7483ddd75e00b30456fa5178a12405605..b790faaf21fde29ab60a721cf594b822891a0df7 100644
--- a/src/include/FluidSimulation.h
+++ b/src/include/FluidSimulation.h
@@ -3,6 +3,7 @@
 #include <vector>
 #include <cuda_runtime.h>
 #include <GL/glew.h>
+#include "SimulationMode.h"
 
 namespace FluidSimulation {
 
@@ -22,6 +23,8 @@ namespace FluidSimulation {
         float m_smoothingRadius;
         int m_simulationStepsPerFrame;
 
+        SimulationMode m_mode;
+
     public:
         FluidSimulation(int numParticles);
         ~FluidSimulation();
@@ -44,13 +47,15 @@ namespace FluidSimulation {
         float getSmoothingRadius() { return m_smoothingRadius; }
 
         int getSimulationStepsPerFrame() { return m_simulationStepsPerFrame; };
-        
+        SimulationMode getSimulationMode() { return m_mode; };
+
         void setParticleMass(float particleMass) { m_particleMass = particleMass; };
         void setRestDensity(float restDensity) { m_restDensity = restDensity; };
         void setViscosity(float visocisty) { m_viscosity = visocisty; };
         void setSmoothingRadius(float smoothRadius) { m_smoothingRadius = smoothRadius; };
     
         void setSimulationStepsPerFrame(int steps) { m_simulationStepsPerFrame = steps; };
+        void setSimulationMode(SimulationMode newMode);
 
         void ResetFluidParameters();
     };
diff --git a/src/include/SimulationMode.h b/src/include/SimulationMode.h
new file mode 100644
index 0000000000000000000000000000000000000000..f41ced573327a8e3e81312832193fd8440a80e85
--- /dev/null
+++ b/src/include/SimulationMode.h
@@ -0,0 +1,13 @@
+#ifndef SIMULATION_MODE_H
+#define SIMULATION_MODE_H
+
+namespace FluidSimulation
+{
+    enum class SimulationMode
+    {
+        OneLargeCube,
+        TwoSmallCubes
+    };
+}
+
+#endif // SIMULATION_MODE_H
diff --git a/src/include/spawn_particles.h b/src/include/spawn_particles.h
index a27939e3dd27ffb61d4aa2f3aa547ded1df6821e..b968e8d489d03784ee3ecf5e7049632f880d35e0 100644
--- a/src/include/spawn_particles.h
+++ b/src/include/spawn_particles.h
@@ -4,6 +4,7 @@
 
 #include <vector>
 #include <cuda_runtime.h> // for float3
+#include "SimulationMode.h"
 
 namespace FluidSimulation {
 
@@ -18,7 +19,8 @@ void initParticles(
     std::vector<float3> &h_positions,
     std::vector<float3> &h_velocities,
     const int NUM_PARTICLES,
-    const float GRID_SPACING
+    const float GRID_SPACING,
+    SimulationMode mode
 );
 }
 #endif // SPAWN_PARTICLES_H
\ No newline at end of file
diff --git a/src/spawn_particles.cpp b/src/spawn_particles.cpp
index 59986b0b7994b2420699443313781a82db6dbc58..9f6672cc496225281d8b1a52023d21df3560e138 100644
--- a/src/spawn_particles.cpp
+++ b/src/spawn_particles.cpp
@@ -13,75 +13,61 @@ void initParticles(
     std::vector<float3> &h_positions,
     std::vector<float3> &h_velocities,
     const int NUM_PARTICLES,
-    const float GRID_SPACING)
+    const float GRID_SPACING,
+    SimulationMode mode
+)
 {
     // Seed the random number generator
     srand(static_cast<unsigned int>(time(0)));
 
-    // Calculate particles per axis
-    int particlesPerAxis = static_cast<int>(ceil(pow(NUM_PARTICLES, 1.0f / 3.0f)));
-
-    int index = 0;
-    for (int i = 0; i < particlesPerAxis && index < NUM_PARTICLES; ++i) {
-        for (int j = 0; j < particlesPerAxis && index < NUM_PARTICLES; ++j) {
-            for (int k = 0; k < particlesPerAxis && index < NUM_PARTICLES; ++k) {
-                // Center particles around the origin
-                float x = (i - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING;
-                float y = (j - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 4.0f;
-                float z = (k - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 0.5;
-
-                // Add random perturbation within [-0.1*dx, 0.1*dx]
-                float perturbation = 0.1f * GRID_SPACING;
-                float rand_x = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-                float rand_y = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-                float rand_z = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-
-                h_positions[index] = make_float3(
-                    x + rand_x,
-                    y + rand_y,
-                    z + rand_z
-                );
-
-                // Initialize velocities
-                h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f);
-
-                ++index;
+    auto initializeCube = [&](int startIndex, int particlesPerAxis, float offsetX, float offsetY, float offsetZ) {
+        int index = startIndex;
+        for (int i = 0; i < particlesPerAxis && index < NUM_PARTICLES; ++i) {
+            for (int j = 0; j < particlesPerAxis && index < NUM_PARTICLES; ++j) {
+                for (int k = 0; k < particlesPerAxis && index < NUM_PARTICLES; ++k) {
+                    float x = (i - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + offsetX;
+                    float y = (j - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + offsetY;
+                    float z = (k - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + offsetZ;
+
+                    float perturbation = 0.1f * GRID_SPACING;
+                    float rand_x = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+                    float rand_y = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+                    float rand_z = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
+
+                    h_positions[index] = make_float3(x + rand_x, y + rand_y, z + rand_z);
+                    h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f);
+                    ++index;
+                }
             }
         }
-    }
-    /*
-    for (int i = 0; i < particlesPerAxis && index < NUM_PARTICLES; ++i) {
-        for (int j = 0; j < particlesPerAxis && index < NUM_PARTICLES; ++j) {
-            for (int k = 0; k < particlesPerAxis && index < NUM_PARTICLES; ++k) {
-                // Center particles around the origin
-                float x = (i - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 9;
-                float y = (j - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING + 2.0f;
-                float z = (k - (particlesPerAxis - 1) / 2.0f) * GRID_SPACING - 0.5;
-
-                // Add random perturbation within [-0.1*dx, 0.1*dx]
-                float perturbation = 0.1f * GRID_SPACING;
-                float rand_x = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-                float rand_y = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-                float rand_z = ((float)rand() / RAND_MAX - 0.5f) * 2.0f * perturbation;
-
-                h_positions[index] = make_float3(
-                    x + rand_x,
-                    y + rand_y,
-                    z + rand_z
-                );
-
-                // Initialize velocities
-                h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f);
-
-                ++index;
-            }
+        return index; // Return the next index
+    };
+
+    int particlesPerAxis = 0;
+    int index = 0;
+
+    switch (mode)
+    {
+        case SimulationMode::OneLargeCube:
+            particlesPerAxis = static_cast<int>(ceil(pow(NUM_PARTICLES, 1.0f / 3.0f)));
+            index = initializeCube(0, particlesPerAxis, 0.0f, 4.0f, 0.5f);
+            break;
+
+        case SimulationMode::TwoSmallCubes:
+        {
+            int halfParticles = NUM_PARTICLES / 2; // Split particles evenly between two cubes
+            particlesPerAxis = static_cast<int>(ceil(pow(halfParticles, 1.0f / 3.0f)));
+
+            index = initializeCube(0, particlesPerAxis, -4.0f, 4.0f, 1.0f);
+            index = initializeCube(index, particlesPerAxis, 8.0f, 4.0f, -1.0f);
+            break;
         }
     }
-    */
-    // If NUM_PARTICLES is not a perfect cube, initialize remaining particles
+
+    // Initialize remaining particles, if any
     for (; index < NUM_PARTICLES; ++index) {
-        h_positions[index] = make_float3(0.0f, 10.0f, 0.0f); // Default position
-        h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f); // Default velocity
+        h_positions[index] = make_float3(0.0f, 10.0f, 0.0f);
+        h_velocities[index] = make_float3(0.0f, 0.0f, 0.0f);
     }
 }