diff --git a/PixelEngine/Fluid.cs b/PixelEngine/Fluid.cs index 5a72700..4a1c0d5 100644 --- a/PixelEngine/Fluid.cs +++ b/PixelEngine/Fluid.cs @@ -1,4 +1,5 @@ using System.Numerics; +using Hybridizer.Runtime.CUDAImports; using Raylib_CsLo; using static Raylib_CsLo.RayMath; @@ -43,6 +44,7 @@ public class Fluid this.Vy0 = new float[N*N]; } + private static dynamic wrapper; public void step() { int N = this.size; float visc = this.visc; @@ -59,14 +61,15 @@ public class Fluid diffuse(2, Vy0, Vy, visc, dt); project(Vx0, Vy0, Vx, Vy); - - advect(1, Vx, Vx0, Vx0, Vy0, dt); - advect(2, Vy, Vy0, Vx0, Vy0, dt); + HybRunner runner = HybRunner.Cuda().SetDistrib(32, 32, 16, 16, 1, 0); + wrapper = runner.Wrap(new Game()); + wrapper.Run(wrapper.advect(this,1, Vx, Vx0, Vx0, Vy0, dt)); + advect(this,2, Vy, Vy0, Vx0, Vy0, dt); project(Vx, Vy, Vx0, Vy0); diffuse(0, s, density, diff, dt); - advect(0, density, s, Vx, Vy, dt); + advect(this,0, density, s, Vx, Vy, dt); } @@ -94,8 +97,8 @@ public class Fluid public int IX(int x, int y) { - var X = Normalize(x, 0, size - 1);; - var Y = Normalize(y, 0, size - 1); + var X = Clamp(x, 0, size - 1);; + var Y = Clamp(y, 0, size - 1); return (int) (X + Y * size); } @@ -105,23 +108,23 @@ public class Fluid } public void reduceDensity(int x, int y, float amount, int rSize) { + + Parallel.For(-rSize, rSize, i => + { + for (int j = -rSize; j < rSize; j++) + { + int index = IX(x + i, y + j); - for (int i = -rSize; i < rSize; i++) - { - for (int j = -rSize; j < rSize; j++) - { - int index = IX(x + i, y + j); - - if (this.density[index] - amount < 0) - { - this.density[index] = 0; - } - else - { - this.density[index] -= amount; - } - } - } + if (this.density[index] - amount < 0) + { + this.density[index] = 0; + } + else + { + this.density[index] -= amount; + } + } + }); } @@ -138,67 +141,77 @@ public class Fluid } void lin_solve(int b, float[] x, float[] x0, float a, float c) { float cRecip = 1.0f / c; - for (int k = 0; k < iter; k++) { + Parallel.For(0, iter, k => + { for (int j = 1; j < size - 1; j++) { for (int i = 1; i < size - 1; i++) { x[IX(i, j)] = (x0[IX(i, j)] - + a*( x[IX(i+1, j)] - +x[IX(i-1, j)] - +x[IX(i, j+1)] - +x[IX(i, j-1)] - )) * cRecip; + + a*( x[IX(i+1, j)] + +x[IX(i-1, j)] + +x[IX(i, j+1)] + +x[IX(i, j-1)] + )) * cRecip; } } set_bnd(b, x); - } + }); } - void project(float[] velocX, float[] velocY, float[] p, float[] div) { - for (int j = 1; j < size - 1; j++) { - for (int i = 1; i < size - 1; i++) { - div[IX(i, j)] = -0.5f*( - velocX[IX(i+1, j)] - -velocX[IX(i-1, j)] - +velocY[IX(i, j+1)] - -velocY[IX(i, j-1)] - )/size; + void project(float[] velocX, float[] velocY, float[] p, float[] div) + { + Parallel.For(1, size, j => + { + for (int i = 1; i < size - 1; i++) + { + div[IX(i, j)] = -0.5f * ( + velocX[IX(i + 1, j)] + - velocX[IX(i - 1, j)] + + velocY[IX(i, j + 1)] + - velocY[IX(i, j - 1)] + ) / size; p[IX(i, j)] = 0; } - } + }); set_bnd(0, div); set_bnd(0, p); lin_solve(0, p, div, 1, 4); - for (int j = 1; j < size - 1; j++) { - for (int i = 1; i < size - 1; i++) { - velocX[IX(i, j)] -= 0.5f * ( p[IX(i+1, j)] - -p[IX(i-1, j)]) * size; - velocY[IX(i, j)] -= 0.5f * ( p[IX(i, j+1)] - -p[IX(i, j-1)]) * size; + Parallel.For(1, size, j => + { + for (int i = 1; i < size - 1; i++) + { + velocX[IX(i, j)] -= 0.5f * (p[IX(i + 1, j)] + - p[IX(i - 1, j)]) * size; + velocY[IX(i, j)] -= 0.5f * (p[IX(i, j + 1)] + - p[IX(i, j - 1)]) * size; } - } + }); set_bnd(1, velocX); set_bnd(2, velocY); } - void advect(int b, float[] d, float[] d0, float[] velocX, float[] velocY, float dt) { + + + + [EntryPoint("run")] + public void advect(Fluid fluid,int b, float[] d, float[] d0, float[] velocX, float[] velocY, float dt) { float i0, i1, j0, j1; - float dtx = dt * (size - 2); - float dty = dt * (size - 2); + float dtx = dt * (fluid.size - 2); + float dty = dt * (fluid.size - 2); float s0, s1, t0, t1; float tmp1, tmp2, x, y; - float Nfloat = size; + float Nfloat = fluid.size; float ifloat, jfloat; int i, j; - for (j = 1, jfloat = 1; j < size - 1; j++, jfloat++) { - for (i = 1, ifloat = 1; i < size - 1; i++, ifloat++) { - tmp1 = dtx * velocX[IX(i, j)]; - tmp2 = dty * velocY[IX(i, j)]; + for (j = 1, jfloat = 1; j < fluid.size - 1; j++, jfloat++) { + for (i = 1, ifloat = 1; i < fluid.size - 1; i++, ifloat++) { + tmp1 = dtx * velocX[fluid.IX(i, j)]; + tmp2 = dty * velocY[fluid.IX(i, j)]; x = ifloat - tmp1; y = jfloat - tmp2; @@ -222,23 +235,27 @@ public class Fluid int j1i = (int)j1; // DOUBLE CHECK THIS!!! - d[IX(i, j)] = - s0 * (t0 * d0[IX(i0i, j0i)] + t1 * d0[IX(i0i, j1i)]) + - s1 * (t0 * d0[IX(i1i, j0i)] + t1 * d0[IX(i1i, j1i)]); + d[fluid.IX(i, j)] = + s0 * (t0 * d0[fluid.IX(i0i, j0i)] + t1 * d0[fluid.IX(i0i, j1i)]) + + s1 * (t0 * d0[fluid.IX(i1i, j0i)] + t1 * d0[fluid.IX(i1i, j1i)]); } } - set_bnd(b, d); + fluid.set_bnd(b, d); } - void set_bnd(int b, float[] x) { - for (int i = 1; i < size - 1; i++) { + void set_bnd(int b, float[] x) + { + Parallel.For(1, size, i => + { x[IX(i, 0 )] = b == 2 ? -x[IX(i, 1 )] : x[IX(i, 1 )]; x[IX(i, size-1)] = b == 2 ? -x[IX(i, size-2)] : x[IX(i, size-2)]; - } - for (int j = 1; j < size - 1; j++) { + }); + + Parallel.For(1, size, j => + { x[IX(0, j)] = b == 1 ? -x[IX(1, j)] : x[IX(1, j)]; - x[IX(size-1, j)] = b == 1 ? -x[IX(size-2, j)] : x[IX(size-2, j)]; - } + x[IX(size - 1, j)] = b == 1 ? -x[IX(size - 2, j)] : x[IX(size - 2, j)]; + }); x[IX(0, 0)] = 0.5f * (x[IX(1, 0)] + x[IX(0, 1)]); x[IX(0, size-1)] = 0.5f * (x[IX(1, size-1)] + x[IX(0, size-2)]); diff --git a/PixelEngine/PixelEngine.csproj b/PixelEngine/PixelEngine.csproj index d17f885..7eba0f1 100644 --- a/PixelEngine/PixelEngine.csproj +++ b/PixelEngine/PixelEngine.csproj @@ -9,6 +9,7 @@ +