|
|
|
@ -5,131 +5,17 @@ |
|
|
|
|
#include <omp.h> |
|
|
|
|
#include <iostream> |
|
|
|
|
#include "FPS.h" |
|
|
|
|
#include "Pendulum.h" |
|
|
|
|
#include <QOpenGLContext> |
|
|
|
|
#include <QOpenGLShaderProgram> |
|
|
|
|
#include <QOffscreenSurface> |
|
|
|
|
|
|
|
|
|
Simulation::Simulation() { |
|
|
|
|
ups = new FPS; |
|
|
|
|
ups->setUpdateInterval(100); |
|
|
|
|
|
|
|
|
|
timer = new QTimer(this); |
|
|
|
|
QTimer::connect(timer, &QTimer::timeout, this, &Simulation::update); |
|
|
|
|
timer->setInterval(updateInterval); |
|
|
|
|
timer->start(); |
|
|
|
|
}; |
|
|
|
|
|
|
|
|
|
void Simulation::initialize() { |
|
|
|
|
|
|
|
|
|
// OpenGL Offscreen functions
|
|
|
|
|
{ |
|
|
|
|
context = new QOpenGLContext; |
|
|
|
|
auto format = QSurfaceFormat::defaultFormat(); |
|
|
|
|
format.setVersion(4, 3); |
|
|
|
|
context->setFormat(format); |
|
|
|
|
context->create(); |
|
|
|
|
context->makeCurrent(offSurface); |
|
|
|
|
f = new QOpenGLFunctions_4_3_Core; |
|
|
|
|
if (!f->initializeOpenGLFunctions()){ |
|
|
|
|
std::string message = "OpenGL 4.3 is required for Compute Shaders! Falling back to CPU-Multithreading."; |
|
|
|
|
printf("%s\n", message.c_str()); |
|
|
|
|
useGPUAcceleration = false; |
|
|
|
|
emit gpuNotSupported(message); |
|
|
|
|
} else if (!context->hasExtension("GL_ARB_gpu_shader_fp64")){ |
|
|
|
|
std::string message = "Double precision not supported by OpenGL! Falling back to CPU-Multithreading."; |
|
|
|
|
printf("%s\n", message.c_str()); |
|
|
|
|
useGPUAcceleration = false; |
|
|
|
|
emit gpuNotSupported(message); |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
if (!useGPUAcceleration){ |
|
|
|
|
context->doneCurrent(); |
|
|
|
|
return; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// Read GPU Limits
|
|
|
|
|
{ |
|
|
|
|
for (int i = 0; i < 3; i++){ |
|
|
|
|
f->glGetIntegeri_v(GL_MAX_COMPUTE_WORK_GROUP_COUNT, i, gpuLimits.maxWGCount + i); |
|
|
|
|
f->glGetIntegeri_v(GL_MAX_COMPUTE_WORK_GROUP_SIZE, i, gpuLimits.maxWGSize + i); |
|
|
|
|
} |
|
|
|
|
f->glGetIntegerv(GL_MAX_COMPUTE_WORK_GROUP_INVOCATIONS, &gpuLimits.maxWGInvocations); |
|
|
|
|
|
|
|
|
|
printf("Max work group count: (%d, %d, %d)\n", gpuLimits.maxWGCount[0], gpuLimits.maxWGCount[1], gpuLimits.maxWGCount[2]); |
|
|
|
|
printf("Max work group size: (%d, %d, %d)\n", gpuLimits.maxWGSize[0], gpuLimits.maxWGSize[1], gpuLimits.maxWGSize[2]); |
|
|
|
|
printf("Max work group invocations (x * y * z): %d\n", gpuLimits.maxWGInvocations); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
// Shader
|
|
|
|
|
{ |
|
|
|
|
program = new QOpenGLShaderProgram; |
|
|
|
|
program->addShaderFromSourceFile(QOpenGLShader::Compute, ":/shaders/compute.glsl"); |
|
|
|
|
program->link(); |
|
|
|
|
|
|
|
|
|
f->glGenBuffers(6, &SSBO.positions); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, SSBO.positions); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, SSBO.velocities); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, SSBO.invMasses); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, SSBO.lengths); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 4, SSBO.indices); |
|
|
|
|
f->glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 5, SSBO.segmentCounts); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
context->doneCurrent(); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
// When the layout changes
|
|
|
|
|
void Simulation::updateGPUData() { |
|
|
|
|
context->makeCurrent(offSurface); |
|
|
|
|
|
|
|
|
|
std::vector<GLdouble> positions; |
|
|
|
|
std::vector<GLdouble> velocities; |
|
|
|
|
std::vector<GLdouble> invMasses; |
|
|
|
|
std::vector<GLdouble> lengths; |
|
|
|
|
std::vector<GLuint> indices; |
|
|
|
|
std::vector<GLuint> segmentCounts; |
|
|
|
|
|
|
|
|
|
GLuint currentIndex = 0; |
|
|
|
|
for (const Pendulum * p : pendula){ |
|
|
|
|
indices.push_back(currentIndex); |
|
|
|
|
currentIndex += p->X.size(); |
|
|
|
|
segmentCounts.push_back(p->X.size()); |
|
|
|
|
for (size_t i = 0; i < p->X.size(); i++){ |
|
|
|
|
positions.push_back(p->X[i].x); |
|
|
|
|
positions.push_back(p->X[i].y); |
|
|
|
|
velocities.push_back(p->V[i].x); |
|
|
|
|
velocities.push_back(p->V[i].y); |
|
|
|
|
invMasses.push_back(1.0 / p->M[i]); |
|
|
|
|
lengths.push_back(p->L[i]); |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
auto posSize = GLsizeiptr(positions.size() * sizeof(GLdouble)); |
|
|
|
|
auto velSize = GLsizeiptr(velocities.size() * sizeof(GLdouble)); |
|
|
|
|
auto invMassSize = GLsizeiptr(invMasses.size() * sizeof(GLdouble)); |
|
|
|
|
auto lengthSize = GLsizeiptr(lengths.size() * sizeof(GLdouble)); |
|
|
|
|
auto indicesSize = GLsizeiptr(indices.size() * sizeof(GLuint)); |
|
|
|
|
auto segmentCountsSize = GLsizeiptr(segmentCounts.size() * sizeof(GLuint)); |
|
|
|
|
|
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.positions); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, posSize, positions.data(), GL_DYNAMIC_DRAW); |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.velocities); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, velSize, velocities.data(), GL_DYNAMIC_DRAW); |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.invMasses); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, invMassSize, invMasses.data(), GL_STATIC_DRAW); |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.lengths); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, lengthSize, lengths.data(), GL_STATIC_DRAW); |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.indices); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, indicesSize, indices.data(), GL_STATIC_DRAW); |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.segmentCounts); |
|
|
|
|
f->glBufferData(GL_SHADER_STORAGE_BUFFER, segmentCountsSize, segmentCounts.data(), GL_STATIC_DRAW); |
|
|
|
|
|
|
|
|
|
context->doneCurrent(); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
// Every few milliseconds
|
|
|
|
|
void Simulation::update() { |
|
|
|
|
if (!isPlaying) |
|
|
|
|
return; |
|
|
|
@ -143,67 +29,16 @@ void Simulation::update() { |
|
|
|
|
double newPotentialEnergy = 0; |
|
|
|
|
double newKineticEnergy = 0; |
|
|
|
|
|
|
|
|
|
if (useGPUAcceleration) { |
|
|
|
|
context->makeCurrent(offSurface); |
|
|
|
|
program->bind(); |
|
|
|
|
f->glUniform1d(program->uniformLocation("h"), h); |
|
|
|
|
f->glUniform1d(program->uniformLocation("gravity"), gravity); |
|
|
|
|
f->glUniform1ui(program->uniformLocation("substeps"), substeps); |
|
|
|
|
|
|
|
|
|
f->glDispatchCompute(pendula.size(), 1, 1); |
|
|
|
|
f->glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT); |
|
|
|
|
|
|
|
|
|
// Read updated positions
|
|
|
|
|
{ |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.positions); |
|
|
|
|
auto * newPositions = (double*) f->glMapBuffer(GL_SHADER_STORAGE_BUFFER, GL_READ_ONLY); |
|
|
|
|
int index = 0; |
|
|
|
|
for (Pendulum * p : pendula){ |
|
|
|
|
for (Vector &point : p->X){ |
|
|
|
|
point.x = newPositions[index++]; |
|
|
|
|
point.y = newPositions[index++]; |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
f->glUnmapBuffer(GL_SHADER_STORAGE_BUFFER); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
// Read updated velocities
|
|
|
|
|
{ |
|
|
|
|
f->glBindBuffer(GL_SHADER_STORAGE_BUFFER, SSBO.velocities); |
|
|
|
|
auto * newVelocities = (double*) f->glMapBuffer(GL_SHADER_STORAGE_BUFFER, GL_READ_ONLY); |
|
|
|
|
int index = 0; |
|
|
|
|
for (Pendulum * p : pendula){ |
|
|
|
|
for (Vector &velocity : p->V){ |
|
|
|
|
velocity.x = newVelocities[index++]; |
|
|
|
|
velocity.y = newVelocities[index++]; |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
f->glUnmapBuffer(GL_SHADER_STORAGE_BUFFER); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
#pragma omp parallel for |
|
|
|
|
for (int i = 0; i < pendula.size(); i++){ // NOLINT(*-loop-convert) not ranged based for msvc
|
|
|
|
|
double localPotentialEnergy = pendula[i]->potentialEnergy(gravity); |
|
|
|
|
double localKineticEnergy = pendula[i]->kineticEnergy(); |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newPotentialEnergy += localPotentialEnergy; |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newKineticEnergy += localKineticEnergy; |
|
|
|
|
} |
|
|
|
|
context->doneCurrent(); |
|
|
|
|
} else { |
|
|
|
|
// CPU-Multithreading
|
|
|
|
|
#pragma omp parallel for |
|
|
|
|
for (int i = 0; i < pendula.size(); i++){ // NOLINT(*-loop-convert) not ranged based for msvc
|
|
|
|
|
for (int k = 0; k < substeps; k++) |
|
|
|
|
pendula[i]->update(h, gravity); |
|
|
|
|
double localPotentialEnergy = pendula[i]->potentialEnergy(gravity); |
|
|
|
|
double localKineticEnergy = pendula[i]->kineticEnergy(); |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newPotentialEnergy += localPotentialEnergy; |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newKineticEnergy += localKineticEnergy; |
|
|
|
|
} |
|
|
|
|
#pragma omp parallel for |
|
|
|
|
for (int i = 0; i < pendula.size(); i++){ // NOLINT(*-loop-convert) // not ranged based for msvc
|
|
|
|
|
for (int k = 0; k < substeps; k++) |
|
|
|
|
pendula[i]->update(h, gravity); |
|
|
|
|
double localPotentialEnergy = pendula[i]->potentialEnergy(gravity); |
|
|
|
|
double localKineticEnergy = pendula[i]->kineticEnergy(); |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newPotentialEnergy += localPotentialEnergy; |
|
|
|
|
#pragma omp atomic |
|
|
|
|
newKineticEnergy += localKineticEnergy; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
potentialEnergy = newPotentialEnergy; |
|
|
|
@ -224,8 +59,6 @@ void Simulation::clearPendula() { |
|
|
|
|
pendula.shrink_to_fit(); |
|
|
|
|
|
|
|
|
|
updateEnergy(); |
|
|
|
|
if (useGPUAcceleration) |
|
|
|
|
updateGPUData(); |
|
|
|
|
|
|
|
|
|
emit layoutChanged(); |
|
|
|
|
} |
|
|
|
@ -237,8 +70,6 @@ void Simulation::addPendula(const std::vector<Pendulum *> &add) { |
|
|
|
|
pendula.push_back(p); |
|
|
|
|
|
|
|
|
|
updateEnergy(); |
|
|
|
|
if (useGPUAcceleration) |
|
|
|
|
updateGPUData(); |
|
|
|
|
|
|
|
|
|
emit layoutChanged(); |
|
|
|
|
} |
|
|
|
@ -251,9 +82,3 @@ void Simulation::updateEnergy() { |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
void Simulation::useGPUChanged(int state) { |
|
|
|
|
useGPUAcceleration = state == Qt::Checked; |
|
|
|
|
if (useGPUAcceleration) |
|
|
|
|
updateGPUData(); |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|