C++风流和MATLAB | Python | CUDA 库埃特流泊肃叶流薄膜流体

🎯要点

🍇C++ GPU加速模拟三维风流

在物理学中，流体是一种在施加的剪切应力或外力作用下可能连续移动和变形（流动）的液体、气体或其他材料。它们的剪切模量为零，或者简单地说，是无法抵抗施加在其上的任何剪切力的物质。

虽然流体这一术语通常包括液相和气相，但其定义在不同的科学分支中有所不同。固体的定义也各不相同，根据领域的不同，一些物质可以同时具有流体和固体的性质。非牛顿流体，如橡皮泥，在突然施加力时似乎表现得与固体相似。粘度极高的物质，如沥青，似乎也表现得像固体（见沥青滴落实验）。在粒子物理学中，这一概念扩展到包括液体或气体以外的流体物质。[4] 医学或生物学中的流体是指身体的任何液体成分（体液），而“液体”不是以这种意义使用的。有时，通过饮用或注射来补充体液的液体也被称为流体（例如“多喝水”）。在液压学中，流体是指具有某些特性的液体，比（液压）油更广泛。

相比之下，固体可以抵抗“变形”（例如，通过产生与剪切力成比例的弹性恢复力），而流体则不会并且只会变形。根据这个定义，空气是一种流体，因为它的行为就像由其变形行为定义的流体一样，因此我们可以将流体力学定律应用于它。

偏微分方程将未知多变量函数的偏导数相互关联，从而描述其在各个维度（通常是空间、时间）的行为并隐式定义该函数。偏微分方程可以解释为一个约束，求解偏微分方程意味着找到满足该约束的函数。

将任何物理守恒量写为时间和空间的函数，该函数必须满足微分守恒定律，偏微分方程约束才能守恒。

最基本的微分守恒定律之一是“平流扩散方程”，它描述了运动流体中扩散的稀释溶质的质量如何守恒。
$\frac{\partial c}{\partial t}=\nabla \cdot(D \nabla c)-\nabla \cdot(v c)+R$
这表明溶质浓度对时间（在固定位置）的导数由三项之和给出：扩散项、平流项和源项。将源 R 设置为 0 意味着质量进入或离开体积的唯一方式是通过扩散或平流穿过微分体积边界，从而守恒质量。

对偏微分方程（尤其是纳维-斯托克斯方程）进行分析积分通常仅适用于具有简单边界和初始条件的少数众所周知的问题陈述。计算流体动力学关注的是流体流动问题的数值求解，并且大量研究集中于提高这些偏微分方程求解方法的速度、精度和稳定性。

因此，粒子分布函数和传播缓冲区的大小为 $N^d*Q$ ，而宏观量和边界条件的大小为 $N^d$ 。对于我们的实现，我们只需在 CPU 上声明具有正确大小的缓冲区。请注意，我们在着色器代码中使用扁平索引对这些缓冲区进行索引。

Buffer f(NX*NY*Q, (float*)NULL);      
Buffer fprop(NX*NY*Q, (float*)NULL);  Buffer rho(NX*NY, (float*)NULL);
Buffer v(NX*NY, (glm::vec4*)NULL);Buffer b(NX*NY, (float*)NULL);

然后，我们将这些缓冲区作为着色器存储缓冲区对象绑定到我们的主着色器：

Compute init("shader/init.cs", {"f", "fprop", "b", "rho", "v"});
init.bind<float>("f", &f);
init.bind<float>("fprop", &fprop);
init.bind<float>("b", &b);
init.bind<float>("rho", &rho);
init.bind<glm::vec4>("v", &dirbuf);Compute collide("shader/collide.cs", {"f", "fprop", "b", "rho", "v"});
collide.bind<float>("f", &f);
collide.bind<float>("fprop", &fprop);
collide.bind<float>("b", &b);
collide.bind<float>("rho", &rho);
collide.bind<glm::vec4>("v", &dirbuf);Compute stream("shader/stream.cs", {"f", "fprop", "b"});
stream.bind<float>("f", &f);
stream.bind<float>("fprop", &fprop);
stream.bind<float>("b", &b);

并使用以下方式在着色器中访问它们：

layout (std430, binding = 0) buffer f {float F[];
};layout (std430, binding = 1) buffer fprop {float FPROP[];
};layout (std430, binding = 2) buffer b {float B[];
};

#version 460 core
layout(local_size_x = 16, local_size_y = 1, local_size_z = 16) in;#include lbm.cslayout (std430, binding = 3) buffer rho {float RHO[];
};layout (std430, binding = 4) buffer v {vec4 V[];
};

#version 460 corelayout(local_size_x = 16, local_size_y = 1, local_size_z = 16) in;#include lbm.cslayout (std430, binding = 3) buffer rho {float RHO[];
};

二维和三维速度集合

uniform int NX = 256;
uniform int NY = 128;
const int Q = 9;
const float w[Q] = {4.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/9.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0};
const int cp[Q] = {0, 3, 4, 1, 2, 7, 8, 5, 6};
const ivec2 c[Q] = {ivec2(0, 0),ivec2(1, 0),ivec2(0, 1),ivec2(-1, 0),ivec2(0, -1),ivec2(1, 1),ivec2(-1, 1),ivec2(-1, -1),ivec2(1, -1)
};

三维速度集

uniform int NX = 64;
uniform int NY = 64;
uniform int NZ = 64;
const int Q = 19;const float w[Q] = {1.0/3.0,1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0, 1.0/18.0,1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0,1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0, 1.0/36.0
};const ivec3 c[Q] = {ivec3( 0,  0,  0),ivec3( 1,  0,  0), ivec3(-1,  0,  0),ivec3( 0,  1,  0), ivec3( 0, -1,  0),ivec3( 0,  0,  1), ivec3( 0,  0, -1),ivec3( 1,  1,  0), ivec3(-1, -1,  0),ivec3( 1,  0,  1), ivec3(-1,  0, -1),ivec3( 0,  1,  1), ivec3( 0, -1, -1),ivec3( 1, -1,  0), ivec3(-1,  1,  0),ivec3( 1,  0, -1), ivec3(-1,  0,  1),ivec3( 0,  1, -1), ivec3( 0, -1,  1)
};const int cp[Q] = {0,2, 1, 4, 3, 6, 5,8, 7, 10, 9, 12, 11,14, 13, 16, 15, 18, 17
};