载入中...
搜索中...
未找到
FluidGpuKernels.h
浏览该文件的文档.
1#pragma once
2
36namespace eve::fluids {
37
39inline const char* kFluidClearGrid = R"GLSL(
40#version 450
42layout(local_size_x = 64) in;
44layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
46layout(push_constant) uniform PC { float d[32]; } pc;
48void main() {
49 uint i = gl_GlobalInvocationID.x;
50 uint resX = uint(pc.d[14]);
51 uint resY = uint(pc.d[15]);
52 uint resZ = uint(pc.d[16]);
53 if (i < resX * resY * resZ) head.h[i] = -1;
54}
55)GLSL";
56
58inline const char* kFluidBuildGrid = R"GLSL(
59#version 450
61layout(local_size_x = 64) in;
63layout(set = 0, binding = 0) buffer Pos { vec4 p[]; } pos;
65layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
67layout(set = 0, binding = 3) buffer CellNext { int nx[]; } next;
69layout(push_constant) uniform PC { float d[32]; } pc;
71void main() {
72 uint i = gl_GlobalInvocationID.x;
73 uint n = uint(pc.d[0]);
74 if (i >= n) return;
75 vec3 origin = vec3(pc.d[17], pc.d[18], pc.d[19]);
76 float cs = pc.d[20];
77 ivec3 dim = ivec3(int(pc.d[14]), int(pc.d[15]), int(pc.d[16]));
78 ivec3 c = ivec3(floor((pos.p[i].xyz - origin) / cs));
79 c = clamp(c, ivec3(0), dim - ivec3(1));
80 int cell = c.x + dim.x * (c.y + dim.y * c.z);
81 int prev = atomicExchange(head.h[cell], int(i));
82 next.nx[i] = prev;
83}
84)GLSL";
85
87inline const char* kFluidDensityLambda = R"GLSL(
88#version 450
90layout(local_size_x = 64) in;
92layout(set = 0, binding = 0) buffer Pos { vec4 p[]; } pos;
94layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
96layout(set = 0, binding = 3) buffer CellNext { int nx[]; } next;
98layout(set = 0, binding = 4) buffer Dens { float d[]; } dens;
100layout(set = 0, binding = 6) buffer Lambda { float l[]; } lam;
102layout(set = 0, binding = 7) buffer Grad { vec4 g[]; } grad;
104layout(push_constant) uniform PC { float d[32]; } pc;
105
107float poly6(float r2, float h) {
108 float h2 = h * h;
109 if (r2 >= h2 || r2 <= 0.0) return 0.0;
110 float q = h2 - r2;
111 return 315.0 / (64.0 * 3.14159265358979 * pow(h, 9.0)) * q * q * q;
112}
113
115vec3 spikyGrad(vec3 dx, float h) {
116 float r2 = dot(dx, dx);
117 if (r2 >= h * h || r2 <= 1e-12) return vec3(0.0);
118 float r = sqrt(r2);
119 float q = h - r;
120 float k = -45.0 / (3.14159265358979 * pow(h, 6.0));
121 return dx * (k * q * q / r);
122}
123
125void main() {
126 uint i = gl_GlobalInvocationID.x;
127 uint n = uint(pc.d[0]);
128 if (i >= n) return;
129 float h = pc.d[3];
130 float h2 = h * h;
131 vec3 pi = pos.p[i].xyz;
132 vec3 origin = vec3(pc.d[17], pc.d[18], pc.d[19]);
133 float cs = pc.d[20];
134 ivec3 dim = ivec3(int(pc.d[14]), int(pc.d[15]), int(pc.d[16]));
135 ivec3 c = ivec3(floor((pi - origin) / cs));
136 c = clamp(c, ivec3(0), dim - ivec3(1));
137 float rho = 0.0;
138 vec3 gradSum = vec3(0.0);
139 for (int oz = -1; oz <= 1; ++oz) {
140 for (int oy = -1; oy <= 1; ++oy) {
141 for (int ox = -1; ox <= 1; ++ox) {
142 ivec3 cc = c + ivec3(ox, oy, oz);
143 if (any(lessThan(cc, ivec3(0))) || any(greaterThanEqual(cc, dim))) continue;
144 int cell = cc.x + dim.x * (cc.y + dim.y * cc.z);
145 for (int j = head.h[cell]; j >= 0; j = next.nx[j]) {
146 vec3 dx = pi - pos.p[j].xyz;
147 float r2 = dot(dx, dx);
148 if (r2 < h2) {
149 rho += poly6(r2, h);
150 gradSum += spikyGrad(dx, h);
151 }
152 }
153 }
154 }
155 }
156 dens.d[i] = rho;
157 grad.g[i] = vec4(gradSum, 0.0);
158 float rho0 = max(pc.d[4], 1e-6);
159 float C = max(0.0, rho / rho0 - 1.0);
160 lam.l[i] = -C / (dot(gradSum, gradSum) + 1e-6);
161}
162)GLSL";
163
165inline const char* kFluidComputeDelta = R"GLSL(
166#version 450
168layout(local_size_x = 64) in;
170layout(set = 0, binding = 0) buffer Pos { vec4 p[]; } pos;
172layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
174layout(set = 0, binding = 3) buffer CellNext { int nx[]; } next;
176layout(set = 0, binding = 6) buffer Lambda { float l[]; } lam;
178layout(set = 0, binding = 7) buffer Grad { vec4 g[]; } grad;
180layout(push_constant) uniform PC { float d[32]; } pc;
181
183vec3 spikyGrad(vec3 dx, float h) {
184 float r2 = dot(dx, dx);
185 if (r2 >= h * h || r2 <= 1e-12) return vec3(0.0);
186 float r = sqrt(r2);
187 float q = h - r;
188 float k = -45.0 / (3.14159265358979 * pow(h, 6.0));
189 return dx * (k * q * q / r);
190}
191
193void main() {
194 uint i = gl_GlobalInvocationID.x;
195 uint n = uint(pc.d[0]);
196 if (i >= n) return;
197 float h = pc.d[3];
198 float h2 = h * h;
199 float rho0 = max(pc.d[4], 1e-6);
200 vec3 pi = pos.p[i].xyz;
201 float li = lam.l[i];
202 vec3 origin = vec3(pc.d[17], pc.d[18], pc.d[19]);
203 float cs = pc.d[20];
204 ivec3 dim = ivec3(int(pc.d[14]), int(pc.d[15]), int(pc.d[16]));
205 ivec3 c = ivec3(floor((pi - origin) / cs));
206 c = clamp(c, ivec3(0), dim - ivec3(1));
207 vec3 delta = vec3(0.0);
208 for (int oz = -1; oz <= 1; ++oz) {
209 for (int oy = -1; oy <= 1; ++oy) {
210 for (int ox = -1; ox <= 1; ++ox) {
211 ivec3 cc = c + ivec3(ox, oy, oz);
212 if (any(lessThan(cc, ivec3(0))) || any(greaterThanEqual(cc, dim))) continue;
213 int cell = cc.x + dim.x * (cc.y + dim.y * cc.z);
214 for (int j = head.h[cell]; j >= 0; j = next.nx[j]) {
215 if (j == int(i)) continue;
216 vec3 dx = pi - pos.p[j].xyz;
217 float r2 = dot(dx, dx);
218 if (r2 < h2) delta += (li + lam.l[j]) * spikyGrad(dx, h);
219 }
220 }
221 }
222 }
223 grad.g[i] = vec4(delta / rho0, 0.0);
224}
225)GLSL";
226
228inline const char* kFluidApplyDelta = R"GLSL(
229#version 450
231layout(local_size_x = 64) in;
233layout(set = 0, binding = 0) buffer Pos { vec4 p[]; } pos;
235layout(set = 0, binding = 1) buffer Vel { vec4 v[]; } vel;
237layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
239layout(set = 0, binding = 3) buffer CellNext { int nx[]; } next;
241layout(set = 0, binding = 5) buffer Sdf { float d[]; } sdf;
243layout(set = 0, binding = 7) buffer Grad { vec4 g[]; } grad;
245layout(push_constant) uniform PC { float d[32]; } pc;
246
248float sdfValue(ivec3 res, int x, int y, int z) {
249 int idx = x + res.x * (y + res.y * z);
250 return sdf.d[idx];
251}
252
254float sdfSample(vec3 p) {
255 ivec3 res = ivec3(int(pc.d[21]), int(pc.d[22]), int(pc.d[23]));
256 vec3 o = vec3(pc.d[24], pc.d[25], pc.d[26]);
257 float cs = pc.d[27];
258 vec3 f = (p - o) / cs;
259 vec3 clamped = clamp(f, vec3(0.0), vec3(res - 1));
260 float outside = length((f - clamped) * cs);
261 ivec3 i0 = ivec3(floor(clamped));
262 ivec3 i1 = min(i0 + ivec3(1), res - ivec3(1));
263 vec3 t = fract(clamped);
264 float v000 = sdfValue(res, i0.x, i0.y, i0.z);
265 float v100 = sdfValue(res, i1.x, i0.y, i0.z);
266 float v010 = sdfValue(res, i0.x, i1.y, i0.z);
267 float v110 = sdfValue(res, i1.x, i1.y, i0.z);
268 float v001 = sdfValue(res, i0.x, i0.y, i1.z);
269 float v101 = sdfValue(res, i1.x, i0.y, i1.z);
270 float v011 = sdfValue(res, i0.x, i1.y, i1.z);
271 float v111 = sdfValue(res, i1.x, i1.y, i1.z);
272 float c00 = v000 + (v100 - v000) * t.x;
273 float c10 = v010 + (v110 - v010) * t.x;
274 float c01 = v001 + (v101 - v001) * t.x;
275 float c11 = v011 + (v111 - v011) * t.x;
276 float c0 = c00 + (c10 - c00) * t.y;
277 float c1 = c01 + (c11 - c01) * t.y;
278 return c0 + (c1 - c0) * t.z + outside;
279}
280
282vec3 sdfGrad(vec3 p) {
283 float e = pc.d[27];
285 return vec3(sdfSample(p + vec3(e, 0.0, 0.0)) - sdfSample(p - vec3(e, 0.0, 0.0)),
287 sdfSample(p + vec3(0.0, e, 0.0)) - sdfSample(p - vec3(0.0, e, 0.0)),
289 sdfSample(p + vec3(0.0, 0.0, e)) - sdfSample(p - vec3(0.0, 0.0, e))) / (2.0 * e);
290}
291
293void main() {
294 uint i = gl_GlobalInvocationID.x;
295 uint n = uint(pc.d[0]);
296 if (i >= n) return;
297 if (pc.d[29] > 0.5) {
298 vec3 v = grad.g[i].xyz;
299 vec3 p = pos.p[i].xyz + v * pc.d[1];
300 float radius = pc.d[2];
301 float d = sdfSample(p);
302 vec3 nn = sdfGrad(p);
303 float nl = length(nn);
304 if (nl > 1e-6) nn /= nl;
305 else nn = vec3(0.0, 1.0, 0.0);
306 p -= nn * (d - radius);
307 v -= nn * dot(v, nn);
308 pos.p[i] = vec4(p, 0.0);
309 vel.v[i] = vec4(v, 0.0);
310 return;
311 }
312 vec3 delta = grad.g[i].xyz;
313 float maxD = 0.0005 * pc.d[3];
314 float dl = length(delta);
315 if (dl > maxD && dl > 1e-9) delta *= maxD / dl;
316 vec3 p = pos.p[i].xyz + delta;
317 float radius = pc.d[2];
318 float d = sdfSample(p);
319 {
320 vec3 nn = sdfGrad(p);
321 float nl = length(nn);
322 if (nl > 1e-6) nn /= nl;
323 else nn = vec3(0.0, 1.0, 0.0);
324 p = p - nn * (d - radius);
325 }
326 pos.p[i] = vec4(p, 0.0);
327}
328)GLSL";
329
331inline const char* kFluidIntegrate = R"GLSL(
332#version 450
334layout(local_size_x = 64) in;
336layout(set = 0, binding = 0) buffer Pos { vec4 p[]; } pos;
338layout(set = 0, binding = 1) buffer Vel { vec4 v[]; } vel;
340layout(set = 0, binding = 2) buffer CellHead { int h[]; } head;
342layout(set = 0, binding = 3) buffer CellNext { int nx[]; } next;
344layout(set = 0, binding = 5) buffer Sdf { float d[]; } sdf;
346layout(set = 0, binding = 7) buffer Grad { vec4 g[]; } grad;
348layout(push_constant) uniform PC { float d[32]; } pc;
349
351float poly6(float r2, float h) {
352 float h2 = h * h;
353 if (r2 >= h2 || r2 <= 0.0) return 0.0;
354 float q = h2 - r2;
355 return 315.0 / (64.0 * 3.14159265358979 * pow(h, 9.0)) * q * q * q;
356}
357
359float cohesionKernel(float r, float h) {
360 if (r >= h || r <= 1e-9) return 0.0;
361 float q = h - r;
362 return 32.0 / (3.14159265358979 * pow(h, 9.0)) * q * q * q * r * r * r;
363}
364
366float sdfValue(ivec3 res, int x, int y, int z) {
367 int idx = x + res.x * (y + res.y * z);
368 return sdf.d[idx];
369}
370
372float sdfSample(vec3 p) {
373 ivec3 res = ivec3(int(pc.d[21]), int(pc.d[22]), int(pc.d[23]));
374 vec3 o = vec3(pc.d[24], pc.d[25], pc.d[26]);
375 float cs = pc.d[27];
376 vec3 f = (p - o) / cs;
377 vec3 clamped = clamp(f, vec3(0.0), vec3(res - 1));
378 float outside = length((f - clamped) * cs);
379 ivec3 i0 = ivec3(floor(clamped));
380 ivec3 i1 = min(i0 + ivec3(1), res - ivec3(1));
381 vec3 t = fract(clamped);
382 float v000 = sdfValue(res, i0.x, i0.y, i0.z);
383 float v100 = sdfValue(res, i1.x, i0.y, i0.z);
384 float v010 = sdfValue(res, i0.x, i1.y, i0.z);
385 float v110 = sdfValue(res, i1.x, i1.y, i0.z);
386 float v001 = sdfValue(res, i0.x, i0.y, i1.z);
387 float v101 = sdfValue(res, i1.x, i0.y, i1.z);
388 float v011 = sdfValue(res, i0.x, i1.y, i1.z);
389 float v111 = sdfValue(res, i1.x, i1.y, i1.z);
390 float c00 = v000 + (v100 - v000) * t.x;
391 float c10 = v010 + (v110 - v010) * t.x;
392 float c01 = v001 + (v101 - v001) * t.x;
393 float c11 = v011 + (v111 - v011) * t.x;
394 float c0 = c00 + (c10 - c00) * t.y;
395 float c1 = c01 + (c11 - c01) * t.y;
396 return c0 + (c1 - c0) * t.z + outside;
397}
398
400vec3 sdfGrad(vec3 p) {
401 float e = pc.d[27];
403 return vec3(sdfSample(p + vec3(e, 0.0, 0.0)) - sdfSample(p - vec3(e, 0.0, 0.0)),
405 sdfSample(p + vec3(0.0, e, 0.0)) - sdfSample(p - vec3(0.0, e, 0.0)),
407 sdfSample(p + vec3(0.0, 0.0, e)) - sdfSample(p - vec3(0.0, 0.0, e))) / (2.0 * e);
408}
409
411void main() {
412 uint i = gl_GlobalInvocationID.x;
413 uint n = uint(pc.d[0]);
414 if (i >= n) return;
415 float dt = pc.d[1];
416 float h = pc.d[3];
417 float h2 = h * h;
418 vec3 p = pos.p[i].xyz;
419 vec3 v = vel.v[i].xyz;
420
421 float visc = pc.d[8];
422 float yield = pc.d[9];
423 float coh = pc.d[10];
424 if (visc > 0.0 || coh > 0.0 || yield > 0.0) {
425 vec3 viscAcc = vec3(0.0);
426 vec3 cohAcc = vec3(0.0);
427 float shearAcc = 0.0;
428 float neighborWeight = 0.0;
429 float cohesionWeight = 0.0;
430 float particleVolume = pow(2.0 * pc.d[2], 3.0);
431 vec3 origin = vec3(pc.d[17], pc.d[18], pc.d[19]);
432 float cs = pc.d[20];
433 ivec3 dim = ivec3(int(pc.d[14]), int(pc.d[15]), int(pc.d[16]));
434 ivec3 c = ivec3(floor((p - origin) / cs));
435 c = clamp(c, ivec3(0), dim - ivec3(1));
436 for (int oz = -1; oz <= 1; ++oz) {
437 for (int oy = -1; oy <= 1; ++oy) {
438 for (int ox = -1; ox <= 1; ++ox) {
439 ivec3 cc = c + ivec3(ox, oy, oz);
440 if (any(lessThan(cc, ivec3(0))) || any(greaterThanEqual(cc, dim))) continue;
441 int cell = cc.x + dim.x * (cc.y + dim.y * cc.z);
442 for (int j = head.h[cell]; j >= 0; j = next.nx[j]) {
443 if (j == int(i)) continue;
444 vec3 dx = p - pos.p[j].xyz;
445 float r2 = dot(dx, dx);
446 if (r2 < h2) {
447 float weight = poly6(r2, h) * particleVolume;
448 if (visc > 0.0 || yield > 0.0) {
449 viscAcc += (vel.v[j].xyz - v) * weight;
450 neighborWeight += weight;
451 }
452 if (yield > 0.0) shearAcc += length(vel.v[j].xyz - v) * weight;
453 if (coh > 0.0) {
454 cohAcc -= dx * weight;
455 cohesionWeight += weight;
456 }
457 }
458 }
459 }
460 }
461 }
462 float effVisc = visc;
463 if (yield > 0.0 && shearAcc > 1e-6) effVisc += yield / shearAcc;
464 if (neighborWeight > 1e-6)
465 v += (viscAcc / neighborWeight) * clamp(effVisc * dt, 0.0, 1.0);
466 vec3 surfaceNormal = sdfGrad(p);
467 if (length(surfaceNormal) > 1e-6) surfaceNormal = normalize(surfaceNormal);
468 else surfaceNormal = vec3(0.0, 1.0, 0.0);
469 cohAcc -= surfaceNormal * dot(cohAcc, surfaceNormal);
470 vec3 cohesionDv = cohesionWeight > 1e-6 ?
471 (cohAcc / cohesionWeight) * (coh * 100.0 * dt) : vec3(0.0);
472 float cohesionSpeed = length(cohesionDv);
473 float maxCohesionDv = pc.d[2];
474 if (cohesionSpeed > maxCohesionDv && cohesionSpeed > 1e-9)
475 cohesionDv *= maxCohesionDv / cohesionSpeed;
476 v += cohesionDv;
477 }
478
479 v += vec3(pc.d[5], pc.d[6], pc.d[7]) * dt;
480 v *= max(0.0, 1.0 - pc.d[12] * dt);
481 if (pc.d[11] > 0.0) v *= exp(-pc.d[11] * 300.0 * dt);
482 float sp = length(v);
483 float vmax = pc.d[13];
484 if (sp > vmax && sp > 1e-9) v *= vmax / sp;
485 grad.g[i] = vec4(v, 0.0);
486}
487)GLSL";
488
489} // namespace eve::fluids
eve::fluids
GLSL compute kernels for the GPU surface-flow solver.
Definition FluidTarget.h:12
eve::fluids::kFluidClearGrid
const char * kFluidClearGrid
Zero the linked-list cell heads.
Definition FluidGpuKernels.h:39
eve::fluids::kFluidComputeDelta
const char * kFluidComputeDelta
Accumulate PBF position deltas into the (reused) grad buffer.
Definition FluidGpuKernels.h:165
eve::fluids::kFluidApplyDelta
const char * kFluidApplyDelta
Apply PBF deltas and re-project onto the SDF surface.
Definition FluidGpuKernels.h:228
eve::fluids::kFluidDensityLambda
const char * kFluidDensityLambda
Density, gradient sum and PBF lambda in one pass (mirror computeDensitiesAndGrads + computeLambdas).
Definition FluidGpuKernels.h:87
eve::fluids::kFluidIntegrate
const char * kFluidIntegrate
Viscosity + cohesion + adhesion + gravity + integration + SDF projection.
Definition FluidGpuKernels.h:331
eve::fluids::kFluidBuildGrid
const char * kFluidBuildGrid
Insert every particle into its grid cell (linked list via atomicExchange).
Definition FluidGpuKernels.h:58