在我的浅水实施中奇怪的振荡波纹

Question

我一直在尝试在Unity中实现浅水方程,但我遇到了一个奇怪的错误.我在水中得到了这些奇怪的摆动涟漪.我做了一些截图:

您可以在此处找到视频:https://www.youtube.com/watch？v = crXLrvETdjA

我的代码基于Xing Mei的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 out?ow ?ux
        float dV;

        for (x=1 ; x <= N ; x++ ) {
                for (y=1 ; y <= N ; y++ ) {
                        //
                        // 3.2.1 Out?ow 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);
}

Answer 1

我自己修好了!虽然是在开车送朋友的时候.问题很简单,我在错误的代码中做的是为每个单元格(或网格点)计算通量,然后是高度然后我去下一个单元格.我应该做的是首先计算所有细胞的通量,然后在所有细胞上再次迭代并计算它们的高度.所以代码变成:

for (x=1 ; x <= N ; x++ ) {
    for (y=1 ; y <= N ; y++ ) {
        //
        // 3.2.1 Out?ow Flux Computation
        // --------------------------------------------------------------
        ***
    }
}

for (x=1 ; x <= N ; x++ ) {
    for (y=1 ; y <= N ; y++ ) {
        //
        // 3.2.2 Water Surface
        // ---------------------------------------------------------------------------
        ***
    }
}
Helpers.Swap(ref _tempFlux, ref _flux);
Helpers.Swap(ref _tempHeight, ref _height);

(当然***成为上述问题的相应代码.)

现在我有一个工作水模拟.