GAMES103 Lab4 Shallow Wave | 路人乙の小窝

路人乙の小窝

我若安好，万事大吉

0%

GAMES103 Lab4 Shallow Wave

发表于

前言

最近在学习流体仿真部分，主要是GAMES103相关部分。这篇文章作为学习 Shallow Wave 内容的总结和作业4的展示。

算法背景

Shallow wave 模型用于模拟水池中涟漪、湖面波动。模型把浅水想象成一个有弹性的水面网格，离散化后存储在二维数组中，每个格子只记录当前时刻和前一时刻的高度，并根据这些高度信息计算出新的高度。

更详细的说，对于每个点 $\mathbf{p_i}$ ，模型用当前点的高度信息 $h_i$ 和历史高度 $old\_h_j$，和邻居点 $\mathbf{p_j}$ 的高度信息 $h_j$ 计算得到 $new\_h_i$ ，进行和刚体的交互更新，最终应用到网格。

实验环境

本实验在unity环境下进行仿真，使用c#编写脚本。课程框架由GAMES103课程官网给出。

实验内容

获取 $h_{ij}$

算法将三维的水面抽象成二维的网格，关于网格中的每个点，我们最关注他的高度信息$h_{ij}$。

在框架代码中，提供了一维数组X存储网格高度，为了方便计算，我们需要定义二维数组h[i, j]，并在每次 Update 中将 X 中的值载入 h，处理后再返还到 X。

根据X的初始化代码：

1
2
3
4
5
6
7
for (int i=0; i<size; i++)
for (int j=0; j<size; j++) 
{
    X[i*size+j].x=i*0.1f-size*0.05f;
    X[i*size+j].y=0;
    X[i*size+j].z=j*0.1f-size*0.05f;
}

可以看到 X 坐标的处理方式，我们只需要照葫芦画瓢即可进行 h 数组的读取：

1
2
3
4
// Load X.y into h.
for (int i = 0; i < size; i++)
    for(int j = 0; j < size; j++)
        h[i, j] = X[i * size + j].y;

生成随机水花

为了体现算法效果，我们需要提供一个方法对水面产生随机扰动，具体方法是在某随机一位置生成随机高度的液体。

为了保证水的总质量不变，我们需要在邻居处减去总共相同高度的水。

具体实现:

1
2
3
4
5
6
7
8
9
10
11
12
if (Input.GetKeyDown ("r")) 
	{
		// Add random water.
		int i = Random.Range(1, size - 1);
		int j = Random.Range(1, size - 1);
		float R = Random.Range(0.1f, 1f);
		h[i, j] += R;
		for(int k = 0; k < 4; ++k)
		{
			h[i + dx[k], j + dy[k]] -= R / 4;	
		}
	}

里面有一个细节，为了防止数组越界，我把 i 和 j 的位置限制在了$[1, size - 2]$，这样邻居的位置就限定在了 $[0, size-1]$ ，正好避免了越界问题。dx 和 dy 是我定义的方向数组，增强访问邻居部分代码的可读性。

更新网格节点坐标

先跳过具体的算法部分。假设我们已经完成了h[i, j] 的计算，下一步要做的就是用其中的值更新X，并将X作为新的网格的mesh.vertices，并重新计算法线：

1
2
3
4
5
6
7
8
9
10
// Store h back into X.y and recalculate normal.
for (int i = 0; i < size; ++i)
{
	for(int j = 0; j < size; ++j)
	{
		X[i * size + j].y = h[i , j];
	}
}
mesh.vertices = X;
mesh.RecalculateNormals();

Shallow Wave 模型

完成了上面部分的工作，我们对Shallow_Wave(float[,] old_h, float[,] h, float [,] new_h)函数的要实现的结果已经比较清楚：计算new_h，用 h 和 new_h 更新 old_h。

计算 new_h

二维波动方程：
$$
\frac{\partial^2h}{\partial t^2} = c^2\nabla^2h
$$
对时间二阶导数进行中心差分：
$$
\frac{\partial^2h}{\partial t^2} = \frac{h^{n+1} - 2h^n + h^{n-1}}{\Delta t^2}
$$
代入：
$$
\frac{h^{n+1} - 2h^n + h^{n-1}}{\Delta t^2} = c^2\nabla^2h^n
$$
整理得到
$$
h^{n+1} = 2h^n - h^{n-1} + c^2\Delta t^2 \nabla^2 h^n
$$
即
$$
h^{n+1} = h^n + (h^n - h^{n-1}) + c^2\Delta t^2 \nabla^2 h^n
$$
由于我们并不希望算法在一次处理后将高度恢复到理想位置，而是希望产生类似涟漪的效果，所以我们添加了阻尼系数damping。同时，我们把常数 $(\frac{c\Delta t}{\Delta x})^2$ 记作 rate ，得到程序中的代码：

1
new_h[i, j] = h[i, j] + damping * (h[i, j] - old_h[i, j]) + rate * laplacian;

由于我的物理功底较差，我并不完全理解这些数学公式。但从直观上看：

  • h[i, j] 为上一时刻速度；
  • damping * (h[i, j] - old_h[i, j]) 可以看做一个惯性项，同时这个惯性随着阻力削减；
  • rate * laplacian 为邻居对当前位置高度的影响。

由于damping和 rate都由程序框架给出，尽管我只有一个直观的感受而没有数学上的理解，我也可以做出比较理想的效果。

其中laplacian为邻居高度之和减去邻居数量乘当前高度。一般来说，邻居数量为四。对于边界处，程序使用Neumann边界条件，认为边界内外高度相同。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// Step 1:
// Compute new_h based on the shallow wave model.
for(int i = 0; i < size; ++i)
{
	for(int j = 0; j < size; ++j)
	{
		float laplacian = 0f;
		for(int k = 0; k < 4; ++k)
		{
			int x = i + dx[k];
			int y = j + dy[k];
			if (x < 0 || y < 0 || x == size || y == size) continue;
			laplacian += h[x, y] - h[i, j];
		}
		new_h[i, j] = h[i, j] + damping * (h[i, j] - old_h[i, j]) + rate * laplacian;
	}
}

Block -> Water coupling

Unity场景中有两个刚体立方体，程序需要考虑方块对水的影响。

首先，两个BoxCollider可以存在一个数组中方便访问，用GameObject.Find("Cube").GetComponent<BoxCollider>()获取。用unity BoxCollider类自带的bounds属性，可以得到立方体影响到的水面范围，这个范围可以参考 Start() 中 X的初始化代码进行坐标变换。得到影响范围可以用它更新cg_mask，方便后续进行求解。

水面高度会被立方体的底部位置限制，对立方体影响到的每一个网格，设其low_h = bound.min.y，表示其水面高度被限制。据此，我们可以计算当前格子需要的高度修正量
$$
b = \frac{new\_h_{i,j} - low\_h_{i, j}}{rate}
$$
需要除以rate是因为在水面更新时会乘一次rate，这里提前除以rate进行抵消。

为了让液体按照我们的预想移动，我们可以为其添加一层“虚拟高度” $vh$ 。

参考前面的公式和代码，最终的水面修正量：
$$
\Delta h = rate \times \nabla_d^2vh = low_h - new_h
$$
两边除以 rate：
$$
\nabla^2_dvh = - \frac{new_h - low_h}{rate} = -b
$$
令$A = -\nabla_d^2$ ，得到$Avh = b$，即Poisson eqiation的离散形式，可以用CG求解器求解vh。

在得到vh后，我们同样不希望水面高度在一瞬间内完成变化，因此在高度修正时添加了一个系数gamma。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
// Step 2: Block->Water coupling
// for block 1, calculate low_h.
// then set up b and cg_mask for conjugate gradient.
// Solve the Poisson equation to obtain vh (virtual height).

// for block 2, calculate low_h.
// then set up b and cg_mask for conjugate gradient.
// Solve the Poisson equation to obtain vh (virtual height).

// Diminish vh.

// Update new_h by vh.
for(int i = 0; i < size; ++i) {
    for(int j = 0; j < size; ++j) {
        low_h[i, j] = float.PositiveInfinity;
        cg_mask[i, j] = false;
        b[i, j] = 0.0f;
        vh[i, j] = 0.0f;
    }
}

for(int temp = 0; temp < 2; ++temp) {
    Bounds bound = cubes[temp].bounds;
    int li = Mathf.CeilToInt(
        (bound.min.x + size * 0.05f) / 0.1f
    );

    int ui = Mathf.FloorToInt(
        (bound.max.x + size * 0.05f) / 0.1f
    );

    int lj = Mathf.CeilToInt(
        (bound.min.z + size * 0.05f) / 0.1f
    );

    int uj = Mathf.FloorToInt(
        (bound.max.z + size * 0.05f) / 0.1f
    );

    li = Mathf.Clamp(li, 0, size - 1);
    ui = Mathf.Clamp(ui, 0, size - 1);
    lj = Mathf.Clamp(lj, 0, size - 1);
    uj = Mathf.Clamp(uj, 0, size - 1);

    float bottom = bound.min.y;

    for (int i = li; i <= ui; ++i) {
        for (int j = lj; j <= uj; ++j) {
            low_h[i, j] = bottom;
            cg_mask[i, j] = true;
            b[i, j] = (new_h[i, j] - low_h[i, j]) / rate;
        }
    }

    Conjugate_Gradient(cg_mask, b, vh, li, ui, lj, uj);
}

for(int i = 0; i < size; ++i)
    for(int j  = 0; j < size; ++j)
        vh[i, j] *= gamma;

for (int i = 0; i < size; ++i) {
    for (int j = 0; j < size; ++j) {
        for (int k = 0; k < 4; ++k) {
            int x = i + dx[k];
            int y = j + dy[k];
            if (x < 0 || y < 0 || x >= size || y >= size) { continue; }
            new_h[i, j] += rate * (vh[x, y] - vh[i, j]);
        }
    }
}

更新高度

在得到new_h后，old_h和h的修改就很简单明了了。

1
2
3
4
5
6
7
8
9
10
11
// Step 3
// old_h <- h; h <- new_h;
// old_h = h; h = new_h;
for (int i = 0; i < size; ++i)
{
    for(int j =  0; j < size; ++j)
    {
        old_h[i, j] = h[i, j];
        h[i, j] = new_h[i, j];
    }
}

TODO: Water->Block coupling

作为作业的bonus部分，要求考虑水面对刚体方块的影响，目前我还没有进行实现。

实验结果

image-20260721204311450

程序没有报错，Unity成功运行，可以用鼠标拖动方块产生涟漪，同样可以点按r键产生水花。

之后的工作

GAMES103的流体部分已经学完，打算之后一段时间继续学习流体有关内容。Bonus部分作业也会早日完成。

我很可爱，请给我钱qwq