diff --git a/.gitignore b/.gitignore index ffa735c..fa71050 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,8 @@ dependencies/OpenCL bubbleSim/data +# CLtracer files (used for profiling OpenCL) +*.cltracer # User-specific files *.rsuser diff --git a/bubbleSim/bubbleSim.vcxproj.filters b/bubbleSim/bubbleSim.vcxproj.filters index a988896..535d8e5 100644 --- a/bubbleSim/bubbleSim.vcxproj.filters +++ b/bubbleSim/bubbleSim.vcxproj.filters @@ -1,14 +1,6 @@  - - {4FC737F1-C7A5-4376-A066-2A32D752A2FF} - cpp;c;cc;cxx;c++;cppm;ixx;def;odl;idl;hpj;bat;asm;asmx - - - {93995380-89BD-4b04-88EB-625FBE52EBFB} - h;hh;hpp;hxx;h++;hm;inl;inc;ipp;xsd - {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms @@ -16,81 +8,90 @@ {ad3d76ba-6f94-404f-bc5d-04c6fe891d8f} - - {ead60fb0-b85b-4979-a0a9-ef637c5cf283} + + {d8b60eea-a10a-4d44-93e7-1c97ff73bee3} + + + {95515f3a-561d-4459-b2c3-c2100ac20678} - - {0187af86-67e0-4ada-9316-2c50a6d1db91} + + {0a343cca-4fdc-401b-a415-4dcaa58ee610} - - {70c9ae99-7764-4421-ad94-5ac76b6b096b} + + {8d4d525e-9435-47ad-9eb5-1de1e4fb1975} - - {a5c559dd-9d80-4aef-a16b-8988be73efcd} + + {3c2ad733-df12-4a24-b5ff-a4632821a808} - - {729ee2fd-1b8b-46ac-87ea-9dccda36f30d} + + {01d28b61-e628-46b8-b3f5-1a908e45d035} + + + {737877a4-c2eb-4180-b746-d325876d14d1} + + + {f596f1c2-9307-4993-9c65-9f20f7cfe7bc} - - Header Files - - - Header Files - - - Header Files - - - Header Files + + Physics\Particle - Header Files\Objects - - - Header Files\Objects + Physics\Bubble - Header Files\Objects + Physics - - Source Files\Generators + + Physics\Dynamics + + + Stream + + + OpenCL - Header Files + Config - - Header Files + + Code - - Header Files\Streamers + + Code + + + Code - Header Files + Code + + + Code - - Source Files - - - Source Files - - - Source Files + + Physics\Particle - Source Files\Objects - - - Source Files\Objects + Physics\Bubble - Source Files + Physics\Dynamics - Source Files\Streamers + Stream + + + OpenCL + + + Code + + + Code diff --git a/bubbleSim/collision.cpp b/bubbleSim/collision.cpp index fb64f8b..a6814a6 100644 --- a/bubbleSim/collision.cpp +++ b/bubbleSim/collision.cpp @@ -5,22 +5,23 @@ CollisionCellCollection::CollisionCellCollection( bool t_doubleCellCount, cl::Context& cl_context) { int openCLerrNum; + /* TODO + * Add label to particles which distinguises if particle is inside or outside the bubble -> collision cell count gets 2x bigger + * + */ + m_doubleCellCount = t_doubleCellCount; + if (t_cellCountInOneAxis % 2 != 1) { std::cerr << "Cell count in one axis must be odd. (" << t_cellCountInOneAxis << ")" << std::endl; std::terminate(); } - - if (t_doubleCellCount) { - m_cellCount = 2 * (unsigned int)std::pow(t_cellCountInOneAxis, 3) + 1; - } else { - m_cellCount = (unsigned int)std::pow(t_cellCountInOneAxis, 3) + 1; - } + m_cellCount = (unsigned int)std::pow(t_cellCountInOneAxis, 3) + 1; m_cellCountBuffer = cl::Buffer(cl_context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(unsigned int), &m_cellCount, &openCLerrNum); - m_doubleCellCount = t_doubleCellCount; + m_meanFreePath = t_meanFreePath; m_cellCountInOneAxis = t_cellCountInOneAxis; @@ -52,29 +53,29 @@ 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 - double phi, theta; + numType phi, theta; - std::vector> frames; - frames.resize(m_cellCount); + std::vector> frames(m_cellCount, { (cl_numType)0,(cl_numType)0,(cl_numType)0,(cl_numType)0,(cl_numType)0,(cl_numType)1 }); /* * 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 @@ -82,16 +83,32 @@ void CollisionCellCollection::recalculate_cells( for (size_t i = 1; i < m_cellCount; i++) { m_collisionCells[i].particle_count = (int)frames[i][4]; if (frames[i][4] > 1) { + /*if (frames[i][4] != 2) { + m_collisionCells[i].particle_count = (int)0; + continue; + } + if (t_rng.generate_number() >= 0.01*0.01 / (9*frames[i][5])) { + m_collisionCells[i].particle_count = (int)0; + continue; + } + */ + + if (t_rng.generate_number() >= std::pow(0.01, frames[i][4]) / (9*frames[i][5])) { + m_collisionCells[i].particle_count = (int)0; + continue; + } + + 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]; + 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]; /* * beta = sqrt(v_x^2 + v_y^2 + v_z^2); @@ -104,16 +121,77 @@ void CollisionCellCollection::recalculate_cells( * 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(); + // Generate parameters for rotation axis + phi = std::acos((numType)1. - (numType)2. * t_rng.generate_number()); + theta = (numType)2. * (numType)M_PI * t_rng.generate_number(); + + // Define rotation angle and axis + 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); } } +} + +void CollisionCellCollection::recalculate_cells2( + std::vector& t_particles, RandomNumberGenerator& t_rng) { + // 0 index cell is for particles outside the collision cell grid + + numType phi, theta; + + std::vector> frames; + frames.resize(m_cellCount); + /* + * Sum up all momentums, energies and masses for each cell + */ + for (Particle& particle : t_particles) { + if (particle.idxCollisionCell == 0) continue; + 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; + } + /* + * 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].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]; + + /* + * 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); + + // Generate parameters for rotation axis + phi = std::acos((numType)1. - (numType)2. * t_rng.generate_number()); + theta = (numType)2. * (numType)M_PI * t_rng.generate_number(); + + // Define rotation angle and axis + 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..ce73e88 100644 --- a/bubbleSim/collision.h +++ b/bubbleSim/collision.h @@ -2,6 +2,7 @@ #include "base.h" #include "objects.h" typedef struct CollisionCell { + cl_numType gamma; cl_numType v_x; cl_numType v_y; cl_numType v_z; @@ -12,13 +13,12 @@ typedef struct CollisionCell { cl_numType theta; - cl_numType p_E; - cl_numType p_x; - cl_numType p_y; - cl_numType p_z; + 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,9 +30,14 @@ class CollisionCellCollection { std::array& getShiftVector() { return m_shiftVector; } + unsigned int getCellCount() { return m_cellCount; } + void recalculate_cells(std::vector& t_particles, RandomNumberGenerator& t_rng); + void recalculate_cells2(std::vector& t_particles, + RandomNumberGenerator& t_rng); + void generateShiftVector(RandomNumberGenerator& t_rng); cl::Buffer& getCellBuffer() { return m_collisionCellsBuffer; } diff --git a/bubbleSim/datastreamer.cpp b/bubbleSim/datastreamer.cpp index 9fd885d..a6db1c2 100644 --- a/bubbleSim/datastreamer.cpp +++ b/bubbleSim/datastreamer.cpp @@ -26,6 +26,34 @@ 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 +198,10 @@ 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 +216,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 +259,12 @@ 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; } @@ -233,8 +276,8 @@ void DataStreamer::stream(Simulation& simulation, if (m_initialized_RadialVelocity && (particleRadius < m_maxRadius_RadialVelocity)) { // Average A(N) = [A(N-1) * (N-1) + a(N) ]/N - particleRadialVelocity = - particleCollection.calculateParticleRadialVelocity(i); + particleRadialVelocity = (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); + //particleCollection.calculateParticleRadialVelocity(i); bins_RadialVelocity[(int)(particleRadius / m_dr_RadialVelocity)] = (bins_RadialVelocity[(int)(particleRadius / m_dr_RadialVelocity)] * bins_RadialVelocityCount[(int)(particleRadius / @@ -249,8 +292,13 @@ 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 +400,19 @@ 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( diff --git a/bubbleSim/datastreamer.h b/bubbleSim/datastreamer.h index 3faff42..872cf79 100644 --- a/bubbleSim/datastreamer.h +++ b/bubbleSim/datastreamer.h @@ -39,7 +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); @@ -85,6 +85,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; @@ -95,9 +98,11 @@ class DataStreamer { std::ofstream m_stream_RadialVelocity; std::ofstream m_stream_TangentialVelocity; + // 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,8 @@ 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 +122,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 +130,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..f9c8ba7 100644 --- a/bubbleSim/opencl_kernels.cpp +++ b/bubbleSim/opencl_kernels.cpp @@ -2,10 +2,10 @@ 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 cellAssignKernelName = "assign_particle_to_collision_cell"; + std::string transformKernelName = "collide_particles"; std::string particleStepKernelName = "particle_step"; - std::string particleBounceKernelName = "particle_bounce"; + std::string particleBounceKernelName = "particle_boundary_check"; std::string particleBubbleBoundaryStepKernelName = "particle_bubble_step_cyclic"; @@ -28,6 +28,10 @@ OpenCLLoader::OpenCLLoader(std::string kernelPath, std::string kernelName) { std::string particleBubbleStepKernelName = "particle_bubble_step"; std::string particleBubbleBoundaryStepKernelName = "particle_bubble_step_cyclic"; + // NB: Collision: change + std::string cellAssignKernelName = "assign_particle_to_collision_cell"; + std::string transformKernelName = "collide_particles"; + std::string particleBounceKernelName = "particle_boundary_check"; createContext(m_devices); createProgram(m_context, m_deviceUsed, kernelPath); @@ -37,6 +41,11 @@ OpenCLLoader::OpenCLLoader(std::string kernelPath, std::string kernelName) { createKernel(m_program, m_particleBubbleBoundaryStepKernel, particleBubbleBoundaryStepKernelName.c_str()); createQueue(m_context, m_deviceUsed); + + createKernel(m_program, m_rotationKernel, transformKernelName.c_str()); + createKernel(m_program, m_cellAssignmentKernel, cellAssignKernelName.c_str()); + createKernel(m_program, m_particleBounceKernel, + particleBounceKernelName.c_str()); } void OpenCLLoader::createContext(std::vector& devices) { diff --git a/bubbleSim/particle.cpp b/bubbleSim/particle.cpp index 0c311e8..1835cf4 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)); @@ -88,8 +88,8 @@ 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(); + 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 @@ -127,8 +127,8 @@ void ParticleGenerator::generatePointInSphere( numType& x, numType& y, numType& z, numType t_maxRadius, RandomNumberGenerator& t_generator) { numType phi = - std::acos(1 - 2 * t_generator.generate_number()); // inclination - numType theta = 2 * M_PI * t_generator.generate_number(); + std::acos((numType)1. - (numType)2. * t_generator.generate_number()); // inclination + numType theta = (numType)2. * (numType)M_PI * t_generator.generate_number(); numType radius = std::cbrt(t_generator.generate_number()) * t_maxRadius; x = radius * std::sin(phi) * std::cos(theta); // x y = radius * std::sin(phi) * std::sin(theta); // y @@ -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++) { @@ -177,7 +177,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}); } @@ -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,7 @@ 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 +352,7 @@ ParticleCollection::ParticleCollection( m_particleCountTotal * sizeof(Particle), m_particles.data(), &openCLerrNum); - m_dP = std::vector(m_particleCountTotal, 0.); + m_dP = std::vector(m_particleCountTotal, 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 +393,28 @@ 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(1 - std::pow(calculateParticleRadialVelocity(i), (numType)2.)); } numType ParticleCollection::calculateNumberDensity(numType t_mass, @@ -423,13 +423,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)(2. * std::pow(M_PI, (numType)2.)); return n; } @@ -439,16 +439,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)(2 * std::pow((numType)M_PI, (numType)2.)); return rho; } @@ -456,7 +456,7 @@ 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 +468,7 @@ 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,14 +481,14 @@ 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; } numType ParticleCollection::countParticleEnergyDensity(numType t_radius1, numType t_radius2) { numType energy = countParticlesEnergy(t_radius1, t_radius2); - 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.; return (numType)energy / volume; } @@ -528,8 +528,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..d539051 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..972fdd2 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -217,31 +217,63 @@ void Simulation::step(ParticleCollection& particles, m_time += m_dt; // Move particles cl_queue.enqueueNDRangeKernel(t_particleStepKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); + cl::NDRange(particles.getParticleCountTotal())); + cl_queue.enqueueNDRangeKernel(t_particleBounceKernel, cl::NullRange, + cl::NDRange(particles.getParticleCountTotal())); + + // std::array particleIdx = { 84522, 355324 , 429220, 538040, 704867 }; + // Generate shift vector if (i % 1 == 0) { - cells.generateShiftVector(generator_collision); - cells.writeShiftVectorBuffer(cl_queue); - - // Assign particles to collision cells - 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 - cl_queue.enqueueNDRangeKernel( - t_rotationKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); + cells.generateShiftVector(generator_collision); + cells.writeShiftVectorBuffer(cl_queue); + + // Assign particles to collision cells + 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); + + /*int cell_idx = 0; + for (unsigned int i = 0; i < cells.getCellCount(); i++) { + if ((cells.getCollisionCells()[i].particle_count >= 5) && (cells.getCollisionCells()[i].particle_count < 10)) { + cell_idx = i; + std::cout << "Found a cell " << std::endl; + break; + } + } + + CollisionCell cell = cells.getCollisionCells()[cell_idx]; + Particle p; + std::cout << "Cell idx: " << cell_idx << ", Gamma: " << cell.gamma << ", " << cell.v_x << ", " << cell.v_y << ", " << cell.v_z << std::endl; + */ + /*for (int i=0; i = 6 T^4 + particleGenerator1 = ParticleGenerator(config.particleMassFalse, + 3*config.particleTemperatureFalse); } else { particleGenerator1 = - ParticleGenerator(config.particleMassFalse, temperatureFalse, - 30 * temperatureFalse, 1e-5 * temperatureFalse); + 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.particleTemperatureFalse, config.particleCountTrue, config.particleCountFalse, config.parameterCoupling, config.bubbleIsTrueVacuum, kernels.getContext()); @@ -136,13 +141,14 @@ int main(int argc, char* argv[]) { numType genreatedParticleEnergy; if (b_collisionDevelopment) { genreatedParticleEnergy = particleGenerator1.generateNParticlesInBox( - config.bubbleInitialRadius, (u_int)particles.getParticleCountTotal(), + std::sqrt(3) * config.bubbleInitialRadius, (u_int)particles.getParticleCountTotal(), rn_generator, particles.getParticles()); } else { genreatedParticleEnergy = particleGenerator1.generateNParticlesInSphere( config.bubbleInitialRadius, (u_int)particles.getParticleCountTotal(), rn_generator, particles.getParticles()); } + numType total_energy = 0; for (unsigned int i = 0; i < config.particleCountFalse; i++) { total_energy += particles.getParticleEnergy(i); @@ -216,8 +222,7 @@ int main(int argc, char* argv[]) { kernels.m_cellAssignmentKernel, kernels.m_rotationKernel); - simulation.set_particle_step_buffers(particles, cells, - kernels.m_particleStepKernel); + simulation.set_particle_step_buffers(particles, cells, *stepKernel); simulation.set_particle_bounce_buffers(particles, cells, kernels.m_particleBounceKernel); @@ -238,8 +243,8 @@ int main(int argc, char* argv[]) { std::filesystem::path filePath = createSimulationFilePath(config.m_dataSavePath, dataFolderName); + DataStreamer streamer(filePath.string()); - if (config.streamDataOn) { streamer.initStream_Data(); } @@ -267,8 +272,13 @@ int main(int argc, char* argv[]) { streamer.initStream_MomentumOut(config.binsCountMomentumOut, config.maxValueMomentumOut); } + streamer.initStream_Momentum(config.binsCountMomentumIn, + config.maxValueMomentumIn); + streamer.stream(simulation, particles, bubble, kernels.getCommandQueue()); + + /* * =============== Display text =============== */ @@ -290,10 +300,23 @@ int main(int argc, char* argv[]) { initialNumberDensityFalse, initialEnergyDensityFalse); + bool streamCellCount = false; + std::ofstream pLocStream; + std::vector pLocArray; + if (streamCellCount) { + std::ofstream pLocStream(filePath / "cellIdx.csv"); + pLocArray.resize(config.collisionCellCount * config.collisionCellCount * config.collisionCellCount + 1, 0); + } /* =============== Run simulation =============== */ + /*std::array particleIdx = { 84522, 355324 , 429220, 538040, 704867 }; + for (int i : particleIdx) { + particles.printParticleInfo(i); + } +*/ + numType simTimeSinceLastStream = 0.; int stepsSinceLastStream = 0; std::cout << "=============== Simulation ===============" << std::endl; @@ -308,23 +331,30 @@ int main(int argc, char* argv[]) { << ", E: " << simulation.getTotalEnergy() / simulation.getInitialTotalEnergy() << std::endl; - + double px = 0; + double py = 0; + double pz = 0; // 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); - }*/ 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()); + simulation.step(particles, cells, rn_generator, i, *stepKernel, + kernels.m_cellAssignmentKernel, kernels.m_rotationKernel, + kernels.m_particleBounceKernel, + kernels.getCommandQueue()); + if (streamCellCount){ + for (Particle& p : particles.getParticles()) { + pLocArray[p.idxCollisionCell] += 1; + } + for (unsigned int loc : pLocArray) { + pLocStream << loc << ","; + } + pLocStream << "\n"; + std::fill(pLocArray.begin(), pLocArray.end(), 0); + } } else { if (config.bubbleInteractionsOn) { simulation.step(particles, bubble, *stepKernel, diff --git a/configs/config.json b/configs/config.json index 281eae7..8c4c1a2 100644 --- a/configs/config.json +++ b/configs/config.json @@ -5,7 +5,7 @@ }, "simulation": { "seed": 1, - "max_steps": 1500, + "max_steps": 100, "dt": 0.001, "max_time": 1.5, "cyclic_boundary_on": true, diff --git a/configs/test_oclgrind.json b/configs/test_oclgrind.json index 4b66aee..96fb5ba 100644 --- a/configs/test_oclgrind.json +++ b/configs/test_oclgrind.json @@ -6,7 +6,7 @@ "seed": 1, "dt": 0.1, "max_steps": 100, - "max_time": 150, + "max_time": 1, "cyclic_boundary_on": true, "cyclic_boundary_radius": 500 }, @@ -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..b167474 100644 --- a/kernels/kernel.cl +++ b/kernels/kernel.cl @@ -1,6 +1,6 @@ #pragma OPENCL EXTENSION cl_khr_fp64 : enable -// Particle type (also defined in c++ code) +// Particle struct (also defined in C++ code) typedef struct Particle { double x; double y; @@ -16,7 +16,8 @@ typedef struct Particle { char b_inBubble; int idxCollisionCell; } Particle; -// Bubble type (also defined in c++ code) + +// Bubble struct (also defined in C++ code) typedef struct Bubble { double radius; double radius2; // Squared @@ -28,22 +29,24 @@ typedef struct Bubble { // In development typedef struct CollisionCell { + // + double gamma; double vX; double vY; double vZ; + // Rotation vector double x; double y; double z; double theta; - double p_E; - double p_x; - double p_y; - double p_z; + double pE; + double pX; + double pY; + double pZ; double v2; // v2 = Sum: v_i^2 - double gamma; double mass; unsigned int particle_count; } CollisionCell; @@ -148,7 +151,7 @@ double calculateEnergy(Particle particle){ } __kernel void particle_bubble_step( - __global Particle *t_particles, + __global Particle *particles, __global double *t_dP, __global char *t_interactedFalse, __global char *t_passedFalse, @@ -169,7 +172,7 @@ __kernel void particle_bubble_step( // Bubble parameters struct Bubble bubble = t_bubble[0]; // Particle - struct Particle particle = t_particles[gid]; + struct Particle particle = particles[gid]; double M_in = t_m_in[0]; double M_out = t_m_out[0]; @@ -320,18 +323,18 @@ __kernel void particle_bubble_step( } // Update particle information - t_particles[gid].x = particle.x; - t_particles[gid].y = particle.y; - t_particles[gid].z = particle.z; - t_particles[gid].pX = particle.pX; - t_particles[gid].pY = particle.pY; - t_particles[gid].pZ = particle.pZ; - t_particles[gid].E = particle.E; - t_particles[gid].m = particle.m; + particles[gid].x = particle.x; + particles[gid].y = particle.y; + particles[gid].z = particle.z; + particles[gid].pX = particle.pX; + particles[gid].pY = particle.pY; + particles[gid].pZ = particle.pZ; + particles[gid].E = particle.E; + particles[gid].m = particle.m; } __kernel void particle_bubble_step_cyclic( - __global Particle *t_particles, + __global Particle *particles, __global double *t_dP, __global char *t_interactedFalse, __global char *t_passedFalse, @@ -341,7 +344,7 @@ __kernel void particle_bubble_step_cyclic( __constant double *t_m_in, __constant double *t_m_out, __constant double *t_delta_m2, - __constant double *t_cycleRadius + __constant double *boundaryRadius ){ unsigned int gid = get_global_id(0); @@ -351,7 +354,7 @@ __kernel void particle_bubble_step_cyclic( // Bubble parameters struct Bubble bubble = t_bubble[0]; // Particle - struct Particle particle = t_particles[gid]; + struct Particle particle = particles[gid]; double M_in = t_m_in[0]; double M_out = t_m_out[0]; @@ -534,30 +537,30 @@ __kernel void particle_bubble_step_cyclic( * If particle.x < -R_cyclic -> Then particle.x = particle.x + 2 * R_cyclic */ - double cyclicRadius = t_cycleRadius[0]; - particle.x = (cyclicRadius > fabs(particle.x)) * particle.x + - (particle.x > cyclicRadius) * (particle.x - 2*cyclicRadius) + - (particle.x < -cyclicRadius) * (particle.x + 2*cyclicRadius); - particle.y = (cyclicRadius > fabs(particle.y)) * particle.y + - (particle.y > cyclicRadius) * (particle.y - 2*cyclicRadius) + - (particle.y < -cyclicRadius) * (particle.y + 2*cyclicRadius); - particle.z = (cyclicRadius > fabs(particle.z)) * particle.z + - (particle.z > cyclicRadius) * (particle.z - 2*cyclicRadius) + - (particle.z < -cyclicRadius) * (particle.z + 2*cyclicRadius); + double r_boundaryRadius = boundaryRadius[0]; + particle.x = (r_boundaryRadius > fabs(particle.x)) * particle.x + + (particle.x > r_boundaryRadius) * (particle.x - 2*r_boundaryRadius) + + (particle.x < -r_boundaryRadius) * (particle.x + 2*r_boundaryRadius); + particle.y = (r_boundaryRadius > fabs(particle.y)) * particle.y + + (particle.y > r_boundaryRadius) * (particle.y - 2*r_boundaryRadius) + + (particle.y < -r_boundaryRadius) * (particle.y + 2*r_boundaryRadius); + particle.z = (r_boundaryRadius > fabs(particle.z)) * particle.z + + (particle.z > r_boundaryRadius) * (particle.z - 2*r_boundaryRadius) + + (particle.z < -r_boundaryRadius) * (particle.z + 2*r_boundaryRadius); - t_particles[gid].x = particle.x; - t_particles[gid].y = particle.y; - t_particles[gid].z = particle.z; + particles[gid].x = particle.x; + particles[gid].y = particle.y; + particles[gid].z = particle.z; - t_particles[gid].pX = particle.pX; - t_particles[gid].pY = particle.pY; - t_particles[gid].pZ = particle.pZ; - t_particles[gid].E = particle.E; - t_particles[gid].m = particle.m; + particles[gid].pX = particle.pX; + particles[gid].pY = particle.pY; + particles[gid].pZ = particle.pZ; + particles[gid].E = particle.E; + particles[gid].m = particle.m; } __kernel void particle_bubble_step_cyclic_mass_inverted( - __global Particle *t_particles, + __global Particle *particles, __global double *t_dP, __global char *t_interactedFalse, __global char *t_passedFalse, @@ -567,7 +570,7 @@ __kernel void particle_bubble_step_cyclic_mass_inverted( __constant double *t_m_in, __constant double *t_m_out, __constant double *t_delta_m2, - __constant double *t_cycleRadius + __constant double *boundaryRadius ){ unsigned int gid = get_global_id(0); @@ -577,7 +580,7 @@ __kernel void particle_bubble_step_cyclic_mass_inverted( // Bubble parameters struct Bubble bubble = t_bubble[0]; // Particle - struct Particle particle = t_particles[gid]; + struct Particle particle = particles[gid]; double M_in = t_m_in[0]; double M_out = t_m_out[0]; @@ -766,54 +769,53 @@ __kernel void particle_bubble_step_cyclic_mass_inverted( * If particle.x < -R_cyclic -> Then particle.x = particle.x + 2 * R_cyclic */ - double cyclicRadius = t_cycleRadius[0]; - particle.x = (cyclicRadius > fabs(particle.x)) * particle.x + - (particle.x > cyclicRadius) * (particle.x - 2*cyclicRadius) + - (particle.x < -cyclicRadius) * (particle.x + 2*cyclicRadius); - particle.y = (cyclicRadius > fabs(particle.y)) * particle.y + - (particle.y > cyclicRadius) * (particle.y - 2*cyclicRadius) + - (particle.y < -cyclicRadius) * (particle.y + 2*cyclicRadius); - particle.z = (cyclicRadius > fabs(particle.z)) * particle.z + - (particle.z > cyclicRadius) * (particle.z - 2*cyclicRadius) + - (particle.z < -cyclicRadius) * (particle.z + 2*cyclicRadius); + double r_boundaryRadius = boundaryRadius[0]; + particle.x = (r_boundaryRadius > fabs(particle.x)) * particle.x + + (particle.x > r_boundaryRadius) * (particle.x - 2*r_boundaryRadius) + + (particle.x < -r_boundaryRadius) * (particle.x + 2*r_boundaryRadius); + particle.y = (r_boundaryRadius > fabs(particle.y)) * particle.y + + (particle.y > r_boundaryRadius) * (particle.y - 2*r_boundaryRadius) + + (particle.y < -r_boundaryRadius) * (particle.y + 2*r_boundaryRadius); + particle.z = (r_boundaryRadius > fabs(particle.z)) * particle.z + + (particle.z > r_boundaryRadius) * (particle.z - 2*r_boundaryRadius) + + (particle.z < -r_boundaryRadius) * (particle.z + 2*r_boundaryRadius); - t_particles[gid].x = particle.x; - t_particles[gid].y = particle.y; - t_particles[gid].z = particle.z; + particles[gid].x = particle.x; + particles[gid].y = particle.y; + particles[gid].z = particle.z; - t_particles[gid].pX = particle.pX; - t_particles[gid].pY = particle.pY; - t_particles[gid].pZ = particle.pZ; - t_particles[gid].E = particle.E; - t_particles[gid].m = particle.m; + particles[gid].pX = particle.pX; + particles[gid].pY = particle.pY; + particles[gid].pZ = particle.pZ; + particles[gid].E = particle.E; + particles[gid].m = particle.m; } - __kernel void particle_step( - __global Particle *t_particles, - __constant double *boundaryDistanceFromCenter, + __global Particle *particles, + __constant double *boundaryRadius, __constant double *t_dt ){ unsigned int gid = get_global_id(0); - Particle particle = t_particles[gid]; - double Vx = particle.pX/particle.E; - double Vy = particle.pY/particle.E; - double Vz = particle.pZ/particle.E; - - double dt = t_dt[0]; + Particle particle = particles[gid]; + // double Vx = particle.pX/particle.E; + // double Vy = particle.pY/particle.E; + // double Vz = particle.pZ/particle.E; + // double dt = t_dt[0]; - moveLinear(&particle, Vx, Vy, Vz, dt); + moveLinear(&particle, particle.pX/particle.E, particle.pY/particle.E, particle.pZ/particle.E, t_dt[0]); - t_particles[gid] = particle; + + particles[gid] = particle; } __kernel void particles_with_false_bubble_step_reflect( - __global Particle *t_particles, __global double *t_dP, + __global Particle *particles, __global double *t_dP, __global char *t_interactedFalse, __global char *t_passedFalse, __global char *t_interactedTrue, __constant Bubble *t_bubble, __constant double *t_dt, __constant double *t_m_in, - __constant double *t_m_out, __constant double *t_delta_m2, __constant double *t_cycleRadius) { + __constant double *t_m_out, __constant double *t_delta_m2, __constant double *boundaryRadius) { unsigned int gid = get_global_id(0); // dE - dP is not actual energy difference. dE = ΔE/R_b -> to avoid // singularities/noise near R_b ~ 0 @@ -821,7 +823,7 @@ __kernel void particles_with_false_bubble_step_reflect( // Bubble parameters struct Bubble bubble = t_bubble[0]; // Particle - struct Particle particle = t_particles[gid]; + struct Particle particle = particles[gid]; double M_in = t_m_in[0]; double M_out = t_m_out[0]; @@ -901,116 +903,199 @@ __kernel void particles_with_false_bubble_step_reflect( particle.z = 100000; } - t_particles[gid].x = particle.x; - t_particles[gid].y = particle.y; - t_particles[gid].z = particle.z; + particles[gid].x = particle.x; + particles[gid].y = particle.y; + particles[gid].z = particle.z; - t_particles[gid].pX = particle.pX; - t_particles[gid].pY = particle.pY; - t_particles[gid].pZ = particle.pZ; - t_particles[gid].E = particle.E; - t_particles[gid].m = particle.m; + particles[gid].pX = particle.pX; + particles[gid].pY = particle.pY; + particles[gid].pZ = particle.pZ; + particles[gid].E = particle.E; + particles[gid].m = particle.m; } -__kernel void assign_cell_index_to_particle( - __global Particle *t_particles, +__kernel void assign_particle_to_collision_cell( + __global Particle *particles, __global const unsigned int *maxCellIndex, __global const double *cellLength, __global const double *cuboidShift ){ unsigned int gid = get_global_id(0); - - Particle particle = t_particles[gid]; - // Find cell numbers - int x_index = (int) ((particle.x + cellLength[0]*maxCellIndex[0]/2 + cuboidShift[0]) / cellLength[0]); - int y_index = (int) ((particle.y + cellLength[0]*maxCellIndex[0]/2 + cuboidShift[1]) / cellLength[0]); - int z_index = (int) ((particle.z + cellLength[0]*maxCellIndex[0]/2 + cuboidShift[2]) / cellLength[0]); - // Idx = 0 -> if particle is outside of the cuboid cell structure - // Convert cell number into 1D vector - if ((x_index < 0) || (x_index >= maxCellIndex[0])){ - t_particles[gid].idxCollisionCell = 0; - } - else if ((y_index < 0) || (y_index >= maxCellIndex[0])){ - t_particles[gid].idxCollisionCell = 0; - } - else if ((z_index < 0) || (z_index >= maxCellIndex[0])){ - t_particles[gid].idxCollisionCell = 0; - } - else { - t_particles[gid].idxCollisionCell = 1 + x_index + y_index * maxCellIndex[0] + z_index * maxCellIndex[0] * maxCellIndex[0]; - } - - //t_particles[gid].idxCollisionCell = x_index + y_index + z_index; -} + // Read variables to local registers + Particle particle = particles[gid]; -__kernel void assign_particle_cell_index_two_phase( - __global Particle *t_particles, - __global const int *maxCellIndex, - __global const double *cellLength, - __global const double *cuboidShift // Random particle location shift - - ){ - unsigned int gid = get_global_id(0); + unsigned int r_maximumCellIndex = maxCellIndex[0]; + double r_cellLength = cellLength[0]; + + // Find cell number in 3D cell + int x_index = (int) ((particle.x + r_cellLength*r_maximumCellIndex/2 + cuboidShift[0]) / r_cellLength); + int y_index = (int) ((particle.y + r_cellLength*r_maximumCellIndex/2 + cuboidShift[1]) / r_cellLength); + int z_index = (int) ((particle.z + r_cellLength*r_maximumCellIndex/2 + cuboidShift[2]) / r_cellLength); - Particle particle = t_particles[gid]; - // Find cell numbers - int a = (int) ((particle.x + cuboidShift[0]) / cellLength[0]); - int b = (int) ((particle.y + cuboidShift[1]) / cellLength[1]); - int c = (int) ((particle.z + cuboidShift[2]) / cellLength[2]); - // Idx = 0 -> if particle is outside of the cuboid cell structure - // Convert cell number into 1D vector. First half of the vector is for outisde the bubble and second half is inside the bubble - if ((a < 0) || (a >= maxCellIndex[0])){ - t_particles[gid].idxCollisionCell = 0; + // Assign particles which don't fit to the cell structure at index = 0. + int isIndexZero = ((x_index < 0) || (x_index >= r_maximumCellIndex) || (y_index < 0) || (y_index >= r_maximumCellIndex) || (z_index < 0) || (z_index >= r_maximumCellIndex)); + // If index not zero -> !isIndexZero + particles[gid].idxCollisionCell = 0 + !isIndexZero*(1 + x_index + y_index * r_maximumCellIndex + z_index * r_maximumCellIndex * r_maximumCellIndex); + + // Old implementation + /* + if ((x_index < 0) || (x_index >= r_maximumCellIndex)){ + particles[gid].idxCollisionCell = 0; } - else if ((b < 0) || (b >= maxCellIndex[1])){ - t_particles[gid].idxCollisionCell = 0; + else if ((y_index < 0) || (y_index >= r_maximumCellIndex)){ + particles[gid].idxCollisionCell = 0; } - else if ((c < 0) || (c >= maxCellIndex[2])){ - t_particles[gid].idxCollisionCell = 0; + else if ((z_index < 0) || (z_index >= r_maximumCellIndex)){ + particles[gid].idxCollisionCell = 0; } else { - t_particles[gid].idxCollisionCell = 1 + a + b * maxCellIndex[0] + c * maxCellIndex[0] * maxCellIndex[1] + particle.b_inBubble * maxCellIndex[0] * maxCellIndex[1]*maxCellIndex[2]; + particles[gid].idxCollisionCell = 1 + x_index + y_index * r_maximumCellIndex + z_index * r_maximumCellIndex * r_maximumCellIndex; } - + */ } -__kernel void label_particles_position_by_coordinate( - __global Particle *t_particles, +__kernel void is_particle_in_bubble( + __global Particle *particles, __global Bubble *t_bubble ){ unsigned int gid = get_global_id(0); // If R_b^2 > R_x^2 then particle is inside the bubble - t_particles[gid].b_inBubble = fma( - t_particles[gid].x, t_particles[gid].x, - fma(t_particles[gid].y, t_particles[gid].y, - t_particles[gid].z * t_particles[gid].z )) < t_bubble[0].radius2; + particles[gid].b_inBubble = fma( + particles[gid].x, particles[gid].x, + fma(particles[gid].y, particles[gid].y, + particles[gid].z * particles[gid].z )) < t_bubble[0].radius2; } - -__kernel void label_particles_position_by_mass( - __global Particle *t_particles, - __global double *mass_in - ){ - unsigned int gid = get_global_id(0); +__kernel void collide_particles( + __global Particle *particles, + __global CollisionCell *cells, + __global unsigned int *number_of_cells - // If R_b^2 > R_x^2 then particle is inside the bubble - t_particles[gid].b_inBubble = t_particles[gid].m == mass_in[0]; + ){ + /* + * Collision algorith based on multi particle collision algorithm. + * + */ + unsigned int gid = get_global_id(0); + Particle particle = particles[gid]; + // If in bubble then cell number is doubled and second half is in bubble cells. + if (particle.idxCollisionCell != 0){ + CollisionCell cell = cells[particle.idxCollisionCell]; + if ((cell.particle_count > 1) && (cell.mass != 0)){ + //printf("siin\n"); + double gamma_minus_one = cell.gamma - 1; + double cos_theta = cos(cell.theta); + double one_minus_cos_theta = 1 - cos_theta; + double sin_theta = sin(cell.theta); + double p0, p1, p2, p3; + // Lorentz boost + //if (gid==0){ + // printf("p0: %.8f\np1: %.8f\np2: %.8f\np3: %.8f\nv1: %.8f\nv2: %.8f\nv3: %.8f\nn1: %.8f\nn2: %.8f\nn3: %.8f\ntheta: %.8f\n", + // particle.E, particle.pX, particle.pY, particle.pZ, cell.vX, cell.vY, cell.vZ, cell.x, cell.y, cell.z, cell.theta + // ); + //} + /* + * Lorentz transformation to zero momentum frame + */ + + p0 = cell.gamma * ( + particle.E + - particle.pX * cell.vX + - particle.pY * cell.vY + - particle.pZ * cell.vZ + ); + p1 = - particle.E * cell.gamma * cell.vX + + particle.pX * (1 + cell.vX * cell.vX * gamma_minus_one/cell.v2) + + particle.pY * cell.vX * cell.vY * gamma_minus_one/cell.v2 + + particle.pZ * cell.vX * cell.vZ * gamma_minus_one/cell.v2; + + p2 = -cell.gamma * particle.E * cell.vY + + particle.pY * (1 + cell.vY * cell.vY * gamma_minus_one/cell.v2) + + particle.pX * cell.vX * cell.vY * gamma_minus_one/cell.v2 + + particle.pZ * cell.vY * cell.vZ * gamma_minus_one/cell.v2; + + p3 = -cell.gamma * particle.E * cell.vZ + + particle.pZ * (1 + cell.vZ * cell.vZ * gamma_minus_one/cell.v2) + + particle.pX * cell.vX * cell.vZ * gamma_minus_one/cell.v2 + + particle.pY * cell.vY * cell.vZ * gamma_minus_one/cell.v2; + + particle.E = p0; + particle.pX = p1; + particle.pY = p2; + particle.pZ = p3; + + /* + * Rotate particle momentum + */ + p1 = particle.pX * (cos_theta + pow(cell.x, 2) * one_minus_cos_theta) + + particle.pY * (cell.x * cell.y * one_minus_cos_theta - cell.z * sin_theta) + + particle.pZ * (cell.x * cell.z * one_minus_cos_theta + cell.y * sin_theta); + p2 = particle.pX * (cell.x * cell.y * one_minus_cos_theta + cell.z * sin_theta) + + particle.pY * (cos_theta + pow(cell.y, 2) * one_minus_cos_theta) + + particle.pZ * (cell.y*cell.z*one_minus_cos_theta - cell.x * sin_theta); + p3 = particle.pX * (cell.x * cell.z * one_minus_cos_theta - cell.y * sin_theta) + + particle.pY * (cell.y * cell.z * one_minus_cos_theta + cell.x * sin_theta) + + particle.pZ * (cos_theta + pow(cell.z, 2) * one_minus_cos_theta); + + particle.pX = p1; + particle.pY = p2; + particle.pZ = p3; + + /* + * Lorentz transformation back to initial momentum frame + */ + p0 = cell.gamma * ( + particle.E + + particle.pX * cell.vX + + particle.pY * cell.vY + + particle.pZ * cell.vZ + ); + p1 = cell.gamma * particle.E * cell.vX + + particle.pX * (1 + cell.vX * cell.vX * gamma_minus_one/cell.v2) + + particle.pY * cell.vX * cell.vY * gamma_minus_one/cell.v2 + + particle.pZ * cell.vX * cell.vZ * gamma_minus_one/cell.v2; + + p2 = cell.gamma * particle.E * cell.vY + + particle.pY * (1 + cell.vY * cell.vY * gamma_minus_one/cell.v2) + + particle.pX * cell.vX * cell.vY * gamma_minus_one/cell.v2 + + particle.pZ * cell.vY * cell.vZ * gamma_minus_one/cell.v2; + + p3 = cell.gamma * particle.E * cell.vZ + + particle.pZ * (1 + cell.vZ * cell.vZ * gamma_minus_one/cell.v2) + + particle.pX * cell.vX * cell.vZ * gamma_minus_one/cell.v2 + + particle.pY * cell.vY * cell.vZ * gamma_minus_one/cell.v2; + + particle.E = p0; + particle.pX = p1; + particle.pY = p2; + particle.pZ = p3; + // if (gid==0){ + // printf("E: %.8f, pX: %.8f, pY: %.8f, pZ: %.8f\n", p0, p1, p2, p3); + // } + + particles[gid] = particle; + } + } } -__kernel void transform_momentum( - __global Particle *t_particles, - __global CollisionCell *t_cells, +__kernel void collide_particles2( + __global Particle *particles, + __global CollisionCell *cells, __global unsigned int *number_of_cells ){ + /* + * Collision algorith based on multi particle collision algorithm. + * + */ unsigned int gid = get_global_id(0); - Particle particle = t_particles[gid]; + Particle particle = particles[gid]; // If in bubble then cell number is doubled and second half is in bubble cells. if (particle.idxCollisionCell != 0){ - CollisionCell cell = t_cells[particle.idxCollisionCell]; + CollisionCell cell = cells[particle.idxCollisionCell]; double gamma_minus_one = cell.gamma - 1; double cos_theta = cos(cell.theta); @@ -1018,11 +1103,17 @@ __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 - + //if (gid==0){ + // printf("p0: %.8f\np1: %.8f\np2: %.8f\np3: %.8f\nv1: %.8f\nv2: %.8f\nv3: %.8f\nn1: %.8f\nn2: %.8f\nn3: %.8f\ntheta: %.8f\n", + // particle.E, particle.pX, particle.pY, particle.pZ, cell.vX, cell.vY, cell.vZ, cell.x, cell.y, cell.z, cell.theta + // ); + //} + /* + * Lorentz transformation to zero momentum frame + */ + p0 = cell.gamma * ( particle.E - particle.pX * cell.vX @@ -1049,7 +1140,9 @@ __kernel void transform_momentum( particle.pY = p2; particle.pZ = p3; - // Rotate momentum + /* + * Rotate particle momentum + */ p1 = particle.pX * (cos_theta + pow(cell.x, 2) * one_minus_cos_theta) + particle.pY * (cell.x * cell.y * one_minus_cos_theta - cell.z * sin_theta) + particle.pZ * (cell.x * cell.z * one_minus_cos_theta + cell.y * sin_theta); @@ -1064,7 +1157,9 @@ __kernel void transform_momentum( particle.pY = p2; particle.pZ = p3; - // Lorentz inverse transformation + /* + * Lorentz transformation back to initial momentum frame + */ p0 = cell.gamma * ( particle.E + particle.pX * cell.vX @@ -1090,35 +1185,43 @@ __kernel void transform_momentum( particle.pX = p1; particle.pY = p2; particle.pZ = p3; + // if (gid==0){ + // printf("E: %.8f, pX: %.8f, pY: %.8f, pZ: %.8f\n", p0, p1, p2, p3); + // } - t_particles[gid] = particle; + particles[gid] = particle; } } - } +} -__kernel void particle_bounce( - __global Particle *t_particles, - __global double *boundaryDistanceFromCenter // [x_delta] +__kernel void particle_boundary_check( + __global Particle *particles, + __global double *boundaryRadius // [x_delta] ){ + /* + * If particle is out of boundaries then update it's location. + * Update coordinate that particle goes to the other side of the simulation space. + */ unsigned int gid = get_global_id(0); - Particle particle = t_particles[gid]; - - // If Particle is inside the boundary leave value same. Otherwise change the sign - - // abs(x) < Boundary -> leave momentum - // abs(x) > Boundary and x < -Boundary -> Set momentum positive - // abs(x) > Boundary and x > Boundary -> Set momentum negative - particle.pX = ( boundaryDistanceFromCenter[0] > fabs(particle.x)) * particle.pX + // If inside, leave momentum alone - ((boundaryDistanceFromCenter[0] < fabs(particle.x)) && (-boundaryDistanceFromCenter[0] > particle.x)) * fabs(particle.pX) - // particle too -, chane to + - ((boundaryDistanceFromCenter[0] < fabs(particle.x)) && (boundaryDistanceFromCenter[0] < particle.x)) * fabs(particle.pX); // particle too +, chane to - - particle.pY = ( boundaryDistanceFromCenter[0] > fabs(particle.y)) * particle.pY + - ((boundaryDistanceFromCenter[0] < fabs(particle.y)) && (-boundaryDistanceFromCenter[0] > particle.y)) * fabs(particle.pY) - - ((boundaryDistanceFromCenter[0] < fabs(particle.y)) && (boundaryDistanceFromCenter[0] < particle.y)) * fabs(particle.pY); - particle.pZ = ( boundaryDistanceFromCenter[0] > fabs(particle.z)) * particle.pZ + - ((boundaryDistanceFromCenter[0] < fabs(particle.z)) && (-boundaryDistanceFromCenter[0] > particle.z)) * fabs(particle.pZ) - - ((boundaryDistanceFromCenter[0] < fabs(particle.z)) && (boundaryDistanceFromCenter[0] < particle.z)) * fabs(particle.pZ); + Particle particle = particles[gid]; + double r_boundaryRadius = boundaryRadius[0]; + double doubleBoundaryRadius = 2 * r_boundaryRadius; + + // If Particle is inside the boundary leave value same. Otherwise change the sign + // abs(x) < r_boundaryRadius -> leave just as it is + // x > r_boundaryRadius -> x - 2 * r_boundaryRadius + // x < r_boundaryRadius -> x + 2 * r_boundaryRadius + particle.x = (r_boundaryRadius > fabs(particle.x)) * particle.x + + (particle.x > r_boundaryRadius) * (particle.x - doubleBoundaryRadius) + + (particle.x < -r_boundaryRadius) * (particle.x + doubleBoundaryRadius); + particle.y = (r_boundaryRadius > fabs(particle.y)) * particle.y + + (particle.y > r_boundaryRadius) * (particle.y - doubleBoundaryRadius) + + (particle.y < -r_boundaryRadius) * (particle.y + doubleBoundaryRadius); + particle.z = (r_boundaryRadius > fabs(particle.z)) * particle.z + + (particle.z > r_boundaryRadius) * (particle.z - doubleBoundaryRadius) + + (particle.z < -r_boundaryRadius) * (particle.z + doubleBoundaryRadius); // Update result - t_particles[gid] = particle; + particles[gid] = particle; }