From fbcef1ebcd760aedafbdef9c358370181ee4d20f Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Mon, 5 Jun 2023 18:04:45 +0300 Subject: [PATCH 1/8] Preparation for collisions. Optimization --- .gitignore | 2 + bubbleSim/bubbleSim.vcxproj.filters | 113 +++++++++--------- bubbleSim/collision.cpp | 17 ++- bubbleSim/datastreamer.cpp | 13 ++- bubbleSim/opencl_kernels.cpp | 9 ++ bubbleSim/particle.cpp | 66 +++++------ bubbleSim/simulation.cpp | 11 +- bubbleSim/simulation.h | 2 +- bubbleSim/source.cpp | 40 +++++-- configs/config.json | 2 +- kernels/kernel.cl | 175 +++++++++++++++++++++++----- 11 files changed, 305 insertions(+), 145 deletions(-) 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..fc948d9 100644 --- a/bubbleSim/collision.cpp +++ b/bubbleSim/collision.cpp @@ -52,23 +52,23 @@ 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); /* * 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; @@ -110,10 +110,15 @@ void CollisionCellCollection::recalculate_cells( /* * Generate rotation axis and angle for CollisionCell */ - m_collisionCells[i].theta = 2 * M_PI * t_rng.generate_number(); + m_collisionCells[i].theta = -M_PI / 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); } } + + /*int particle1_idx = t_particles[0].idxCollisionCell; + std::cout << t_particles[0].E << ", " << t_particles[0].p_x << ", " + << t_particles[0].p_y << ", " << t_particles[0].p_z << std::endl; + std::cout << m_collisionCells[particle1_idx].theta << std::endl;*/ } \ No newline at end of file diff --git a/bubbleSim/datastreamer.cpp b/bubbleSim/datastreamer.cpp index 9fd885d..c9dac17 100644 --- a/bubbleSim/datastreamer.cpp +++ b/bubbleSim/datastreamer.cpp @@ -233,8 +233,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].p_x + particleCollection.getParticles()[i].y * particleCollection.getParticles()[i].p_y + particleCollection.getParticles()[i].z * particleCollection.getParticles()[i].p_z) / (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 +249,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)] * diff --git a/bubbleSim/opencl_kernels.cpp b/bubbleSim/opencl_kernels.cpp index 8e998b6..1c08812 100644 --- a/bubbleSim/opencl_kernels.cpp +++ b/bubbleSim/opencl_kernels.cpp @@ -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_cell_index_to_particle"; + std::string transformKernelName = "transform_momentum_massive"; + std::string particleBounceKernelName = "particle_bounce"; 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..a0e128d 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); @@ -414,7 +414,7 @@ numType ParticleCollection::calculateParticleRadialVelocity(u_int 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/simulation.cpp b/bubbleSim/simulation.cpp index d7af051..c63eab7 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -216,22 +216,17 @@ void Simulation::step(ParticleCollection& particles, cl::CommandQueue& cl_queue) { m_time += m_dt; // Move particles - cl_queue.enqueueNDRangeKernel(t_particleStepKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); + // Generate shift vector if (i % 1 == 0) { + particles.readParticlesBuffer(cl_queue); 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 @@ -242,6 +237,8 @@ void Simulation::step(ParticleCollection& particles, t_rotationKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); } + cl_queue.enqueueNDRangeKernel(t_particleStepKernel, cl::NullRange, + cl::NDRange(particles.getParticleCountTotal())); cl_queue.enqueueNDRangeKernel(t_particleBounceKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); } diff --git a/bubbleSim/simulation.h b/bubbleSim/simulation.h index 6b7548c..5f98657 100644 --- a/bubbleSim/simulation.h +++ b/bubbleSim/simulation.h @@ -51,7 +51,7 @@ class Simulation { } }; - int getStep() { return m_step; } + size_t getStep() { return m_step; } void step(PhaseBubble& bubble, numType t_dP); /* diff --git a/bubbleSim/source.cpp b/bubbleSim/source.cpp index 973aee4..42fe6ff 100644 --- a/bubbleSim/source.cpp +++ b/bubbleSim/source.cpp @@ -76,12 +76,14 @@ 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" std::string s_kernelPath = argv[2]; // "kernel.cl"; + s_configPath = "D:\\dev\\bubbleSim\\configs\\test_collision.json"; + std::cout << "Config path: " << s_configPath << std::endl; std::cout << "Kernel path: " << s_kernelPath << std::endl; @@ -91,7 +93,7 @@ int main(int argc, char* argv[]) { using std::chrono::milliseconds; auto programStartTime = high_resolution_clock::now(); ConfigReader config(s_configPath); - /* + /* =============== =============== */ // If seed = 0 then it generates random seed. @@ -119,7 +121,8 @@ int main(int argc, char* argv[]) { 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, @@ -143,6 +146,7 @@ int main(int argc, char* argv[]) { 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 +220,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 +241,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(); } @@ -269,6 +272,8 @@ int main(int argc, char* argv[]) { } streamer.stream(simulation, particles, bubble, kernels.getCommandQueue()); + + /* * =============== Display text =============== */ @@ -308,16 +313,15 @@ 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, @@ -325,6 +329,22 @@ int main(int argc, char* argv[]) { kernels.getCommandQueue()); particles.readParticlesBuffer(kernels.getCommandQueue()); + /*std::cout << particles.getParticles()[0].E << ", " + << particles.getParticles()[0].p_x << ", " + << particles.getParticles()[0].p_y << ", " + << particles.getParticles()[0].p_z << std::endl;*/ + /*px = 0; + py = 0; + pz = 0; + for (Particle& p : particles.getParticles()) { + px += p.p_x; + py += p.p_y; + pz += p.p_z; + } + std::cout << px / 750000. << std::endl; + std::cout << py / 750000. << std::endl; + std::cout << pz / 750000. << std::endl;*/ + } 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/kernels/kernel.cl b/kernels/kernel.cl index a9dfbe5..767abe3 100644 --- a/kernels/kernel.cl +++ b/kernels/kernel.cl @@ -788,23 +788,22 @@ __kernel void particle_bubble_step_cyclic_mass_inverted( t_particles[gid].m = particle.m; } - __kernel void particle_step( __global Particle *t_particles, - __constant double *boundaryDistanceFromCenter, + __constant double *t_cycleRadius, __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]; + // 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; } @@ -920,7 +919,6 @@ __kernel void assign_cell_index_to_particle( __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]); @@ -928,6 +926,7 @@ __kernel void assign_cell_index_to_particle( 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; } @@ -989,7 +988,6 @@ __kernel void label_particles_position_by_coordinate( t_particles[gid].z * t_particles[gid].z )) < t_bubble[0].radius2; } - __kernel void label_particles_position_by_mass( __global Particle *t_particles, __global double *mass_in @@ -1018,11 +1016,16 @@ __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 + // ); + //} + + p0 = cell.gamma * ( particle.E - particle.pX * cell.vX @@ -1090,35 +1093,153 @@ __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; } } - } +} + +__kernel void transform_momentum_massive( + __global Particle *t_particles, + __global CollisionCell *t_cells, + __global unsigned int *number_of_cells + ){ + unsigned int gid = get_global_id(0); + // Copy object to register memory. Improves performance. + Particle particle = t_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]; + //if (gid==0){ + // printf("Cell count: %i ,",cell.particle_count); + // } + 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; + double p0result, p1result, p2result, p3result; + double new_gamma; + double mgamma; + + if ((cell.particle_count > 1) && (cell.mass != 0)){ + // Lorentz boost + //if (gid==0){ + // printf("x: %.4f, %.4f, %.4f, cell: %i\n", particle.x, particle.y, particle.z, particle.idxCollisionCell); + //} + //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 + // ); + //} + + + + 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 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 inverse transformation + 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; + + + t_particles[gid] = particle; + } + } +} + + + + __kernel void particle_bounce( __global Particle *t_particles, - __global double *boundaryDistanceFromCenter // [x_delta] + __global double *t_cycleRadius // [x_delta] ){ unsigned int gid = get_global_id(0); Particle particle = t_particles[gid]; - + double cyclicRadius = t_cycleRadius[0]; + // 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); - + // abs(x) > Boundary and x < -Boundary + // abs(x) > Boundary and x > Boundary + 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); // Update result t_particles[gid] = particle; } From 029c0f3dce445ac268c41153140949fa65ac0f75 Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Fri, 9 Jun 2023 16:30:25 +0300 Subject: [PATCH 2/8] Refactor. Checkpoint. --- bubbleSim/collision.cpp | 105 +++++++--- bubbleSim/collision.h | 15 +- bubbleSim/datastreamer.cpp | 58 +++++- bubbleSim/datastreamer.h | 11 +- bubbleSim/opencl_kernels.cpp | 12 +- bubbleSim/particle.cpp | 18 +- bubbleSim/particle.h | 22 +- bubbleSim/simulation.cpp | 80 +++++-- bubbleSim/source.cpp | 44 ++-- configs/test_oclgrind.json | 2 +- kernels/kernel.cl | 389 +++++++++++++++++------------------ 11 files changed, 459 insertions(+), 297 deletions(-) diff --git a/bubbleSim/collision.cpp b/bubbleSim/collision.cpp index fc948d9..5d42812 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; @@ -70,9 +71,9 @@ void CollisionCellCollection::recalculate_cells( */ 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; } @@ -88,10 +89,10 @@ void CollisionCellCollection::recalculate_cells( 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,21 +105,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 = -M_PI / 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); } } +} - /*int particle1_idx = t_particles[0].idxCollisionCell; - std::cout << t_particles[0].E << ", " << t_particles[0].p_x << ", " - << t_particles[0].p_y << ", " << t_particles[0].p_z << std::endl; - std::cout << m_collisionCells[particle1_idx].theta << std::endl;*/ +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 c9dac17..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,7 +276,7 @@ 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.getParticles()[i].x * particleCollection.getParticles()[i].p_x + particleCollection.getParticles()[i].y * particleCollection.getParticles()[i].p_y + particleCollection.getParticles()[i].z * particleCollection.getParticles()[i].p_z) / (particleCollection.getParticles()[i].E * particleRadius); + 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)] * @@ -357,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..4021c61 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 = true; 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 1c08812..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"; @@ -29,9 +29,9 @@ OpenCLLoader::OpenCLLoader(std::string kernelPath, std::string kernelName) { std::string particleBubbleBoundaryStepKernelName = "particle_bubble_step_cyclic"; // NB: Collision: change - std::string cellAssignKernelName = "assign_cell_index_to_particle"; - std::string transformKernelName = "transform_momentum_massive"; - std::string particleBounceKernelName = "particle_bounce"; + 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); diff --git a/bubbleSim/particle.cpp b/bubbleSim/particle.cpp index a0e128d..1835cf4 100644 --- a/bubbleSim/particle.cpp +++ b/bubbleSim/particle.cpp @@ -393,23 +393,23 @@ 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)); } 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 c63eab7..39ade76 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -216,29 +216,65 @@ void Simulation::step(ParticleCollection& particles, cl::CommandQueue& cl_queue) { m_time += m_dt; // Move particles - - // Generate shift vector - if (i % 1 == 0) { - particles.readParticlesBuffer(cl_queue); - cells.generateShiftVector(generator_collision); - cells.writeShiftVectorBuffer(cl_queue); - // Assign particles to collision cells - cl_queue.enqueueNDRangeKernel( - t_cellAssignmentKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); - // 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())); - } cl_queue.enqueueNDRangeKernel(t_particleStepKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); cl_queue.enqueueNDRangeKernel(t_particleBounceKernel, cl::NullRange, - cl::NDRange(particles.getParticleCountTotal())); + 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); + + 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 pLocArray; + 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; @@ -328,23 +341,14 @@ int main(int argc, char* argv[]) { kernels.m_particleBounceKernel, kernels.getCommandQueue()); - particles.readParticlesBuffer(kernels.getCommandQueue()); - /*std::cout << particles.getParticles()[0].E << ", " - << particles.getParticles()[0].p_x << ", " - << particles.getParticles()[0].p_y << ", " - << particles.getParticles()[0].p_z << std::endl;*/ - /*px = 0; - py = 0; - pz = 0; - for (Particle& p : particles.getParticles()) { - px += p.p_x; - py += p.p_y; - pz += p.p_z; + /*for (Particle& p : particles.getParticles()) { + pLocArray[p.idxCollisionCell] += 1; } - std::cout << px / 750000. << std::endl; - std::cout << py / 750000. << std::endl; - std::cout << pz / 750000. << std::endl;*/ - + 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/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 767abe3..5be0d37 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,36 +769,36 @@ __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 *t_cycleRadius, + __global Particle *particles, + __constant double *boundaryRadius, __constant double *t_dt ){ unsigned int gid = get_global_id(0); - Particle particle = t_particles[gid]; + Particle particle = particles[gid]; // double Vx = particle.pX/particle.E; // double Vy = particle.pY/particle.E; // double Vz = particle.pZ/particle.E; @@ -804,15 +807,15 @@ __kernel void particle_step( 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 @@ -820,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]; @@ -900,115 +903,93 @@ __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; + + // Read variables to local registers + Particle particle = particles[gid]; + + if (gid == 0){ + printf("%p \n", &particle); } - else { - t_particles[gid].idxCollisionCell = 1 + x_index + y_index * maxCellIndex[0] + z_index * maxCellIndex[0] * maxCellIndex[0]; + if (gid == 100){ + printf("%p \n", &particle.x); } - - //t_particles[gid].idxCollisionCell = x_index + y_index + z_index; - -} + unsigned int r_maximumCellIndex = maxCellIndex[0]; + double r_cellLength = cellLength[0]; -__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); + // 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); - - // 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]; -} - -__kernel void transform_momentum( - __global Particle *t_particles, - __global CollisionCell *t_cells, +__kernel void collide_particles( + __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,14 +999,15 @@ __kernel void transform_momentum( 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 @@ -1052,7 +1034,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); @@ -1067,7 +1051,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 @@ -1097,52 +1083,48 @@ __kernel void transform_momentum( // printf("E: %.8f, pX: %.8f, pY: %.8f, pZ: %.8f\n", p0, p1, p2, p3); // } - t_particles[gid] = particle; + particles[gid] = particle; } } } - -__kernel void transform_momentum_massive( - __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); - // Copy object to register memory. Improves performance. - 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]; - //if (gid==0){ - // printf("Cell count: %i ,",cell.particle_count); - // } + CollisionCell cell = cells[particle.idxCollisionCell]; + 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; - double p0result, p1result, p2result, p3result; - double new_gamma; - double mgamma; if ((cell.particle_count > 1) && (cell.mass != 0)){ // Lorentz boost //if (gid==0){ - // printf("x: %.4f, %.4f, %.4f, cell: %i\n", particle.x, particle.y, particle.z, particle.idxCollisionCell); - //} - //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.E + - particle.pX * cell.vX + - particle.pY * cell.vY - particle.pZ * cell.vZ ); p1 = - particle.E * cell.gamma * cell.vX @@ -1155,7 +1137,7 @@ __kernel void transform_momentum_massive( + 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 + 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; @@ -1165,7 +1147,9 @@ __kernel void transform_momentum_massive( 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); @@ -1180,7 +1164,9 @@ __kernel void transform_momentum_massive( 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 @@ -1206,40 +1192,43 @@ __kernel void transform_momentum_massive( 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 *t_cycleRadius // [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]; - double cyclicRadius = t_cycleRadius[0]; + 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 + // If Particle is inside the boundary leave value same. Otherwise change the sign - // abs(x) < Boundary -> leave momentum - // abs(x) > Boundary and x < -Boundary - // abs(x) > Boundary and x > Boundary - 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); + // 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; } From 1f6871eb269fbc635549f87981d74c42887fff6d Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Fri, 9 Jun 2023 16:34:00 +0300 Subject: [PATCH 3/8] Test fix. --- bubbleSim/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bubbleSim/source.cpp b/bubbleSim/source.cpp index 3b3d075..8b79acd 100644 --- a/bubbleSim/source.cpp +++ b/bubbleSim/source.cpp @@ -82,7 +82,7 @@ int main(int argc, char* argv[]) { std::string s_configPath = argv[1]; // "config.json" std::string s_kernelPath = argv[2]; // "kernel.cl"; - s_configPath = "D:\\dev\\bubbleSim\\configs\\test_collision.json"; + // s_configPath = "D:\\dev\\bubbleSim\\configs\\test_collision.json"; std::cout << "Config path: " << s_configPath << std::endl; std::cout << "Kernel path: " << s_kernelPath << std::endl; From 67fa89a3c76d12c15925e8ac4d8b989776cb4803 Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Wed, 14 Jun 2023 15:07:19 +0300 Subject: [PATCH 4/8] Collision fix. --- bubbleSim/collision.cpp | 20 ++++++++++++++-- bubbleSim/simulation.cpp | 17 +++++++------- bubbleSim/source.cpp | 49 ++++++++++++++++++++++------------------ kernels/kernel.cl | 19 +++++----------- 4 files changed, 59 insertions(+), 46 deletions(-) diff --git a/bubbleSim/collision.cpp b/bubbleSim/collision.cpp index 5d42812..8e1d64c 100644 --- a/bubbleSim/collision.cpp +++ b/bubbleSim/collision.cpp @@ -64,8 +64,7 @@ void CollisionCellCollection::recalculate_cells( 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 */ @@ -76,6 +75,7 @@ void CollisionCellCollection::recalculate_cells( 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 @@ -83,6 +83,22 @@ 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/3, frames[i][4]) / (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]; diff --git a/bubbleSim/simulation.cpp b/bubbleSim/simulation.cpp index 39ade76..b45a6a6 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -221,7 +221,7 @@ void Simulation::step(ParticleCollection& particles, cl_queue.enqueueNDRangeKernel(t_particleBounceKernel, cl::NullRange, cl::NDRange(particles.getParticleCountTotal())); - std::array particleIdx = { 84522, 355324 , 429220, 538040, 704867 }; + // std::array particleIdx = { 84522, 355324 , 429220, 538040, 704867 }; // Generate shift vector if (i % 1 == 0) { @@ -240,7 +240,7 @@ void Simulation::step(ParticleCollection& particles, cells.recalculate_cells(particles.getParticles(), generator_collision); - int cell_idx = 0; + /*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; @@ -252,11 +252,12 @@ void Simulation::step(ParticleCollection& particles, 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, - config.particleTemperatureFalse); + 3*config.particleTemperatureFalse); } else { particleGenerator1 = ParticleGenerator(config.particleMassFalse, config.particleTemperatureFalse, - 30 * temperatureFalse, 1e-5 * temperatureFalse); + 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()); @@ -298,19 +299,22 @@ int main(int argc, char* argv[]) { initialNumberDensityFalse, initialEnergyDensityFalse); - - std::ofstream pLocStream("cellIdx.csv"); + bool streamCellCount = false; + std::ofstream pLocStream; std::vector pLocArray; - pLocArray.resize(config.collisionCellCount* config.collisionCellCount* config.collisionCellCount + 1, 0); + 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 }; + /*std::array particleIdx = { 84522, 355324 , 429220, 538040, 704867 }; for (int i : particleIdx) { particles.printParticleInfo(i); } - +*/ numType simTimeSinceLastStream = 0.; int stepsSinceLastStream = 0; @@ -336,19 +340,20 @@ int main(int argc, char* argv[]) { break; } if (b_collisionDevelopment) { - simulation.step(particles, cells, rn_generator, i, *stepKernel, - kernels.m_cellAssignmentKernel, kernels.m_rotationKernel, - kernels.m_particleBounceKernel, - kernels.getCommandQueue()); - - /*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);*/ + 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/kernels/kernel.cl b/kernels/kernel.cl index 5be0d37..b167474 100644 --- a/kernels/kernel.cl +++ b/kernels/kernel.cl @@ -926,12 +926,6 @@ __kernel void assign_particle_to_collision_cell( // Read variables to local registers Particle particle = particles[gid]; - if (gid == 0){ - printf("%p \n", &particle); - } - if (gid == 100){ - printf("%p \n", &particle.x); - } unsigned int r_maximumCellIndex = maxCellIndex[0]; double r_cellLength = cellLength[0]; @@ -990,14 +984,13 @@ __kernel void collide_particles( // If in bubble then cell number is doubled and second half is in bubble cells. if (particle.idxCollisionCell != 0){ CollisionCell cell = cells[particle.idxCollisionCell]; - - 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; - 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", From 2bbacf2a00d321525c8274b2b64d708550288201 Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Tue, 27 Jun 2023 14:04:15 +0300 Subject: [PATCH 5/8] Fix max_step parameter usage. --- bubbleSim/collision.cpp | 2 +- bubbleSim/simulation.cpp | 2 +- configs/test_oclgrind.json | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/bubbleSim/collision.cpp b/bubbleSim/collision.cpp index 8e1d64c..a6814a6 100644 --- a/bubbleSim/collision.cpp +++ b/bubbleSim/collision.cpp @@ -93,7 +93,7 @@ void CollisionCellCollection::recalculate_cells( } */ - if (t_rng.generate_number() >= std::pow(0.01/3, frames[i][4]) / (frames[i][5])) { + if (t_rng.generate_number() >= std::pow(0.01, frames[i][4]) / (9*frames[i][5])) { m_collisionCells[i].particle_count = (int)0; continue; } diff --git a/bubbleSim/simulation.cpp b/bubbleSim/simulation.cpp index b45a6a6..972fdd2 100644 --- a/bubbleSim/simulation.cpp +++ b/bubbleSim/simulation.cpp @@ -274,6 +274,6 @@ void Simulation::step(ParticleCollection& particles, particles.printParticleInfo(i); } }*/ - } + m_step += 1; } diff --git a/configs/test_oclgrind.json b/configs/test_oclgrind.json index 8b7167c..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 }, From a12b5621565ff5701b76461bee1e55e16ba637dd Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Wed, 28 Jun 2023 10:38:00 +0300 Subject: [PATCH 6/8] Bug fix --- bubbleSim/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bubbleSim/source.cpp b/bubbleSim/source.cpp index 224ea80..ac5fa73 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 = true; + bool b_collisionDevelopment = false; // Read configs and kernel file std::string s_configPath = argv[1]; // "config.json" From 45c0ab17e42fa3872535c6a449b72f1621950825 Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Wed, 28 Jun 2023 10:42:11 +0300 Subject: [PATCH 7/8] Bug fix --- bubbleSim/source.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/bubbleSim/source.cpp b/bubbleSim/source.cpp index ac5fa73..f31d8d4 100644 --- a/bubbleSim/source.cpp +++ b/bubbleSim/source.cpp @@ -129,6 +129,7 @@ int main(int argc, char* argv[]) { 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, From 2c21862fe855d7c7742f867a722aabbfe4601413 Mon Sep 17 00:00:00 2001 From: Martin Vasar Date: Wed, 28 Jun 2023 16:52:57 +0300 Subject: [PATCH 8/8] Memeory fix --- bubbleSim/datastreamer.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bubbleSim/datastreamer.h b/bubbleSim/datastreamer.h index 4021c61..872cf79 100644 --- a/bubbleSim/datastreamer.h +++ b/bubbleSim/datastreamer.h @@ -102,7 +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 = true; + bool m_initialized_Momentum = false; bool m_initialized_MomentumIn = false; bool m_initialized_MomentumOut = false; // Energy and number desnity profiles