除第8节：“质量守恒”外，大多数情况下所有方法都有效。在那里，我们建立了一个方程矩阵，以计算使液体自由扩散所需的压力。
我相信我的代码与本文相符，但是在质量步长守恒过程中得到了不可解的矩阵。
这是生成矩阵A的步骤：
将对角线条目$ A_ {i，i} $设置为与单元格i相邻的液体单元数的负数。
设置条目$ A_ {i，j} $和$ A_ {j，i} $设为1（如果单元格i和j都为液体）。
请注意，在我的实现中，液体网格中的单元格$ i $，$ j $对应于行矩阵中的$ i + $ gridWidth $ * j $。
论文提到：“静态对象和空单元格不会破坏此结构。在这种情况下，压力和速度项可以从两侧消失
”，所以我删除了没有液体的单元格的行和列。
所以我的问题是：为什么矩阵是奇异的？我是否在论文的其他地方缺少某种边界条件？我的实现是2D的事实吗？
这是我的实现中的一个2x2网格示例矩阵，其中0,0处的像元没有液体：
编辑
我的研究使我相信我没有正确处理边界条件。
首先，此时我可以说我的矩阵代表离散压力泊松方程。它是应用拉普拉斯算子将局部压力变化耦合到细胞发散的离散模拟。
据我了解，由于我们正在处理压力差，因此需要边界条件来将压力“锚固”到绝对参考值。否则，方程组可能有无数个解。
在我的理解中，给出了三种不同的边界条件应用方法：
Dirichlet-指定边界处的绝对值。
Neummann-指定边界处的导数。
Robin-指定边界处的绝对值和导数的某种线性组合。
Foster和Fedki的论文没有提及其中任何一个，但我相信它们会强制执行Dirichlet边界条件，因为在7.1.2末尾，设置为大气压力。“我们如何严格执行这些边界条件？看看其他实现，似乎有某种“ Ghost”单元位于边界的概念。
在这里，我已经链接到一些可能对其他人有所帮助的资料。关于Neumann边界条件的Science StackExchange文章
关于Poisson解算器的计算科学StackExchange文章
水物理棒的实现
这是代码I用于生成矩阵。请注意，我没有显式删除列和行，而是生成并使用从液体单元格索引到最终矩阵列/行的映射。
-1 0 1
0 -1 1
1 1 -2
#1 楼
从您的代码片段和2x2示例的结果中，我可以看到您实际上是在模拟仅包含Neumann边界条件（滑移墙）的域。在这种情况下，系统包含一个空空间，并且您的矩阵是奇异的。解决方案中的空间。如果您按照该论文的建议使用共轭梯度（CG），这将很简单。在CG迭代的每个迭代中，只需取当前解向量$ \ mathbf {x} $，然后做< > $$
\ begin {align}
\ mathbf {x}'＆=（\ mathbf {I}-\ mathbf {\ hat {u}} \ mathbf {\ hat {u}} ^ T）\ mathbf {x} \\
＆= \ mathbf {x}-（\ mathbf {\ hat {u}} \ cdot \ mathbf {x}）\ mathbf {\ hat {u}}
\ end {align}
$$
其中$ \ mathbf {\ hat {u}} $是梯度算子的规范化零空间：
$ \ mathbf {u} =（ 1,1，\ dots，1）$，$ \ mathbf {\ hat {u}} = \ frac {\ mathbf {u}} {\ lVert \ mathbf {u} \ rVert} $。
否则，如果您想模拟空气（自由边界或Dirichlet BC），则需要区分墙壁和空气室（即仅具有布尔
hasLiquid是不够的），并对其进行正确的离散化（请参见下文） 。最后一点，对角线输入为负。您可能需要翻转符号以便CG方法起作用。
下面，我想显示更多详细信息。
考虑压力投影过程。将压力投影之前的速度表示为$ \ mathbf {v} ^ * $。它可能是发散的，因此我们计算压力以对其进行校正，并获得发散自由速度$ \ mathbf {v} ^ {n + 1} $。即，
$$
\开始{align}
\ mathbf {v} ^ {n + 1}＆= \ mathbf {v} ^ *-\ frac {\ Delta t} {\ rho} \ nabla P
\ end {align}
$$
求散度，由于$ \ mathbf {v} ^ {n + 1} $是无散度的，
$$
\开始{align}
\ nabla \ cdot \ nabla P＆= \ nabla \ cdot \ mathbf {v} ^ *
\ end {align}
$$
假设没有压力Dirichlet BC存在并且我们有一个解$ P_0 $，那么对于任何常数$ c $来说$ P_0 + c $也是一个解，因为$ \ nabla \ cdot \ nabla（P_0 + c）= \ nabla \ cdot \ nabla P_0 = \ nabla \ cdot \ mathbf {v} ^ * $。
$ c $是我们要投影的空空间。
要处理Dirichlet BC，让我们以一维情况为例。假设您使用的是交错网格，其中压力$ p_i $位于网格中心，而速度$ v_ {i + 1/2} $位于节点$ i $和$ i + 1 $之间的面上。那么一个单元格的一般离散化为
$$
\ frac {p_ {i + 1} -p_i-（p_i-p_ {i-1}）} {\ Delta x ^ 2} = \文本{rhs}
$$
假设$ p_ {i + 1} $是一个气囊，也就是说，它的压力已指定， } $ term移到右侧，并从矩阵中消失。请注意，对角项$ p_i $的计数仍为2。这就是为什么我说您的2x2示例不包含Dirichlet BC。
无论使用Dirichlet还是Neumann BC，矩阵始终是对称正定的。这就是为什么作者说
Static object and empty cells don’t disrupt this structure.
In that case pressure and velocity terms can disappear from both sides
评论
目前尚不清楚为什么静态对象和空单元格应该允许删除行和列。您是将这些行和列设置为零还是完全删除它们以提供较小的矩阵？如果问题出在您猜不到的地方，那么如果您乐于分享，那将有助于查看代码。理想的是MCVE
嘿trichoplax。据我所知，一个全零行或全零的矩阵将是奇异的，因此我将其从矩阵中删除以形成一个较小的矩阵（以及它们在b向量中的对应条目）。
今晚在靠近源的计算机旁边，我将编辑MCVE。
我还怀疑我可能在代码的其他地方做出了错误的假设，但这仅与矩阵结构有关（以及是否为单数形式）。我唯一想到的是什么才算是“表面电池”与空气电池或液体电池。如果这是与空气室相邻的液体室，那么我应该如何处理其对应的列/行？