Collision bubble simulation

.github/workflows/ci.yaml

@@ -14,13 +14,15 @@ env:
 
 jobs:
   build-and-run:
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-22.04
     steps:
       - uses: actions/checkout@v3
       - name: install deps
         run: |
+          sudo apt install -y libpocl2
+          sudo apt install -y mesa-opencl-icd
           sudo apt install -y gcc g++ clang
-          sudo apt install -y opencl-c-headers opencl-clhpp-headers ocl-icd-opencl-dev
+          sudo apt install -y opencl-c-headers opencl-clhpp-headers ocl-icd-opencl-dev opencl-headers
           sudo apt install -y pocl-opencl-icd
           sudo apt install -y oclgrind
 
@@ -30,8 +32,15 @@ jobs:
       - name: run
         run: |
           mkdir data
-          POCL_DEVICES=pthread ./build/bubbleSim.exe $TEST_CONFIG $KERNEL_CL &> log_pthread.txt
+          echo "Starting simulations..."
+          POCL_DEVICES=basic ./build/bubbleSim.exe $TEST_CONFIG $KERNEL_CL
+
           POCL_DEVICES=basic ./build/bubbleSim.exe $TEST_CONFIG $KERNEL_CL &> log_basic.txt
+          echo "Finished basic..."
+          cat log_basic.txt
+          POCL_DEVICES=pthread ./build/bubbleSim.exe $TEST_CONFIG $KERNEL_CL &> log_pthread.txt
+          cat log_pthread.txt
+          echo "Finished pthread..."
 
           diff log_basic.txt log_pthread.txt > diff1.txt || true
           export NUMDIFF=$(wc -l < diff1.txt)

.gitignore

@@ -7,6 +7,9 @@
 dependencies/OpenCL
 bubbleSim/data
 
+# External tool files
+*.cltracer
+
 
 # User-specific files
 *.rsuser
@@ -17,6 +20,9 @@ bubbleSim/data
 
 # User-specific files (MonoDevelop/Xamarin Studio)
 *.userprefs
+tools
+traces
+
 
 # Mono auto generated files
 mono_crash.*

Makefile

@@ -12,7 +12,7 @@ ifndef OPTFLAGS
 endif
 
 #link against common OCL headers
-CXXFLAGS=-std=c++2a -I./dependencies/json/ $(OPTFLAGS)
+CXXFLAGS=-std=c++2a -I./dependencies/json/ -I./dependencies/ankerl/ $(OPTFLAGS)
 LDFLAGS=-lOpenCL -lm -lstdc++ $(OPTFLAGS)
 
 # The final build step of the executable

bubbleSim.sln

@@ -7,12 +7,18 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "bubbleSim", "bubbleSim\bubb
 EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		bubbleSim test|x64 = bubbleSim test|x64
+		bubbleSim test|x86 = bubbleSim test|x86
 		Debug|x64 = Debug|x64
 		Debug|x86 = Debug|x86
 		Release|x64 = Release|x64
 		Release|x86 = Release|x86
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{82AF250D-0A19-4018-B80B-85345DF10DF0}.bubbleSim test|x64.ActiveCfg = bubbleSim test|x64
+		{82AF250D-0A19-4018-B80B-85345DF10DF0}.bubbleSim test|x64.Build.0 = bubbleSim test|x64
+		{82AF250D-0A19-4018-B80B-85345DF10DF0}.bubbleSim test|x86.ActiveCfg = bubbleSim test|Win32
+		{82AF250D-0A19-4018-B80B-85345DF10DF0}.bubbleSim test|x86.Build.0 = bubbleSim test|Win32
 		{82AF250D-0A19-4018-B80B-85345DF10DF0}.Debug|x64.ActiveCfg = Debug|x64
 		{82AF250D-0A19-4018-B80B-85345DF10DF0}.Debug|x64.Build.0 = Debug|x64
 		{82AF250D-0A19-4018-B80B-85345DF10DF0}.Debug|x86.ActiveCfg = Debug|Win32

bubbleSim/base.h

@@ -1,22 +1,110 @@
 #pragma once
 #define _USE_MATH_DEFINES
 
+#define CL_MINIMUM_OPENCL_VERSION 120
 #define CL_TARGET_OPENCL_VERSION 120
 #define CL_HPP_MINIMUM_OPENCL_VERSION 120
 #define CL_HPP_TARGET_OPENCL_VERSION 120
 
-#include <CL/cl.hpp>
+#if __has_include("CL/opencl.hpp")
+# include <CL/opencl.hpp>
+#elif __has_include("CL/cl2.hpp")
+# include <CL/cl2.hpp>
+#else
+# include <CL/cl.hpp>
+#endif
+
 #include <array>
+#include <chrono>
 #include <cmath>
 #include <cstdio>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
+#include <map>
 #include <sstream>
+#include <unordered_map>
+#include <unordered_set>
 #include <vector>
 
+#include "ankerl/ankerl_map.hpp"
+
+// #define TIME_COLLIDE_DEBUG
 // #define LOG_DEBUG
+// #define DEBUG_OPENCL_KERNEL_RUNTIME_PROFILE
 
 typedef cl_double cl_numType;
 typedef double numType;
-typedef unsigned int u_int;
+typedef unsigned int u_int;
+
+enum SimulationFlags : std::uint32_t {
+  BUBBLE_ON = 0b00000000000000000000000000000001,
+  TRUE_VACUUM_BUBBLE_ON = 0b00000000000000000000000000000010,
+  BUBBLE_INTERACTION_ON = 0b00000000000000000000000000000100,
+  COLLISION_ON = 0b00000000000000000000000000001000,
+  SIMULATION_BOUNDARY_ON = 0b00000000000000000000000000010000,
+  COLLISION_MASS_STATE_ON = 0b00000000000000000000000000100000,
+};
+
+inline SimulationFlags operator|(SimulationFlags a, SimulationFlags b) {
+  return static_cast<SimulationFlags>(static_cast<std::uint32_t>(a) |
+                                      static_cast<std::uint32_t>(b));
+}
+
+enum StreamFlags : std::uint32_t {
+  STREAM_ON = 0b00000000000000000000000000000001,
+  STREAM_TIMESERIES = 0b00000000000000000000000000000010,
+  STREAM_PROFILE = 0b00000000000000000000000000000100,
+  STREAM_MOMENTUM = 0b00000000000000000000000000001000,
+  STREAM_MOMENTUM_PROFILE = 0b00000000000000000000000000010000,
+};
+
+inline StreamFlags operator|(StreamFlags a, StreamFlags b) {
+  return static_cast<StreamFlags>(static_cast<std::uint32_t>(a) |
+                                  static_cast<std::uint32_t>(b));
+}
+
+enum SimulationBufferFlags : std::uint64_t {
+  PARTICLE_X_BUFFER = (std::uint64_t)1 << 0,
+  PARTICLE_Y_BUFFER = (std::uint64_t)1 << 1,
+  PARTICLE_Z_BUFFER = (std::uint64_t)1 << 2,
+  PARTICLE_E_BUFFER = (std::uint64_t)1 << 3,
+  PARTICLE_PX_BUFFER = (std::uint64_t)1 << 4,
+  PARTICLE_PY_BUFFER = (std::uint64_t)1 << 5,
+  PARTICLE_PZ_BUFFER = (std::uint64_t)1 << 6,
+  PARTICLE_M_BUFFER = (std::uint64_t)1 << 7,
+  PARTICLE_COLLIDE_BUFFER = (std::uint64_t)1 << 8,
+  PARTICLE_IN_BUBBLE_BUFFER = (std::uint64_t)1 << 9,
+  PARTICLE_COLLISION_CELL_IDX_BUFFER = (std::uint64_t)1 << 10,
+  PARTICLE_dP_BUFFER = (std::uint64_t)1 << 11,
+  PARTICLE_INTERACTED_FALSE_BUFFER = (std::uint64_t)1 << 12,
+  PARTICLE_INTERACTED_TRUE_BUFFER = (std::uint64_t)1 << 13,
+  PARTICLE_PASSED_FALSE_BUFFER = (std::uint64_t)1 << 14,
+  CELL_THETA_ROTATION_BUFFER = (std::uint64_t)1 << 15,
+  CELL_THETA_AXIS_BUFFER = (std::uint64_t)1 << 16,
+  CELL_PHI_AXIS_BUFFER = (std::uint64_t)1 << 17,
+  CELL_E_BUFFER = (std::uint64_t)1 << 18,
+  CELL_LOGE_BUFFER = (std::uint64_t)1 << 19,
+  CELL_PX_BUFFER = (std::uint64_t)1 << 20,
+  CELL_PY_BUFFER = (std::uint64_t)1 << 21,
+  CELL_PZ_BUFFER = (std::uint64_t)1 << 22,
+  CELL_COLLIDE_BUFFER = (std::uint64_t)1 << 23,
+  CELL_PARTICLE_COUNT_BUFFER = (std::uint64_t)1 << 24,
+  CELL_LENGTH_BUFFER = (std::uint64_t)1 << 25,
+  CELL_COUNT_IN_ONE_AXIS_BUFFER = (std::uint64_t)1 << 26,
+  CELL_SHIFT_VECTOR_BUFFER = (std::uint64_t)1 << 27,
+  CELL_SEED_INT64_BUFFER = (std::uint64_t)1 << 28,
+  CELL_NO_COLLISION_PROBABILITY_BUFFER = (std::uint64_t)1 << 29,
+  BUBBLE_BUFFER = (std::uint64_t)1 << 30,
+  SIMULATION_DT_BUFFER = (std::uint64_t)1 << 31,
+  SIMULATION_BOUNDARY_BUFFER = (std::uint64_t)1 << 32,
+  SIMULATION_MASS_IN_BUFFER = (std::uint64_t)1 << 33,
+  SIMULATION_MASS_OUT_BUFFER = (std::uint64_t)1 << 34,
+  SIMULATION_DELTA_MASS_BUFFER = (std::uint64_t)1 << 35
+};
+
+inline SimulationBufferFlags operator|(SimulationBufferFlags a,
+                                       SimulationBufferFlags b) {
+  return static_cast<SimulationBufferFlags>(static_cast<std::uint64_t>(a) |
+                                            static_cast<std::uint64_t>(b));
+}

bubbleSim/bubble.cpp

@@ -1,7 +1,45 @@
 #include "bubble.h"
 
+PhaseBubble::PhaseBubble(numType t_initialRadius, numType t_initialSpeed,
+                         numType t_dV, numType t_sigma) {
+  /*
+  cl_double radius;
+   cl_double radius2;  // Squared
+   cl_double radiusAfterStep2;  // (radius + speed * dt)^2
+
+   cl_double speed;
+
+   cl_double gamma;
+   cl_double gammaXspeed;  // gamma * speed
+  */
+  // m_gamma = 1/sqrt(1-v^2) = 1/exp(log1p(-v^2)*0.5)
+  numType radius2 = std::pow(t_initialRadius, 2);
+  numType gamma =
+      1.0 / std::exp((std::log1p((-t_initialSpeed * t_initialSpeed)) * 0.5));
+  numType gammaXspeed = t_initialSpeed * gamma;
+  m_bubble =
+      Bubble{t_initialRadius, radius2, t_initialSpeed, gamma, gammaXspeed};
+
+  m_dV = t_dV;
+  m_sigma = t_sigma;
+  m_initialRadius = t_initialRadius;
+  m_energy = calculateEnergy();
+  m_initialEnergy = calculateEnergy();
+
+  if (m_sigma < 0) {
+    std::cerr << "sigma < 0" << std::endl;
+    std::terminate();
+  }
+
+  // m_area = 4 * M_PI * std::pow(t_initialRadius, 2);
+  // m_volume = 4 * (pow(t_initialRadius, 3) * M_PI) / 3;
+  // m_energy = calculateEnergy();
+}
+
 PhaseBubble::PhaseBubble(numType t_initialRadius, numType t_initialSpeed,
                          numType t_dV, numType t_sigma,
+                         numType t_critical_radius,
+                         std::uint64_t& t_buffer_flags,
                          cl::Context& cl_context) {
   /*
   cl_double radius;
@@ -20,15 +58,23 @@ PhaseBubble::PhaseBubble(numType t_initialRadius, numType t_initialSpeed,
   numType gamma =
       1.0 / std::exp((std::log1p((-t_initialSpeed * t_initialSpeed)) * 0.5));
   numType gammaXspeed = t_initialSpeed * gamma;
-  m_bubble = Bubble{t_initialRadius, radius2,     radius2, t_initialSpeed,
-                    gamma,           gammaXspeed, gamma};
+
+  m_bubble =
+      Bubble{t_initialRadius, radius2, t_initialSpeed, gamma, gammaXspeed};
   m_bubble_copy = m_bubble;
   m_bubbleBuffer =
       cl::Buffer(cl_context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                  sizeof(Bubble), &m_bubble, &openCLerrNum);
+  t_buffer_flags |= BUBBLE_BUFFER;
+
   m_dV = t_dV;
   m_sigma = t_sigma;
+
+  m_energy = calculateEnergy();
+
   m_initialRadius = t_initialRadius;
+  m_initialEnergy = calculateEnergy();
+  m_criticalRadius = t_critical_radius;
 
   if (m_sigma < 0) {
     std::cerr << "sigma < 0" << std::endl;
@@ -48,48 +94,89 @@ numType PhaseBubble::calculateVolume() {
 }
 
 numType PhaseBubble::calculateEnergy() {
-  return calculateArea() * m_sigma / sqrt(1 - m_bubble.speed * m_bubble.speed) -
-         calculateVolume() * m_dV;
-}
-numType PhaseBubble::calculateRadiusAfterStep2(numType dt) {
-  m_bubble.radiusAfterStep2 =
-      pow(std::fma(m_bubble.speed, dt, m_bubble.radius), 2);
-  return m_bubble.radiusAfterStep2;
+  return calculateArea() * m_sigma * m_bubble.gamma - calculateVolume() * m_dV;
 }
 
 void PhaseBubble::evolveWall(numType dt, numType dP) {
-  numType newRadius;  // R_(i+1)
-  numType newSpeed;   // V_(i+1)
-  numType newGamma;   // gamma_(i+1)
-  numType gammaChange;
+  numType new_radius;  // R_(i+1)
+  numType new_speed;   // V_(i+1)
+  numType new_gamma;   // gamma_(i+1)
+  numType gamma_change;
   numType sgn = ((0 < m_bubble.speed) - (m_bubble.speed < 0));  // sign of speed
 
-  if (m_bubble.gamma >= 10) {
-    newRadius = m_bubble.radius + dt * m_bubble.speed;
+  new_radius = m_bubble.radius + dt * m_bubble.speed;
 
-    gammaChange = (std::fma(m_dV, dt, dP) / m_sigma *
-                       std::sqrt((m_bubble.gamma - 1) / m_bubble.gamma) -
-                   2 * std::sqrt((m_bubble.gamma - 1) * m_bubble.gamma) /
-                       m_bubble.radius * dt) *
-                  sgn;
-    newGamma = m_bubble.gamma + gammaChange;
-    newSpeed = std::sqrt(1 - 1 / std::pow(newGamma, 2)) * sgn;
+  if (m_bubble.gamma >= 20 && false) {
+    gamma_change = (std::fma(m_dV, dt, -dP) / m_sigma *
+                        std::sqrt((m_bubble.gamma - 1) / m_bubble.gamma) -
+                    2 * std::sqrt((m_bubble.gamma - 1) * m_bubble.gamma) /
+                        m_bubble.radius * dt) *
+                   sgn;
 
+    new_gamma = m_bubble.gamma + gamma_change;
+    new_speed = std::sqrt(1 - 1 / std::pow(new_gamma, 2.)) * sgn;
   } else {
-    newRadius = m_bubble.radius + dt * m_bubble.speed;
+    numType velocityElement = 1 - std::pow(m_bubble.speed, 2.);
 
-    numType velocityElement = std::fma(-m_bubble.speed, m_bubble.speed, 1);
+    new_speed = m_bubble.speed +
+                pow(velocityElement, 1.5) * std::fma(m_dV, dt, -dP) / m_sigma -
+                2 * velocityElement / m_bubble.radius * dt;
 
-    newSpeed =
-        m_bubble.speed +
-        std::sqrt(pow(velocityElement, 3)) * std::fma(m_dV, dt, dP) / m_sigma -
-        2 * velocityElement * dt / m_bubble.radius;
-    newGamma = 1.0 / std::exp((std::log1p((-newSpeed * newSpeed)) * 0.5));
+    new_gamma = 1.0 / std::sqrt(1. - std::pow(new_speed, 2.));
   }
-  m_bubble.radius = newRadius;
-  m_bubble.speed = newSpeed;
-  m_bubble.gamma = newGamma;
 
+  m_bubble.radius = new_radius;
+  m_bubble.speed = new_speed;
+  m_bubble.gamma = new_gamma;
+  m_bubble.radius2 = std::pow(new_radius, 2);
+  m_bubble.gammaXspeed = new_speed * new_gamma;
+}
+
+void PhaseBubble::evolveWall2(numType dt, numType dE) {
+  /*
+  Bubble evolution using energy conservation. Energy change can be calculated
+  from collisions. R_(i+1) = R_i + V_i * dt V_(i+1) = V_(R, E)
+  */
+
+  numType new_radius;  // R_(i+1)
+  numType new_speed;   // V_(i+1)
+  numType new_gamma;
+  numType new_energy;
+
+  new_radius = m_bubble.radius + dt * m_bubble.speed;
+  new_energy = m_energy + dE;
+
+  new_speed = std::sqrt(
+      1 - std::pow(
+              (4 * M_PI * m_sigma * std::pow(new_radius, 2.) /
+               (new_energy + 4. * M_PI / 3. * m_dV * std::pow(new_radius, 3.))),
+              2.));
+  new_gamma = 1.0 / std::sqrt(1. - std::pow(new_speed, 2.));
+
+  // Approximation for speed change -> Not accurate -> Need another approach
+  // Needed when bubble wall speed value changes sign
+  /*numType initialRadius = m_bubble.radius;
+  numType dR = m_bubble.speed * dt;
+  numType speedChange1 = 144 * std::pow(M_PI, 2.) *
+                         std::pow(initialRadius, 3.) * std::pow(m_sigma, 2.);
+  numType speedChange2 =
+      (3 * initialRadius * dE +
+       4 * M_PI * std::pow(initialRadius, 3.) * m_dV * dR - 6 * m_energy * dR);
+  numType speedChange3 = std::pow(
+      4 * M_PI * std::pow(initialRadius, 3.) * m_dV + 3 * m_energy, 3.);
+  numType speedChange4 =
+      4 * M_PI * std::pow(initialRadius, 2.) * m_sigma /
+      (4 * M_PI / 3 * std::pow(initialRadius, 3.) * m_dV + m_energy);
+
+  numType speedChange5 = std::sqrt(1 - speedChange4);
+  numType speedChange =
+      speedChange1 * speedChange2 / (speedChange3 * speedChange5);
+  std::cout << speedChange << ", " << new_speed - m_bubble.speed << std::endl;*/
+
+  m_bubble.radius = new_radius;
+  m_bubble.speed = new_speed;
+  m_bubble.gamma = new_gamma;
+  m_energy = new_energy;
   m_bubble.radius2 = std::pow(m_bubble.radius, 2);
   m_bubble.gammaXspeed = m_bubble.speed * m_bubble.gamma;
 }