diff --git a/.gitignore b/.gitignore index ffa735c..66f7ac7 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,9 @@ dependencies/OpenCL bubbleSim/data +# External tool files +*.cltracer + # User-specific files *.rsuser diff --git a/bubbleSim.sln b/bubbleSim.sln index 3846372..45ef9df 100644 --- a/bubbleSim.sln +++ b/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 diff --git a/bubbleSim/bubbleSim.vcxproj b/bubbleSim/bubbleSim.vcxproj index e8f3e43..d4998d0 100644 --- a/bubbleSim/bubbleSim.vcxproj +++ b/bubbleSim/bubbleSim.vcxproj @@ -1,6 +1,14 @@  + + bubbleSim test + Win32 + + + bubbleSim test + x64 + Debug Win32 @@ -65,6 +73,12 @@ v143 Unicode + + Application + true + v143 + Unicode + Application false @@ -78,6 +92,12 @@ v143 Unicode + + Application + true + v143 + Unicode + Application false @@ -93,12 +113,18 @@ + + + + + + @@ -108,6 +134,11 @@ $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ false + + $(SolutionDir)\bin\$(Platform)\$(Configuration)\ + $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ + false + $(SolutionDir)\bin\$(Platform)\$(Configuration)\ $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ @@ -118,6 +149,11 @@ $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ false + + $(SolutionDir)\bin\$(Platform)\$(Configuration)\ + $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ + false + $(SolutionDir)\bin\$(Platform)\$(Configuration)\ $(SolutionDir)\bin\intermediates\$(Platform)\$(Configuration)\ @@ -138,6 +174,21 @@ $(SolutionDir)\dependencies\OpenCL\lib; + + + Level3 + true + WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + $(SolutionDir)\dependencies\OpenCL\include; $(SolutionDir)\dependencies\json + stdcpp17 + + + Console + true + $(SolutionDir)\dependencies\OpenCL\lib; + + Level3 @@ -172,6 +223,21 @@ $(SolutionDir)\dependencies\OpenCL\lib; + + + Level3 + true + _DEBUG;_CONSOLE;%(PreprocessorDefinitions) + true + $(SolutionDir)\dependencies\OpenCL\include; $(SolutionDir)\dependencies\json + stdcpp17 + + + Console + true + $(SolutionDir)\dependencies\OpenCL\lib; + + Level3 diff --git a/bubbleSim/bubbleSim.vcxproj.filters b/bubbleSim/bubbleSim.vcxproj.filters index a988896..6a98c73 100644 --- a/bubbleSim/bubbleSim.vcxproj.filters +++ b/bubbleSim/bubbleSim.vcxproj.filters @@ -16,81 +16,90 @@ {ad3d76ba-6f94-404f-bc5d-04c6fe891d8f} - - {ead60fb0-b85b-4979-a0a9-ef637c5cf283} - - - {0187af86-67e0-4ada-9316-2c50a6d1db91} - {70c9ae99-7764-4421-ad94-5ac76b6b096b} - - {a5c559dd-9d80-4aef-a16b-8988be73efcd} + + {065cebb1-beb8-4e96-9821-6f828b7743fa} + + + {8271bb38-eced-4c45-837f-66505dc396ac} + + + {25677fb9-0349-4027-b7fb-4c5eac0b861c} + + + {6c4c2689-9937-4d27-8cef-7ffb46d3e1fa} + + + {fc07c3d0-554a-4313-8c45-36b30388084b} - - {729ee2fd-1b8b-46ac-87ea-9dccda36f30d} + + {e2bc3991-b791-4914-8ca2-bb079843d4d7} + + + {d00f1d28-f39c-485c-bc09-b2ad361e7fee} - - Header Files - Header Files Header Files - + + Source Files\Generators + + Header Files + + Physics + - Header Files\Objects + Physics\Bubble - - Header Files\Objects + + Physics\Dynamics - - Header Files\Objects + + Physics\Particles - - Source Files\Generators + + Stream - Header Files - - - Header Files + Config - - Header Files\Streamers + + OpenCL - - Header Files + + Physics\Dynamics - - Source Files - Source Files - - Source Files - - Source Files\Objects - - - Source Files\Objects + Physics\Bubble - Source Files + Physics\Dynamics + + + Physics\Particles - Source Files\Streamers + Stream + + + OpenCL + + + Physics\Dynamics diff --git a/bubbleSim/collision.cpp b/bubbleSim/collision.cpp index fb64f8b..3b498f8 100644 --- a/bubbleSim/collision.cpp +++ b/bubbleSim/collision.cpp @@ -3,6 +3,12 @@ CollisionCellCollection::CollisionCellCollection( numType t_meanFreePath, unsigned int t_cellCountInOneAxis, bool t_doubleCellCount, cl::Context& cl_context) { + /* + TODO: + Differentiate if particle is inside or outside the bubble. (Different + collision cells) Change name of t_doubleCellCount to more accurate name. + + */ int openCLerrNum; if (t_cellCountInOneAxis % 2 != 1) { @@ -52,68 +58,87 @@ CollisionCellCollection::CollisionCellCollection( void CollisionCellCollection::generateShiftVector( RandomNumberGenerator& t_rng) { - m_shiftVector = {m_cellLength / 2. - t_rng.generate_number() * m_cellLength, - m_cellLength / 2. - t_rng.generate_number() * m_cellLength, - m_cellLength / 2. - t_rng.generate_number() * m_cellLength}; + m_shiftVector = { + m_cellLength / (numType)2. - t_rng.generate_number() * m_cellLength, + m_cellLength / (numType)2. - t_rng.generate_number() * m_cellLength, + m_cellLength / (numType)2. - t_rng.generate_number() * m_cellLength}; } void CollisionCellCollection::recalculate_cells( std::vector& t_particles, RandomNumberGenerator& t_rng) { - // 0 index cell is for particles outside the collision cell grid + numType phi, theta; - double phi, theta; + /* + * 0 index: Sum pX + * 1 index: Sum pY + * 2 index: Sum pZ + * 3 index: Sum E + * 4 index: Particle count in a cell + * 5 index: Prod(E_i) in a cell + */ + std::vector> frames( + m_cellCount, {(cl_numType)0., (cl_numType)0., (cl_numType)0., + (cl_numType)0., (cl_numType)0., (cl_numType)1.}); - std::vector> frames; - frames.resize(m_cellCount); /* * Sum up all momentums, energies and masses for each cell */ - for (Particle particle : t_particles) { + for (Particle& particle : t_particles) { if (particle.idxCollisionCell == 0) continue; - frames[particle.idxCollisionCell][0] += particle.p_x; - frames[particle.idxCollisionCell][1] += particle.p_y; - frames[particle.idxCollisionCell][2] += particle.p_z; + frames[particle.idxCollisionCell][0] += particle.pX; + frames[particle.idxCollisionCell][1] += particle.pY; + frames[particle.idxCollisionCell][2] += particle.pZ; frames[particle.idxCollisionCell][3] += particle.E; frames[particle.idxCollisionCell][4] += 1; + frames[particle.idxCollisionCell][5] *= particle.E; } /* * Calculate velocities for each collision cell */ for (size_t i = 1; i < m_cellCount; i++) { - m_collisionCells[i].particle_count = (int)frames[i][4]; - if (frames[i][4] > 1) { - m_collisionCells[i].v_x = frames[i][0] / frames[i][3]; - m_collisionCells[i].v_y = frames[i][1] / frames[i][3]; - m_collisionCells[i].v_z = frames[i][2] / frames[i][3]; - m_collisionCells[i].total_mass = - std::sqrt(std::pow(frames[i][3], 2) - std::pow(frames[i][0], 2) - - std::pow(frames[i][1], 2) - std::pow(frames[i][2], 2)); - m_collisionCells[i].p_x = frames[i][0] / frames[i][4]; - m_collisionCells[i].p_y = frames[i][1] / frames[i][4]; - m_collisionCells[i].p_z = frames[i][2] / frames[i][4]; - m_collisionCells[i].p_E = frames[i][3] / frames[i][4]; - - /* - * beta = sqrt(v_x^2 + v_y^2 + v_z^2); - */ - m_collisionCells[i].v2 = - std::fma(m_collisionCells[i].v_x, m_collisionCells[i].v_x, - std::fma(m_collisionCells[i].v_y, m_collisionCells[i].v_y, - m_collisionCells[i].v_z * m_collisionCells[i].v_z)); - /* - * gamma = 1/sqrt(1-beta^2) - */ - m_collisionCells[i].gamma = 1 / std::sqrt(1 - m_collisionCells[i].v2); - phi = std::acos(1 - 2 * t_rng.generate_number()); - theta = 2 * M_PI * t_rng.generate_number(); - - /* - * Generate rotation axis and angle for CollisionCell - */ - m_collisionCells[i].theta = 2 * M_PI * t_rng.generate_number(); - m_collisionCells[i].x = std::sin(phi) * std::cos(theta); - m_collisionCells[i].y = std::sin(phi) * std::sin(theta); - m_collisionCells[i].z = std::cos(phi); + if ((cl_uint)frames[i][4] < 2) { + m_collisionCells[i].particle_count = (cl_uint)0; + continue; } + /* + * If T ^ N / prod_N(E_i) <= RNG then skip rotation + * T - temperature, prod(E_i) is product of particles' energies in a cell, + * N is number of particles in a cell + */ + + if (t_rng.generate_number() >= + std::pow(0.01, frames[i][4]) / frames[i][5]) { + m_collisionCells[i].particle_count = (cl_uint)0; + continue; + } + m_collisionCells[i].particle_count = (int)frames[i][4]; + m_collisionCells[i].vX = frames[i][0] / frames[i][3]; + m_collisionCells[i].vY = frames[i][1] / frames[i][3]; + m_collisionCells[i].vZ = frames[i][2] / frames[i][3]; + m_collisionCells[i].total_mass = + std::sqrt(std::pow(frames[i][3], 2) - std::pow(frames[i][0], 2) - + std::pow(frames[i][1], 2) - std::pow(frames[i][2], 2)); + m_collisionCells[i].pX = frames[i][0] / frames[i][4]; + m_collisionCells[i].pY = frames[i][1] / frames[i][4]; + m_collisionCells[i].pZ = frames[i][2] / frames[i][4]; + m_collisionCells[i].pE = frames[i][3] / frames[i][4]; + + m_collisionCells[i].v2 = + std::fma(m_collisionCells[i].vX, m_collisionCells[i].vX, + std::fma(m_collisionCells[i].vY, m_collisionCells[i].vY, + m_collisionCells[i].vZ * m_collisionCells[i].vZ)); + + m_collisionCells[i].gamma = 1 / std::sqrt(1 - m_collisionCells[i].v2); + phi = std::acos((numType)1. - (numType)2. * t_rng.generate_number()); + theta = (numType)2 * (numType)M_PI * t_rng.generate_number(); + + /* + * Generate rotation axis and angle for CollisionCell + */ + m_collisionCells[i].theta = + (numType)2. * (numType)M_PI * t_rng.generate_number(); + m_collisionCells[i].x = std::sin(phi) * std::cos(theta); + m_collisionCells[i].y = std::sin(phi) * std::sin(theta); + m_collisionCells[i].z = std::cos(phi); } } \ No newline at end of file diff --git a/bubbleSim/collision.h b/bubbleSim/collision.h index 0488c3d..f8da4c3 100644 --- a/bubbleSim/collision.h +++ b/bubbleSim/collision.h @@ -2,23 +2,25 @@ #include "base.h" #include "objects.h" typedef struct CollisionCell { - cl_numType v_x; - cl_numType v_y; - cl_numType v_z; + // Lorentz transformation + cl_numType gamma; + cl_numType vX; + cl_numType vY; + cl_numType vZ; + // Rotation matrix cl_numType x; cl_numType y; cl_numType z; - cl_numType theta; - cl_numType p_E; - cl_numType p_x; - cl_numType p_y; - cl_numType p_z; + // Average 4-momentum values + cl_numType pE; + cl_numType pX; + cl_numType pY; + cl_numType pZ; cl_numType v2; - cl_numType gamma; cl_numType total_mass; cl_uint particle_count; } CollisionCell; @@ -30,6 +32,8 @@ class CollisionCellCollection { std::array& getShiftVector() { return m_shiftVector; } + u_int getCellCount() { return m_cellCount; } + void recalculate_cells(std::vector& t_particles, RandomNumberGenerator& t_rng); diff --git a/bubbleSim/config_reader.hpp b/bubbleSim/config_reader.hpp index fe71473..97d3447 100644 --- a/bubbleSim/config_reader.hpp +++ b/bubbleSim/config_reader.hpp @@ -7,7 +7,7 @@ class ConfigReader { public: std::string kernelName; int m_seed; - int m_maxSteps; + u_int m_maxSteps; numType dt; numType maxTime; bool cyclicBoundaryOn; diff --git a/bubbleSim/datastreamer.cpp b/bubbleSim/datastreamer.cpp index 9fd885d..003acae 100644 --- a/bubbleSim/datastreamer.cpp +++ b/bubbleSim/datastreamer.cpp @@ -26,6 +26,33 @@ void DataStreamer::initStream_Data() { << std::endl; } +void DataStreamer::initStream_Momentum(size_t t_binsCount, + numType t_maxMomentumValue) { + m_initialized_Momentum = true; + m_binsCount_Momentum = t_binsCount; + m_maxMomentum_Momentum = t_maxMomentumValue; + m_dp_Momentum = t_maxMomentumValue / t_binsCount; + + m_stream_MomentumX.open(m_filePath / "pX.csv", std::ios::out); + m_stream_MomentumY.open(m_filePath / "pY.csv", std::ios::out); + m_stream_MomentumZ.open(m_filePath / "pZ.csv", std::ios::out); + + m_stream_MomentumX << t_binsCount << "," << t_maxMomentumValue << "\n"; + m_stream_MomentumY << t_binsCount << "," << t_maxMomentumValue << "\n"; + m_stream_MomentumZ << t_binsCount << "," << t_maxMomentumValue << "\n"; + for (size_t i = 1; i <= t_binsCount; i++) { + if (i == t_binsCount) { + m_stream_MomentumX << i << std::endl; + m_stream_MomentumY << i << std::endl; + m_stream_MomentumZ << i << std::endl; + } else { + m_stream_MomentumX << i << ","; + m_stream_MomentumY << i << ","; + m_stream_MomentumZ << i << ","; + } + } +} + void DataStreamer::initStream_MomentumIn(size_t t_binsCount, numType t_maxMomentumValue) { m_initialized_MomentumIn = true; @@ -170,6 +197,9 @@ void DataStreamer::stream(Simulation& simulation, numType totalEnergy; // Save momentum data + std::vector bins_MomentumX; + std::vector bins_MomentumY; + std::vector bins_MomentumZ; std::vector bins_MomentumIn; std::vector bins_MomentumOut; @@ -184,6 +214,11 @@ void DataStreamer::stream(Simulation& simulation, std::vector bins_TangentialVelocityCount; // Initialize variables + if (m_initialized_Momentum) { + bins_MomentumX.resize(m_binsCount_Momentum, 0); + bins_MomentumY.resize(m_binsCount_Momentum, 0); + bins_MomentumZ.resize(m_binsCount_Momentum, 0); + } if (m_initialized_MomentumIn) { bins_MomentumIn.resize(m_binsCount_MomentumIn, 0); } @@ -222,6 +257,17 @@ void DataStreamer::stream(Simulation& simulation, particleRadius = particleCollection.calculateParticleRadius(i); particleMomentum = particleCollection.calculateParticleMomentum(i); + if (m_initialized_Momentum) { + bins_MomentumX[std::clamp( + (int)(abs(particleCollection.getParticles()[i].pX / m_dp_Momentum)), + 0, (int)m_binsCount_Momentum - 1)] += 1; + bins_MomentumY[std::clamp( + (int)(abs(particleCollection.getParticles()[i].pY / m_dp_Momentum)), + 0, (int)m_binsCount_Momentum - 1)] += 1; + bins_MomentumZ[std::clamp( + (int)(abs(particleCollection.getParticles()[i].pZ / m_dp_Momentum)), + 0, (int)m_binsCount_Momentum - 1)] += 1; + } if (m_initialized_Density && (particleRadius < m_maxRadius_Density)) { bins_Density[(int)(particleRadius / m_dr_Density)] += 1; } @@ -234,7 +280,13 @@ void DataStreamer::stream(Simulation& simulation, (particleRadius < m_maxRadius_RadialVelocity)) { // Average A(N) = [A(N-1) * (N-1) + a(N) ]/N particleRadialVelocity = - particleCollection.calculateParticleRadialVelocity(i); + (particleCollection.getParticles()[i].x * + particleCollection.getParticles()[i].pX + + particleCollection.getParticles()[i].y * + particleCollection.getParticles()[i].pY + + particleCollection.getParticles()[i].z * + particleCollection.getParticles()[i].pZ) / + (particleCollection.getParticles()[i].E * particleRadius); bins_RadialVelocity[(int)(particleRadius / m_dr_RadialVelocity)] = (bins_RadialVelocity[(int)(particleRadius / m_dr_RadialVelocity)] * bins_RadialVelocityCount[(int)(particleRadius / @@ -249,8 +301,12 @@ void DataStreamer::stream(Simulation& simulation, } if (m_initialized_TangentialVelocity && (particleRadius < m_maxRadius_TangentialVelocity)) { - particleTangentialVelocity = - particleCollection.calculateParticleTangentialVelocity(i); + if (m_initialized_RadialVelocity) { + particleTangentialVelocity = std::sqrt(1 - particleRadialVelocity); + } else { + particleTangentialVelocity = + particleCollection.calculateParticleTangentialVelocity(i); + } bins_TangentialVelocity[(int)(particleRadius / m_dr_TangentialVelocity)] = (bins_TangentialVelocity[(int)(particleRadius / m_dr_TangentialVelocity)] * @@ -352,6 +408,16 @@ void DataStreamer::stream(Simulation& simulation, m_stream_TangentialVelocity << bins_TangentialVelocity[m_binsCount_TangentialVelocity - 1] << "\n"; } + if (m_initialized_Momentum) { + for (size_t i = 0; i < m_binsCount_Momentum - 1; i++) { + m_stream_MomentumX << bins_MomentumX[i] << ","; + m_stream_MomentumY << bins_MomentumY[i] << ","; + m_stream_MomentumZ << bins_MomentumZ[i] << ","; + } + m_stream_MomentumX << bins_MomentumX[m_binsCount_Momentum - 1] << "\n"; + m_stream_MomentumY << bins_MomentumY[m_binsCount_Momentum - 1] << "\n"; + m_stream_MomentumZ << bins_MomentumZ[m_binsCount_Momentum - 1] << "\n"; + } /*auto programEndTime = std::chrono::high_resolution_clock::now(); std::cout << "Time taken (stream): " << std::chrono::duration_cast( @@ -568,7 +634,7 @@ void DataStreamer::streamTangentialVelocity( t_stream << average_velocity[t_binsCount - 1] << "\n"; } -void DataStreamer::StreamRadialMomentumProfile( +void DataStreamer::streamRadialMomentumProfile( std::ofstream& t_stream, size_t t_binsCountRadius, size_t t_binsCountMomentum, numType t_minRadiusValue, numType t_maxRadiusValue, numType t_minMomentumValue, diff --git a/bubbleSim/datastreamer.h b/bubbleSim/datastreamer.h index 3faff42..825f281 100644 --- a/bubbleSim/datastreamer.h +++ b/bubbleSim/datastreamer.h @@ -39,6 +39,7 @@ class DataStreamer { void initStream_TangentialVelocity(size_t t_binsCount, numType t_maxRadiusValue); void initStream_radialMomentum(size_t t_binsCount, numType t_maxRadiusValue); + void initStream_Momentum(size_t t_binsCount, numType t_maxMomentumValue); void stream(Simulation& simulation, ParticleCollection& particleCollection, PhaseBubble& bubble, cl::CommandQueue& cl_queue); @@ -74,7 +75,7 @@ class DataStreamer { numType t_maxRadiusValue, ParticleCollection& particleCollection, cl::CommandQueue& cl_queue); - void StreamRadialMomentumProfile( + void streamRadialMomentumProfile( std::ofstream& t_stream, size_t t_binsCountRadius, size_t t_binsCountMomentum, numType t_minRadiusValue, numType t_maxRadiusValue, numType t_minMomentumValue, @@ -85,6 +86,9 @@ class DataStreamer { std::filesystem::path m_filePath; std::ofstream m_stream_Data; + std::ofstream m_stream_MomentumX; + std::ofstream m_stream_MomentumY; + std::ofstream m_stream_MomentumZ; std::ofstream m_stream_MomentumIn; std::ofstream m_stream_MomentumOut; @@ -98,6 +102,7 @@ class DataStreamer { // General data about simulation state bool m_initialized_Data = false; // Momentum profiles in and outside the bubble + bool m_initialized_Momentum = false; bool m_initialized_MomentumIn = false; bool m_initialized_MomentumOut = false; // Energy and number desnity profiles @@ -108,6 +113,7 @@ class DataStreamer { bool m_initialized_TangentialVelocity = false; bool m_initialized_RadialMomentum = false; + size_t m_binsCount_Momentum; size_t m_binsCount_MomentumIn; size_t m_binsCount_MomentumOut; size_t m_binsCount_Density; @@ -115,6 +121,7 @@ class DataStreamer { size_t m_binsCount_RadialVelocity; size_t m_binsCount_TangentialVelocity; + numType m_maxMomentum_Momentum; numType m_maxMomentum_MomentumIn; numType m_maxMomentum_MomentumOut; numType m_maxRadius_Density; @@ -122,6 +129,7 @@ class DataStreamer { numType m_maxRadius_RadialVelocity; numType m_maxRadius_TangentialVelocity; + numType m_dp_Momentum; numType m_dp_MomentumIn; numType m_dp_MomentumOut; numType m_dr_Density; diff --git a/bubbleSim/opencl_kernels.cpp b/bubbleSim/opencl_kernels.cpp index 8e998b6..04c3aa3 100644 --- a/bubbleSim/opencl_kernels.cpp +++ b/bubbleSim/opencl_kernels.cpp @@ -1,23 +1,31 @@ #include "opencl_kernels.h" OpenCLLoader::OpenCLLoader(std::string kernelPath) { - std::string particleBubbleStepKernelName = "particle_bubble_step"; - std::string cellAssignKernelName = "assign_cell_index_to_particle"; - std::string transformKernelName = "transform_momentum"; - std::string particleStepKernelName = "particle_step"; - std::string particleBounceKernelName = "particle_bounce"; - std::string particleBubbleBoundaryStepKernelName = - "particle_bubble_step_cyclic"; + createContext(m_devices); createProgram(m_context, m_deviceUsed, kernelPath); + + + std::string particleBubbleStepKernelName = "particle_bubble_step"; createKernel(m_program, m_particleBubbleStepKernel, particleBubbleStepKernelName.c_str()); // cl::Kernel + + std::string transformKernelName = "rotate_momentum"; createKernel(m_program, m_rotationKernel, transformKernelName.c_str()); + + std::string cellAssignKernelName = "assign_particle_to_collision_cell"; createKernel(m_program, m_cellAssignmentKernel, cellAssignKernelName.c_str()); + + std::string particleStepKernelName = "particle_step"; createKernel(m_program, m_particleStepKernel, particleStepKernelName.c_str()); + + std::string particleBounceKernelName = "particle_boundary_check"; createKernel(m_program, m_particleBounceKernel, particleBounceKernelName.c_str()); + + std::string particleBubbleBoundaryStepKernelName = + "particle_bubble_step_cyclic"; createKernel(m_program, m_particleBubbleBoundaryStepKernel, particleBubbleBoundaryStepKernelName.c_str()); @@ -25,17 +33,31 @@ OpenCLLoader::OpenCLLoader(std::string kernelPath) { } OpenCLLoader::OpenCLLoader(std::string kernelPath, std::string kernelName) { - std::string particleBubbleStepKernelName = "particle_bubble_step"; - std::string particleBubbleBoundaryStepKernelName = - "particle_bubble_step_cyclic"; - createContext(m_devices); createProgram(m_context, m_deviceUsed, kernelPath); + + // For collision + std::string cellAssignKernelName = "assign_particle_to_collision_cell"; + createKernel(m_program, m_cellAssignmentKernel, cellAssignKernelName.c_str()); + + std::string transformKernelName = "rotate_momentum"; + createKernel(m_program, m_rotationKernel, transformKernelName.c_str()); + + std::string particleBounceKernelName = "particle_boundary_check"; + createKernel(m_program, m_particleBounceKernel, + particleBounceKernelName.c_str()); + createKernel(m_program, m_kernel, kernelName.c_str()); // cl::Kernel + + std::string particleBubbleStepKernelName = "particle_bubble_step"; createKernel(m_program, m_particleBubbleStepKernel, particleBubbleStepKernelName.c_str()); // cl::Kernel + + std::string particleBubbleBoundaryStepKernelName = + "particle_bubble_step_cyclic"; createKernel(m_program, m_particleBubbleBoundaryStepKernel, particleBubbleBoundaryStepKernelName.c_str()); + createQueue(m_context, m_deviceUsed); } diff --git a/bubbleSim/particle.cpp b/bubbleSim/particle.cpp index 0c311e8..2db072f 100644 --- a/bubbleSim/particle.cpp +++ b/bubbleSim/particle.cpp @@ -21,7 +21,7 @@ void ParticleGenerator::calculateCPD(numType t_temperature, numType t_pMax, m_cumulativeProbabilityFunction[0].reserve(vectorSize); m_cumulativeProbabilityFunction[1].reserve(vectorSize); - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType lastCPFValue = 0.; numType lastMomentumValue = 0.; @@ -32,7 +32,7 @@ void ParticleGenerator::calculateCPD(numType t_temperature, numType t_pMax, // Integral of: dp * p^2 * exp(-sqrt(p^2+m^2)/T) m_cumulativeProbabilityFunction[0].push_back( lastCPFValue + - t_dp * std::pow(lastMomentumValue, 2) * + t_dp * std::pow(lastMomentumValue, (numType)2.) * std::exp( -std::sqrt(std::fma(lastMomentumValue, lastMomentumValue, m2)) / t_temperature)); @@ -87,9 +87,9 @@ numType ParticleGenerator::interp(numType t_xValue, void ParticleGenerator::generateRandomDirection( numType& x, numType& y, numType& z, numType t_radius, RandomNumberGenerator& t_generator) { - numType phi = - std::acos(1 - 2 * t_generator.generate_number()); // inclination - numType theta = 2 * M_PI * t_generator.generate_number(); + numType phi = std::acos( + (numType)1. - (numType)2. * t_generator.generate_number()); // inclination + numType theta = (numType)2. * (numType)M_PI * t_generator.generate_number(); x = t_radius * std::sin(phi) * std::cos(theta); // x y = t_radius * std::sin(phi) * std::sin(theta); // y z = t_radius * std::cos(phi); // z @@ -143,7 +143,7 @@ numType ParticleGenerator::generateNParticlesInBox( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; for (u_int i = 0; i < t_N; i++) { @@ -152,7 +152,7 @@ numType ParticleGenerator::generateNParticlesInBox( // Generates 3D space coordinates and pushes to m_P vector generateParticleMomentum(p_x, p_y, p_z, pValue, t_generator); - E = std::sqrt(m2 + pow(pValue, 2)); + E = std::sqrt(m2 + pow(pValue, (numType)2.)); totalEnergy += E; t_particles.push_back(Particle{x, y, z, p_x, p_y, p_z, E, m_mass}); } @@ -167,7 +167,7 @@ numType ParticleGenerator::generateNParticlesInBox( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; for (u_int i = 0; i < t_N; i++) { @@ -192,7 +192,7 @@ numType ParticleGenerator::generateNParticlesInBox( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; numType radius; for (u_int i = 0; i < t_N; i++) { @@ -203,7 +203,7 @@ numType ParticleGenerator::generateNParticlesInBox( } while (radius < t_radiusIn); // Generates 3D space coordinates and pushes to m_P vector generateParticleMomentum(p_x, p_y, p_z, pValue, t_generator); - E = std::sqrt(m2 + pow(pValue, 2)); + E = std::sqrt(m2 + pow(pValue, (numType)2.)); totalEnergy += E; t_particles.push_back(Particle{x, y, z, p_x, p_y, p_z, E, m_mass}); } @@ -219,7 +219,7 @@ numType ParticleGenerator::generateNParticlesInBox( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; numType radius; for (u_int i = 0; i < t_N; i++) { @@ -231,7 +231,7 @@ numType ParticleGenerator::generateNParticlesInBox( } while (radius < t_radiusIn); // Generates 3D space coordinates and pushes to m_P vector generateParticleMomentum(p_x, p_y, p_z, pValue, t_generator); - E = std::sqrt(m2 + pow(pValue, 2)); + E = std::sqrt(m2 + pow(pValue, (numType)2.)); totalEnergy += E; t_particles.push_back(Particle{x, y, z, p_x, p_y, p_z, E, m_mass}); } @@ -245,7 +245,7 @@ numType ParticleGenerator::generateNParticlesInSphere( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; numType radius; @@ -258,7 +258,7 @@ numType ParticleGenerator::generateNParticlesInSphere( // Generates 3D space coordinates and pushes to m_P vector generateParticleMomentum(p_x, p_y, p_z, pValue, t_generator); - E = std::sqrt(m2 + pow(pValue, 2)); + E = std::sqrt(m2 + pow(pValue, (numType)2.)); totalEnergy += E; t_particles.push_back(Particle{x, y, z, p_x, p_y, p_z, E, m_mass}); } @@ -272,7 +272,7 @@ numType ParticleGenerator::generateNParticlesInSphere( numType x, y, z; numType p_x, p_y, p_z; numType E; - numType m2 = std::pow(m_mass, 2); + numType m2 = std::pow(m_mass, (numType)2.); numType pValue; numType radius; @@ -285,7 +285,7 @@ numType ParticleGenerator::generateNParticlesInSphere( // Generates 3D space coordinates and pushes to m_P vector generateParticleMomentum(p_x, p_y, p_z, pValue, t_generator); - E = std::sqrt(m2 + pow(pValue, 2)); + E = std::sqrt(m2 + pow(pValue, (numType)2.)); totalEnergy += E; t_particles.push_back(Particle{x, y, z, p_x, p_y, p_z, E, m_mass}); } @@ -315,7 +315,8 @@ ParticleCollection::ParticleCollection( m_particleCountIn = t_particleCountFalse; m_particleCountOut = t_particleCountTrue; } - m_massDelta2 = std::abs(std::pow(t_massTrue, 2) - std::pow(t_massFalse, 2)); + m_massDelta2 = std::abs(std::pow(t_massTrue, (numType)2.) - + std::pow(t_massFalse, (numType)2.)); m_massInBuffer = cl::Buffer(cl_context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, @@ -352,7 +353,7 @@ ParticleCollection::ParticleCollection( m_particleCountTotal * sizeof(Particle), m_particles.data(), &openCLerrNum); - m_dP = std::vector(m_particleCountTotal, 0.); + m_dP = std::vector(m_particleCountTotal, (numType)0.); m_dPBuffer = cl::Buffer(cl_context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, m_particleCountTotal * sizeof(numType), m_dP.data(), &openCLerrNum); @@ -393,28 +394,29 @@ numType ParticleCollection::calculateParticleMomentum(u_int i) { // return std::sqrt(m_P[3 * i] * m_P[3 * i] + m_P[3 * i + 1] * m_P[3 * i + 1] // + m_P[3 * i + 2] * m_P[3 * i + 2]); - return std::sqrt(std::fma(m_particles[i].p_x, m_particles[i].p_x, - std::fma(m_particles[i].p_y, m_particles[i].p_y, - m_particles[i].p_z * m_particles[i].p_z))); + return std::sqrt(std::fma(m_particles[i].pX, m_particles[i].pX, + std::fma(m_particles[i].pY, m_particles[i].pY, + m_particles[i].pZ * m_particles[i].pZ))); } numType ParticleCollection::calculateParticleEnergy(u_int i) { return std::sqrt( - std::fma(m_particles[i].p_x, m_particles[i].p_x, - std::fma(m_particles[i].p_y, m_particles[i].p_y, - std::fma(m_particles[i].p_z, m_particles[i].p_z, + std::fma(m_particles[i].pX, m_particles[i].pX, + std::fma(m_particles[i].pY, m_particles[i].pY, + std::fma(m_particles[i].pZ, m_particles[i].pZ, m_particles[i].m * m_particles[i].m)))); } numType ParticleCollection::calculateParticleRadialVelocity(u_int i) { - return std::fma(m_particles[i].p_x, m_particles[i].x, - std::fma(m_particles[i].p_y, m_particles[i].y, - m_particles[i].p_z * m_particles[i].z)) / + return std::fma(m_particles[i].pX, m_particles[i].x, + std::fma(m_particles[i].pY, m_particles[i].y, + m_particles[i].pZ * m_particles[i].z)) / (m_particles[i].E * calculateParticleRadius(i)); } numType ParticleCollection::calculateParticleTangentialVelocity(u_int i) { - return std::sqrt(1 - std::pow(calculateParticleRadialVelocity(i), 2)); + return std::sqrt((numType)1. - + std::pow(calculateParticleRadialVelocity(i), (numType)2.)); } numType ParticleCollection::calculateNumberDensity(numType t_mass, @@ -423,13 +425,13 @@ numType ParticleCollection::calculateNumberDensity(numType t_mass, numType t_pMax) { numType n = 0; numType p = 0; - numType m2 = std::pow(t_mass, 2); + numType m2 = std::pow(t_mass, (numType)2.); for (; p <= t_pMax; p += t_dp) { - n += t_dp * std::pow(p, 2) * + n += t_dp * std::pow(p, (numType)2.) * std::exp(-std::sqrt(std::fma(p, p, m2)) / t_temperature); } - n = n / (numType)(2. * std::pow(M_PI, 2)); + n = n / (numType)((numType)2. * std::pow((numType)M_PI, (numType)2.)); return n; } @@ -439,16 +441,16 @@ numType ParticleCollection::calculateEnergyDensity(numType t_mass, numType t_pMax) { numType rho = 0; numType p = 0; - numType m2 = std::pow(t_mass, 2); + numType m2 = std::pow(t_mass, (numType)2.); numType sqrt_p2_m2; for (; p <= t_pMax; p += t_dp) { sqrt_p2_m2 = std::sqrt(std::fma(p, p, m2)); - rho += t_dp * std::pow(p, 2) * sqrt_p2_m2 * + rho += t_dp * std::pow(p, (numType)2.) * sqrt_p2_m2 * std::exp(-sqrt_p2_m2 / t_temperature); } - rho = rho / (numType)(2 * std::pow(M_PI, 2)); + rho = rho / (numType)((numType)2. * std::pow(M_PI, (numType)2.)); return rho; } @@ -456,7 +458,8 @@ numType ParticleCollection::calculateEnergyDensity(numType t_mass, numType ParticleCollection::countParticleNumberDensity(numType t_radius1) { u_int counter = 0; - numType volume = 4 * (pow(t_radius1, 3) * M_PI) / 3; + numType volume = + (numType)4. * (pow(t_radius1, (numType)3.) * (numType)M_PI) / (numType)3.; for (int i = 0; i < m_particleCountTotal; i++) { if (calculateParticleRadius(i) < t_radius1) { counter += 1; @@ -468,7 +471,10 @@ numType ParticleCollection::countParticleNumberDensity(numType t_radius1) { numType ParticleCollection::countParticleNumberDensity(numType t_radius1, numType t_radius2) { u_int counter = 0; - numType volume = 4 * ((pow(t_radius2, 3) - pow(t_radius1, 3)) * M_PI) / 3; + numType volume = (numType)4. * + ((pow(t_radius2, (numType)3.) - pow(t_radius1, (numType)3.)) * + (numType)M_PI) / + (numType)3.; for (int i = 0; i < m_particleCountTotal; i++) { if ((calculateParticleRadius(i) > t_radius1) && (calculateParticleRadius(i) < t_radius2)) { @@ -481,7 +487,8 @@ numType ParticleCollection::countParticleNumberDensity(numType t_radius1, numType ParticleCollection::countParticleEnergyDensity(numType t_radius1) { numType energy = countParticlesEnergy(t_radius1); std::cout << energy << std::endl; - numType volume = 4 * (pow(t_radius1, 3) * M_PI) / 3; + numType volume = + (numType)4. * (pow(t_radius1, (numType)3.) * (numType)M_PI) / (numType)3.; return energy / volume; } @@ -528,8 +535,8 @@ void ParticleCollection::print_info(ConfigReader t_config, // exp(mu/T) -> ~ chemical potential -> std::cout << std::setprecision(6); numType exp_mu_T = - std::log((m_particleCountIn * std::pow(M_PI, 2)) / - (t_bubble.calculateVolume() * std::pow(m_temperatureFalse, 3))); + std::log((m_particleCountIn * std::pow((numType)M_PI, (numType)2.)) / + (t_bubble.calculateVolume() * std::pow(m_temperatureFalse, (numType)3))); numType nFromParameters = m_particleCountIn / t_bubble.calculateVolume(); // Assuming m = 0 numType rhoFromParameters = diff --git a/bubbleSim/particle.h b/bubbleSim/particle.h index b0abb40..bd26fe4 100644 --- a/bubbleSim/particle.h +++ b/bubbleSim/particle.h @@ -9,9 +9,9 @@ typedef struct Particle { cl_numType y; cl_numType z; - cl_numType p_x; - cl_numType p_y; - cl_numType p_z; + cl_numType pX; + cl_numType pY; + cl_numType pZ; cl_numType E; cl_numType m; @@ -343,9 +343,9 @@ class ParticleCollection { numType getParticleMomentum(u_int i) { return std::sqrt( - std::fma(m_particles[i].p_x, m_particles[i].p_x, - std::fma(m_particles[i].p_y, m_particles[i].p_y, - m_particles[i].p_z * m_particles[i].p_z))); + std::fma(m_particles[i].pX, m_particles[i].pX, + std::fma(m_particles[i].pY, m_particles[i].pY, + m_particles[i].pZ * m_particles[i].pZ))); } numType getParticleMass(u_int i) { return m_particles[i].m; } @@ -392,8 +392,8 @@ class ParticleCollection { void print_particle_info(unsigned int i) { Particle particle = m_particles[i]; std::cout << "(" << particle.x << ", " << particle.y << ", " << particle.z - << ") (" << particle.E << ", " << particle.p_x << ", " - << particle.p_y << ", " << particle.p_z << ")" << std::endl; + << ") (" << particle.E << ", " << particle.pX << ", " + << particle.pY << ", " << particle.pZ << ")" << std::endl; } void revertParticlesToLastStep(cl::CommandQueue& cl_queue) { @@ -452,4 +452,10 @@ class ParticleCollection { makeInteractedBubbleTrueStateCopy(); makePassedBubbleFalseStateCopy(); } + + void printParticleInfo(size_t i) { + Particle p = m_particles[i]; + std::printf("X: (%.7f, %.7f, %.7f), P: (%.7f, %.7f, %.7f, %.7f), m: %.7f\n", + p.x, p.y, p.z, p.E, p.pX, p.pY, p.pZ, p.m); + } }; diff --git a/bubbleSim/simulation.cpp b/bubbleSim/simulation.cpp index d7af051..a762e71 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -205,7 +205,7 @@ void Simulation::step(PhaseBubble& bubble, numType t_dP) { bubble.evolveWall(m_dt, t_dP); } -// In development: collisions +// Collision step void Simulation::step(ParticleCollection& particles, CollisionCellCollection& cells, RandomNumberGenerator& generator_collision, int i, @@ -214,12 +214,19 @@ void Simulation::step(ParticleCollection& particles, cl::Kernel& t_rotationKernel, cl::Kernel& t_particleBounceKernel, cl::CommandQueue& cl_queue) { - m_time += m_dt; + /* + * 1) Move particles + * 2) Solve boundaries + */ // Move particles cl_queue.enqueueNDRangeKernel(t_particleStepKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); - // Generate shift vector + cl_queue.enqueueNDRangeKernel( + t_particleBounceKernel, cl::NullRange, + cl::NDRange(particles.getParticleCountTotal())); + if (i % 1 == 0) { + // Generate shift vector cells.generateShiftVector(generator_collision); cells.writeShiftVectorBuffer(cl_queue); @@ -227,21 +234,24 @@ void Simulation::step(ParticleCollection& particles, cl_queue.enqueueNDRangeKernel( t_cellAssignmentKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); + // Update particle data on CPU particles.readParticlesBuffer(cl_queue); // Calculate COM and genrate rotation matrix for each cell - cells.recalculate_cells(particles.getParticles(), generator_collision); + // Update data on GPU particles.writeParticlesBuffer(cl_queue); cells.writeCollisionCellBuffer(cl_queue); - // Update momentum + // Rotate momentum cl_queue.enqueueNDRangeKernel( t_rotationKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); + //particles.readParticlesBuffer(cl_queue); + // Update simulation time + m_time += m_dt; + m_step += 1; } - cl_queue.enqueueNDRangeKernel(t_particleBounceKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); } diff --git a/bubbleSim/simulation.h b/bubbleSim/simulation.h index 6b7548c..db9d115 100644 --- a/bubbleSim/simulation.h +++ b/bubbleSim/simulation.h @@ -51,7 +51,7 @@ class Simulation { } }; - int getStep() { return m_step; } + unsigned int getStep() { return m_step; } void step(PhaseBubble& bubble, numType t_dP); /* @@ -100,11 +100,12 @@ class Simulation { } private: - // Sim time paramters: - // Cumulative time int m_seed; + + // Simulation time state numType m_time = 0.; - size_t m_step = 0; + u_int m_step = 0; + // One step time length numType m_dt; numType m_step_dt; diff --git a/bubbleSim/source.cpp b/bubbleSim/source.cpp index 973aee4..1e01cd1 100644 --- a/bubbleSim/source.cpp +++ b/bubbleSim/source.cpp @@ -76,7 +76,7 @@ int main(int argc, char* argv[]) { } // Collision is in development - bool b_collisionDevelopment = false; + bool b_collisionDevelopment = true; // Read configs and kernel file std::string s_configPath = argv[1]; // "config.json" @@ -106,29 +106,25 @@ int main(int argc, char* argv[]) { // 2) Initialize openCL (kernels, commandQueues) OpenCLLoader kernels(s_kernelPath, config.kernelName); - // 3) Define required physical parameters - numType temperatureFalse = - std::sqrt(std::abs(std::pow(config.particleMassTrue, 2) - - std::pow(config.particleMassFalse, 2))) / - config.parameterEta; - - numType temperatureTrue = 0; // -> No particles generated in true vacuum, - // thus set temperature in true vacuum 0 - // 4) Generate particles ParticleGenerator particleGenerator1; // 4.1) Create generator (calculates distribution) + if (b_collisionDevelopment) { - particleGenerator1 = ParticleGenerator(config.particleMassFalse, 3.); + particleGenerator1 = ParticleGenerator(config.particleMassFalse, 3.*config.particleTemperatureFalse); } else { - particleGenerator1 = - ParticleGenerator(config.particleMassFalse, temperatureFalse, - 30 * temperatureFalse, 1e-5 * temperatureFalse); + particleGenerator1 = ParticleGenerator( + config.particleMassFalse, config.particleTemperatureFalse, + 30 * config.particleTemperatureFalse, + 1e-5 * config.particleTemperatureFalse); } + // 4.2) Create arrays for particles which hold the data of the particles ParticleCollection particles( - config.particleMassTrue, config.particleMassFalse, temperatureTrue, - temperatureFalse, config.particleCountTrue, config.particleCountFalse, + config.particleMassTrue, config.particleMassFalse, + config.particleTemperatureTrue, + config.particleTemperatureFalse, config.particleCountTrue, + config.particleCountFalse, config.parameterCoupling, config.bubbleIsTrueVacuum, kernels.getContext()); @@ -162,24 +158,10 @@ int main(int argc, char* argv[]) { config.cyclicBoundaryRadius, kernels.getContext()); } - // 6) Define which openCL kernel to use. Kernel defines calculation process on - // the GPU. - // cl::Kernel* stepKernel; // In development: set kernel name in config file. // NB! Different kernels might need different input - /*if (!config.bubbleInteractionsOn) { - stepKernel = &kernels.m_particleStepKernel; - } else if ((config.bubbleInteractionsOn) && (!config.cyclicBoundaryOn)) { - stepKernel = &kernels.m_particleBubbleStepKernel; - } else if ((config.bubbleInteractionsOn) && (config.cyclicBoundaryOn)) { - stepKernel = &kernels.m_particleBubbleBoundaryStepKernel; - } else { - std::cerr << "Kernel for current configuration is not available" - << std::endl; - std::terminate(); - }*/ stepKernel = &kernels.m_kernel; // 7) Initialize bubble object @@ -214,14 +196,10 @@ int main(int argc, char* argv[]) { if (b_collisionDevelopment) { simulation.set_particle_interaction_buffers(particles, cells, kernels.m_cellAssignmentKernel, - kernels.m_rotationKernel); - - simulation.set_particle_step_buffers(particles, cells, - kernels.m_particleStepKernel); + kernels.m_rotationKernel); simulation.set_particle_bounce_buffers(particles, cells, kernels.m_particleBounceKernel); - - } else if (config.bubbleInteractionsOn) { // Set up buffers for GPU + } else { simulation.set_particle_step_buffers(particles, bubble, *stepKernel); } @@ -311,20 +289,16 @@ int main(int argc, char* argv[]) { // auto streamEndTime = high_resolution_clock::now(); // auto streamStartTime = high_resolution_clock::now(); - for (int i = 1; (simulation.getTime() <= config.maxTime); i++) { - if (config.m_maxSteps > 0 && simulation.getStep() > config.m_maxSteps) { - break; - } - /* - simulation.step(bubble, 0); - }*/ + for (int i = 1; + (simulation.getTime() <= config.maxTime) && + (config.m_maxSteps > 0 && simulation.getStep() < config.m_maxSteps); + i++) { + if (b_collisionDevelopment) { simulation.step(particles, cells, rn_generator, i, *stepKernel, kernels.m_cellAssignmentKernel, kernels.m_rotationKernel, kernels.m_particleBounceKernel, kernels.getCommandQueue()); - - particles.readParticlesBuffer(kernels.getCommandQueue()); } else { if (config.bubbleInteractionsOn) { simulation.step(particles, bubble, *stepKernel, diff --git a/configs/config.json b/configs/config.json index 281eae7..c7e112d 100644 --- a/configs/config.json +++ b/configs/config.json @@ -1,19 +1,19 @@ { "kernel":{ - "name": "particle_bubble_step_cyclic" + "name": "particle_step" }, "simulation": { "seed": 1, - "max_steps": 1500, + "max_steps": 10000, "dt": 0.001, - "max_time": 1.5, + "max_time": 15, "cyclic_boundary_on": true, "cyclic_boundary_radius": 500 }, "parameters": { "alpha": 12.24727746868043, - "eta": 17.3898723, + "eta": 100, "upsilon": 0.01, "coupling": 0.0, "dV": 21670.590989998676 @@ -25,26 +25,27 @@ "interaction_on": true }, "particles": { - "mass_false": 0.1, - "mass_true": 0.05576625911018531, - "T_false": 0.003206823957539085, + "mass_false": 0, + "mass_true": 1, + "T_false": 0.01, "T_true": 0.0, "N_false": 750000, "N_true": 0 }, "collision": { "collision_on": false, - "N_cells": 41, - "cell_length": 0.5 + "N_cells": 71, + "cell_length": 0.03 }, "stream": { "stream": true, - "stream_time": 0.01, - "stream_step": 100, + "stream_time": 0.05, + "stream_step": 500, "data_save_path": "data", "stream_data": true, "stream_density_profile": true, "steam_energy_profile": true, + "stream_momentum": false, "stream_momentumIn_profile": true, "stream_momentumOut_profile": true, "stream_radial_velocity": true, @@ -53,13 +54,13 @@ "bins_count_energy": 500, "bins_count_radial_velocity": 500, "bins_count_tangential_velocity": 500, - "bins_count_momentumIn": 500, - "bins_count_momentumOut": 500, + "bins_count_momentumIn": 1000, + "bins_count_momentumOut": 1000, "max_value_density": 1.1, "max_value_energy": 1.1, "max_value_radial_velocity": 1.1, "max_value_tangential_velocity": 1.1, - "max_value_momentumIn": 0.2, - "max_value_momentumOut": 0.2 + "max_value_momentumIn": 0.45, + "max_value_momentumOut": 0.45 } } \ No newline at end of file diff --git a/configs/test.json b/configs/test.json index e1d86a5..1515352 100644 --- a/configs/test.json +++ b/configs/test.json @@ -1,6 +1,6 @@ { "kernel":{ - "name": "particle_bubble_step_cyclic" + "name": "particle_step" }, "simulation": { "seed": 1, @@ -26,7 +26,7 @@ "particles": { "mass_false": 0, "mass_true": 1, - "T_false": 0.0, + "T_false": 0.5, "T_true": 0.0, "N_false": 5000, "N_true": 0 diff --git a/configs/test_oclgrind.json b/configs/test_oclgrind.json index 4b66aee..8b7167c 100644 --- a/configs/test_oclgrind.json +++ b/configs/test_oclgrind.json @@ -37,7 +37,7 @@ "cell_length": 0.5 }, "stream": { - "stream": true, + "stream": false, "stream_time": 0.01, "stream_step": 100, "data_save_path": "data", diff --git a/kernels/kernel.cl b/kernels/kernel.cl index a9dfbe5..62c2c9f 100644 --- a/kernels/kernel.cl +++ b/kernels/kernel.cl @@ -28,6 +28,7 @@ typedef struct Bubble { // In development typedef struct CollisionCell { + double gamma; double vX; double vY; double vZ; @@ -43,7 +44,6 @@ typedef struct CollisionCell { double p_z; double v2; // v2 = Sum: v_i^2 - double gamma; double mass; unsigned int particle_count; } CollisionCell; @@ -913,7 +913,7 @@ __kernel void particles_with_false_bubble_step_reflect( } -__kernel void assign_cell_index_to_particle( +__kernel void assign_particle_to_collision_cell( __global Particle *t_particles, __global const unsigned int *maxCellIndex, __global const double *cellLength, @@ -1000,11 +1000,10 @@ __kernel void label_particles_position_by_mass( t_particles[gid].b_inBubble = t_particles[gid].m == mass_in[0]; } -__kernel void transform_momentum( +__kernel void rotate_momentum( __global Particle *t_particles, __global CollisionCell *t_cells, __global unsigned int *number_of_cells - ){ unsigned int gid = get_global_id(0); Particle particle = t_particles[gid]; @@ -1018,8 +1017,6 @@ __kernel void transform_momentum( double sin_theta = sin(cell.theta); double p0, p1, p2, p3; - - if ((cell.particle_count > 1) && (cell.mass != 0)){ // Lorentz boost @@ -1096,7 +1093,7 @@ __kernel void transform_momentum( } } -__kernel void particle_bounce( +__kernel void particle_boundary_check( __global Particle *t_particles, __global double *boundaryDistanceFromCenter // [x_delta] ){