Skip to content

Instantly share code, notes, and snippets.

@dextorious
Created May 20, 2017
Embed
What would you like to do?
#include <chrono>
#include <iostream>
#include <memory>
#include <random>
static void collision_kernel(double* n, double* ux, double* uy, double* rho,
const double* ex, const double* ey, const double* w,
const int Q, const int NX, const int NY, const double omega) {
for (int i = 0; i < NX; ++i)
for (int j = 0; j < NY; ++j) {
const int idx = i*NY + j;
ux[idx] = 0.0;
uy[idx] = 0.0;
rho[idx] = 0.0;
for (int q = 0; q < Q; ++q) {
const double nq = n[idx * Q + q];
rho[idx] += nq;
ux[idx] += ex[q] * nq;
uy[idx] += ey[q] * nq;
}
const double rho_inv = 1. / rho[idx];
ux[idx] *= rho_inv;
uy[idx] *= rho_inv;
const double usqr = ux[idx] * ux[idx] + uy[idx] * uy[idx];
for (int q = 0; q < Q; ++q) {
const double eu = 3 * (ex[q] * ux[idx] + ey[q] * uy[idx]);
const double neq = rho[idx] * w[q] * (1.0 + eu + 0.5*eu*eu - 1.5*usqr);
n[idx*Q + q] = (1 - omega)*n[idx*Q + q] + omega*neq;
}
}
}
int main() {
const unsigned int NX = 1000;
const unsigned int NY = 1000;
const unsigned int N = NX * NY;
const unsigned int Q = 9;
const double ex[9] = { 0., 1., 0., -1., 0., 1., -1., -1., 1. };
const double ey[9] = { 0., 0., 1., 0., -1., 1., 1., -1., -1. };
const double w[9] = { 4. / 9., 1. / 9., 1. / 9., 1. / 9., 1. / 9., 1. / 36., 1. / 36., 1. / 36., 1. / 36. };
double* n = new double[N*Q];
double* ux = new double[N];
double* uy = new double[N];
double* rho = new double[N];
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<double> dist(0.0, 1.0);
for (auto i = 0; i < N; ++i) {
ux[i] = 0.0;
uy[i] = 0.0;
rho[i] = 1.0;
for (auto q = 0; q < 9; ++q) {
n[9 * i + q] = dist(mt);
}
}
const unsigned int benchmarkCount = 100;
auto start = std::chrono::steady_clock::now();
for (unsigned int t = 0; t < benchmarkCount; ++t)
collision_kernel(n, ux, uy, rho, ex, ey, w, Q, NX, NY, 0.6);
auto end = std::chrono::steady_clock::now();
auto diff = end - start;
std::cout << "n[0] = " << n[0] << std::endl;
std::cout << "Time = " << std::chrono::duration<double, std::milli>(diff).count() / benchmarkCount << " ms" << std::endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment