elias/FluidSimulation
1
0
Fork 0
mirror of https://github.com/eliasstepanik/FluidSimulation.git synced 2026-01-19 09:28:29 +00:00
Code Issues Packages Projects Releases Wiki Activity

asd

Browse Source
This commit is contained in:
Elias Stepanik 2022-10-28 01:07:57 +02:00
parent 6ccbbdfba3
commit 653e85ce83
2 changed files with 82 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

143
PixelEngine/Fluid.cs
View File

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

1
PixelEngine/PixelEngine.csproj
View File

@ -9,6 +9,7 @@
</PropertyGroup> </PropertyGroup>
<ItemGroup> <ItemGroup>
<PackageReference Include="Altimesh.Hybridizer.Runtime.CUDAImports" Version="2.0.0" />
<PackageReference Include="Raylib-CsLo" Version="4.2.0.3" /> <PackageReference Include="Raylib-CsLo" Version="4.2.0.3" />
<PackageReference Include="System.Drawing.Common" Version="6.0.0" /> <PackageReference Include="System.Drawing.Common" Version="6.0.0" />
</ItemGroup> </ItemGroup>