Skip to content

Instantly share code, notes, and snippets.

@dextorious
Created May 20, 2017 01:59
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save dextorious/d987865a7da147645ae34cc17a87729d to your computer and use it in GitHub Desktop.
import std.datetime;
import std.stdio;
import std.random;
import mir.ndslice;
static void collide_kernel(Tensor, Matrix)
(Tensor n, Matrix ux, Matrix uy, Matrix rho,
immutable double[] ex, immutable double[] ey, immutable double[] w,
immutable uint Q, immutable uint NX, immutable uint NY, immutable double omega)
{
for (int i = 0; i < NX; ++i)
for (int j = 0; j < NY; ++j) {
uint idx = i*NY + j;
ux[i,j] = 0.0;
uy[i,j] = 0.0;
rho[i,j] = 0.0;
for (uint q = 0; q < Q; ++q) {
double nq = n[i,j,q];
rho[i,j] += nq;
ux[i,j] += ex[q] * nq;
uy[i,j] += ey[q] * nq;
}
double rho_inv = 1. / rho[i,j];
ux[i,j] *= rho_inv;
uy[i,j] *= rho_inv;
double usqr = ux[i,j] * ux[i,j] + uy[i,j] * uy[i,j];
for (int q = 0; q < Q; ++q) {
double eu = 3 * (ex[q] * ux[i,j] + ey[q] * uy[i,j]);
double neq = rho[i,j] * w[q] * (1.0 + eu + 0.5*eu*eu - 1.5*usqr);
n[i,j,q] = (1 - omega)*n[i,j,q] + omega*neq;
}
}
}
int main(string[] args) {
immutable uint NX = 1000;
immutable uint NY = 1000;
immutable uint Q = 9;
auto n = new double[NX*NY*Q].sliced(NX,NY,Q);
auto ux = new double[NX*NY].sliced(NX,NY);
auto uy = new double[NX*NY].sliced(NX,NY);
auto rho = new double[NX*NY].sliced(NX,NY);
immutable double[Q] ex = [0., 1., 0., -1., 0., 1., -1., -1., 1.];
immutable double[Q] ey = [0., 0., 1., 0., -1., 1., 1., -1., -1.];
immutable double[Q] w = [4. / 9.,
1. / 9., 1. / 9., 1. / 9., 1. / 9.,
1. / 36., 1. / 36., 1. / 36., 1. / 36.];
Random gen;
for (uint i = 0; i < NX; ++i)
for (uint j = 0; j < NY; ++j) {
ux[i,j] = 0;
uy[i,j] = 0;
rho[i,j] = 1;
for (uint q = 0; q < Q; ++q)
n[i,j,q] = uniform(0.0, 1.0, gen);
}
uint benchmarkCount = 100u;
StopWatch sw;
alias Tensor = typeof(n);
alias Matrix = typeof(ux);
sw.start();
for (uint t=0; t<benchmarkCount; ++t)
collide_kernel!(Tensor, Matrix)(n, ux, uy, rho, ex, ey, w, Q, NX, NY, 0.6);
sw.stop();
writeln(n[0,0,0]);
double t1 = sw.peek().usecs;
writeln(t1 / 1000 / benchmarkCount);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment