Select Git revision
Particle.cpp NaN GiB
/// @file Particle.cpp
/// @author Armin Co
///
#include <complex>
#include "Particle.hpp"
#include "Utils.hpp"
#include "Log.hpp"
#include "Worker.hpp"
#include <imgui.h>
#include <GuiCore/Timer.hpp>
using namespace QPong;
double QPong::posAt(int pos, int arraySize)
{
return 2.0 * pos / arraySize - 1.0;
}
double QPong::xAtIndex(int i, int arraySize)
{
return posAt(i % arraySize, arraySize);
}
double QPong::yAtIndex(int i, int arraySize)
{
return posAt(i / arraySize, arraySize);
}
double QPong::potentialAt(int pos, const int arraySize, const double hbar)
{
if (pos > arraySize / 2)
{
pos -= arraySize;
}
constexpr double pi{3.14159265358979323846264};
return pos * pi * hbar;
}
double QPong::pxAtIndex(int i, int arraySize, double hbar)
{
return potentialAt(i % arraySize, arraySize, hbar);
}
double QPong::pyAtIndex(int i, int arraySize, double hbar)
{
return potentialAt(i / arraySize, arraySize, hbar);
}
Particle::Particle(Properties properties, GameOptions &options, Player &playerOne, Player &playerTwo, Bat &leftBat, Bat &rightBat, float *staticBorders, float *dynamicPotential)
: m_properties{properties}
, m_options{options}
, m_playerOne{playerOne}
, m_playerTwo{playerTwo}
, m_leftBat{leftBat}
, m_rightBat{rightBat}
, m_staticBorders{staticBorders}
, m_dynamicPotential{dynamicPotential}
{
Log::get()->debug("Init Particle");
m_xAt = new double[m_properties.elements()];
m_yAt = new double[m_properties.elements()];
m_E_At = new double[m_properties.elements()];
for (int i{0}; i < m_properties.elements(); ++i)
{
m_xAt[i] = xAtIndex(i, m_properties.size);
m_yAt[i] = yAtIndex(i, m_properties.size);
m_E_At[i] = (pow2(pxAtIndex(i, m_properties.size, m_properties.hbar)) + pow2(pyAtIndex(i, m_properties.size, m_properties.hbar))) / (2.0 * 1.0);
}
int memorySize = sizeof(fftw_complex) * m_properties.elements();
m_psi = static_cast<std::complex<double> *>(fftw_malloc(memorySize));
m_momentum = static_cast<std::complex<double> *>(fftw_malloc(memorySize));
initialiseParticleMomentum();
if (options.game == true)
{
calculateBorderPotential(m_staticBorders, m_yAt, m_properties.elements());
}
else
{
// "Simulation"
calculateGaussBox(m_staticBorders, m_xAt, m_yAt, m_properties.elements());
}
fftw_init_threads();
fftw_plan_with_nthreads(m_properties.threads);
m_planForward = fftw_plan_dft_2d(m_properties.size, m_properties.size, (fftw_complex *)m_psi, (fftw_complex *)m_psi, FFTW_FORWARD, FFTW_MEASURE);
fftw_plan_with_nthreads(m_properties.threads);
m_planBackwards = fftw_plan_dft_2d(m_properties.size, m_properties.size, (fftw_complex *)m_psi, (fftw_complex *)m_psi, FFTW_BACKWARD, FFTW_MEASURE);
int numberOfThreads = m_properties.threads;
for (int i = 0; i < numberOfThreads; i++)
{
const int from = i * (m_properties.elements() / numberOfThreads);
int to = (i+1) * (m_properties.elements() / numberOfThreads);
if (i+1 == numberOfThreads)
{
to = m_properties.elements();
}
m_applyPotentialWorkers.emplace_back(std::make_unique<WorkerApplyPotential>(this, from, to));
m_moveForwardWorkers.emplace_back(std::make_unique<WorkerMoveForward>(this, from, to));
}
Log::get()->debug("Init Particle with {} elements done.", m_properties.elements());
}
Particle::~Particle()
{
Log::get()->debug("Deleting Particle");
for (auto &worker : m_applyPotentialWorkers)
{
worker->quit();
}
for (auto &worker : m_moveForwardWorkers)
{
worker->quit();
}
std::this_thread::sleep_for(std::chrono::milliseconds(500));
fftw_free(m_psi);
delete (m_xAt);
delete (m_yAt);
delete (m_E_At);
delete (m_momentum);
}
std::complex<double> *Particle::getPsi() const
{
return m_psi;
}
std::complex<double> *Particle::getMomentum() const
{
return m_momentum;
}
int Particle::arrayElementCount() const
{
return pow2(m_properties.size);
}
void Particle::initialiseParticleMomentum()
{
constexpr std::complex<double> ci{0.0, 1.0}; // complex number
const double hbar = m_properties.hbar;
const double A0 = m_properties.startValue;
const double px0 = m_properties.momentum.x / 20.0;
const double py0 = m_properties.momentum.y / 20.0;
const double x0 = m_properties.position.x;
const double y0 = m_properties.position.y;
const double lambda = m_properties.sharpness;
double sum{0.0};
double m_particleAbsoluteSumStart = 0.0;
double y{0.0};
double x{0.0};
for (auto i{0}; i < arrayElementCount(); ++i)
{
x = m_xAt[i];
y = m_yAt[i];
m_psi[i] = A0 * exp(ci / hbar * (x * px0 + y * py0) - lambda * (pow2(x - x0) + pow2(y - y0)));
sum += pow2(std::real(m_psi[i])) + pow2(std::imag(m_psi[i]));
}
for (auto i{0}; i < arrayElementCount(); ++i)
{
m_psi[i] *= 10.0 / sqrt(sum);
m_particleAbsoluteSumStart += pow2(std::real(m_psi[i])) + pow2(std::imag(m_psi[i]));
}
Log::get()->trace("Particle initialised");
}
double Particle::absorbParticle(int i, double c)
{
double normPrevious = pow2(std::real(m_psi[i])) + pow2(std::imag(m_psi[i]));
double absorb = pow2(std::cos(c * 10));
m_psi[i] *= absorb;
double delta = normPrevious - (pow2(std::real(m_psi[i])) + pow2(std::imag(m_psi[i])));
return delta;
}
void Particle::applyPotentials(int indexBegin, int indexEnd, double dt)
{
constexpr std::complex<double> ci{0.0, 1.0};
const double hbar = m_properties.hbar;
double x{0.0};
double y{0.0};
const double potential = 1.0;
double particleAbsorbtAtPlayerOne = 0.0;
double particleAbsorbtAtPlayerTwo = 0.0;
for (int i{indexBegin}; i < indexEnd; ++i)
{
double potSum = 1.0;
if (m_options.potentialsEnabled)
{
potSum = 0.0;
x = m_xAt[i];
y = m_yAt[i];
if (m_options.absorbtionEnabled)
{
if (x < m_leftBat.getX())
{
double distanceToEdge = -m_leftBat.getX() + x;
double c = distanceToEdge * (M_PI / 2.0) / (1.0 + m_leftBat.getX());
particleAbsorbtAtPlayerOne += absorbParticle(i, c);
}
if (x > m_rightBat.getX())
{
double distanceToEdge = x - m_rightBat.getX();
double c = distanceToEdge * (M_PI / 2.0) / (1.0 + m_rightBat.getX());
particleAbsorbtAtPlayerTwo += absorbParticle(i, c);
}
}
double potLeftBat = m_leftBat.getPotential(x, y);
double potRightBat = m_rightBat.getPotential(x, y);
potSum = potLeftBat + potRightBat + m_staticBorders[i];
}
auto psi = m_psi[i] * exp(-ci * dt / hbar * potential * potSum);
m_psi[i] = psi;
if (m_options.potentialsEnabled)
{
m_dynamicPotential[i] = potSum;
}
}
m_playerOne.addScore(particleAbsorbtAtPlayerTwo);
m_playerTwo.addScore(particleAbsorbtAtPlayerOne);
}
void Particle::moveForward(int indexBegin, int indexEnd, double dt)
{
constexpr std::complex<double> ci{0.0, 1.0};
const double hbar = m_properties.hbar;
const double N = 1.0 / static_cast<double>(arrayElementCount());
for (int i{indexBegin}; i < indexEnd; ++i)
{
m_psi[i] *= exp(-ci * dt / hbar * m_E_At[i]) * N;
}
}
void Particle::updateMomentumView(int indexBegin, int indexEnd)
{
GuiCore::Timer timer("updateMomentumView");
int x{0};
int y{0};
for (int i{indexBegin}; i < indexEnd; ++i)
{
x = i % m_properties.size;
y = i / m_properties.size;
auto halfArraySize = m_properties.size / 2;
if (y > halfArraySize)
{
y = y - halfArraySize;
}
else
{
y = y + halfArraySize;
}
if (x > halfArraySize)
{
x = x - halfArraySize;
}
else
{
x = x + halfArraySize;
}
m_momentum[i] = m_psi[(y * m_properties.size) - x] * 100.0;
}
}
void Particle::runWorkers(std::vector<std::unique_ptr<WorkerParticle>> &workers, double dt)
{
for (auto &worker : workers)
{
worker->start(dt);
}
bool working {true};
while (working)
{
working = false;
for (auto &worker : workers)
{
if (worker->isWorking())
{
std::this_thread::yield();
working = true;
break;
}
}
}
}
void Particle::propagateStep(double dt)
{
runWorkers(m_applyPotentialWorkers, dt);
GuiCore::Timer timerFFT {"",false};
m_timeFFT = 0.0;
fftw_execute(m_planForward);
m_timeFFT = timerFFT.stop();
runWorkers(m_moveForwardWorkers, dt);
if (m_options.momentumEnabled)
{
GuiCore::Timer timer{"momentum view update"};
updateMomentumView(0, m_properties.elements());
}
timerFFT.start();
fftw_execute(m_planBackwards);
m_timeFFT += timerFFT.stop();
}
void Particle::update(double dt)
{
propagateStep(dt);
}
double Particle::getTimeFFT() const
{
return m_timeFFT;
}