collide_bench.cpp

#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;
}