我一直在尝试在 Unity 中实现浅水方程,但我遇到了一个奇怪的错误。我的水中出现了这些奇怪的振荡涟漪。我做了一些截图:
您可以在这里找到一个视频:https ://www.youtube.com/watch?v=crXLrvETdjA
我的代码基于 Xing Mei 的论文Fast Hydraulic Erosion Simulation and Visualization on GPU。您可以在此处找到完整的求解器代码:http: //pastebin.com/JktpizHW(或见下文。)每次我使用论文中的公式时,我都会添加其编号作为注释。
我尝试了不同的时间步长,对于我使用 0.02 的视频,降低它只会使其振荡更慢。我还尝试了一个更大的网格(视频使用 100,我尝试了 200,但后来涟漪变小了。)我检查了所有公式几次,找不到任何错误。
这里有谁能弄清楚出了什么问题?
额外信息:
正如您从 pastebin 中看到的那样,我用 c# 对其进行了编程。我使用 Unity 作为我的可视化引擎,我只是使用网格来可视化水。我改变网格的顶点 y 分量以匹配我计算的高度。
DoUpdate 方法获取一个float[][] lowerLayersHeight
参数,它基本上是水下地形的高度。在视频中,这就是全部0
。
public override void DoUpdate(float dt, float dx, float[][] lowerLayersHeight) {
int x, y;
float totalHeight, dhL, dhR, dhT, dhB;
float dt_A_g_l = dt * _A * g / dx; //all constants for equation 2
float K; // scaling factor for the outflow flux
float dV;
for (x=1 ; x <= N ; x++ ) {
for (y=1 ; y <= N ; y++ ) {
//
// 3.2.1 Outflow Flux Computation
// --------------------------------------------------------------
totalHeight = lowerLayersHeight[x][y] + _height[x][y];
dhL = totalHeight - lowerLayersHeight[x-1][y] - _height[x-1][y]; //(3)
dhR = totalHeight - lowerLayersHeight[x+1][y] - _height[x+1][y];
dhT = totalHeight - lowerLayersHeight[x][y+1] - _height[x][y+1];
dhB = totalHeight - lowerLayersHeight[x][y-1] - _height[x][y-1];
_tempFlux[x][y].left = Mathf.Max(0.0f, _flux[x][y].left + dt_A_g_l * dhL ); //(2)
_tempFlux[x][y].right = Mathf.Max(0.0f, _flux[x][y].right + dt_A_g_l * dhR );
_tempFlux[x][y].top = Mathf.Max(0.0f, _flux[x][y].top + dt_A_g_l * dhT );
_tempFlux[x][y].bottom = Mathf.Max(0.0f, _flux[x][y].bottom + dt_A_g_l * dhB );
float totalFlux = _tempFlux[x][y].left + _tempFlux[x][y].right + _tempFlux[x][y].top + _tempFlux[x][y].bottom;
if (totalFlux > 0) {
K = Mathf.Min(1.0f, _height[x][y] * dx * dx / totalFlux / dt); //(4)
_tempFlux[x][y].left = K * _tempFlux[x][y].left; //(5)
_tempFlux[x][y].right = K * _tempFlux[x][y].right;
_tempFlux[x][y].top = K * _tempFlux[x][y].top;
_tempFlux[x][y].bottom = K * _tempFlux[x][y].bottom;
}
//swap temp and the real one after the for-loops
//
// 3.2.2 Water Surface
// ----------------------------------------------------------------------------------------
dV = dt * (
//sum in
_tempFlux[x-1][y].right + _tempFlux[x][y-1].top + _tempFlux[x+1][y].left + _tempFlux[x][y+1].bottom
//minus sum out
- _tempFlux[x][y].right - _tempFlux[x][y].top - _tempFlux[x][y].left - _tempFlux[x][y].bottom
); //(6)
_tempHeight[x][y] = _height[x][y] + dV / (dx*dx); //(7)
//swap temp and the real one after the for-loops
}
}
Helpers.Swap(ref _tempFlux, ref _flux);
Helpers.Swap(ref _tempHeight, ref _height);
}