用病态稀疏矩阵求解 Newton-Raphson 步

计算科学 矩阵 线性求解器 稀疏矩阵
2021-11-27 09:00:34

我正在尝试为多孔介质中的质量和热量传输构建一个复杂的模拟器。我目前粗略地遵循旧软件中制定的算法并让模拟器工作，但我使用的迭代求解器在大时间步长上遇到问题，并阻止系统收敛到稳定状态。我想知道是否有人可以帮助我改进我的方法。

模拟器使用双精度数通过数值微分构建雅可比矩阵。所需的最终矩阵大小约为 500,000x500,000，每行大约有 15-20 个条目。雅可比的结构如下：

  • 对角线条目源于系统中元素的热力学性质。
  • 非对角线元素描述了元素之间的质量/热量传输。它们随时间步长缩放。由于目前使用的配方，它们是非对称的，但仅在元素之间的体积比上有所不同；如果需要，它们可以制成对称的。

系统 J * (-dx) = r 使用迭代求解器求解；特别是带有 ILU0 的 BICGSTAB 工作正常，但我也尝试过其他的。

当时间步长变大时，问题就出现了。对角线元素保持不变，但非对角线元素会随着时间步长缩放，并且一旦时间步长达到 ~1e+8 秒，使用的求解器似乎就会遇到问题。

最好的方法是什么？我正在考虑以下一项或多项：

  • 使用精度更高的求解器。但是，考虑到我的导数仅以约 8 位有效数字进行，这会起作用吗？
  • 分解矩阵是一种选择吗？如果是这样，什么分解会起作用？
  • 其他求解器是否可以选择？

非常感谢任何帮助。

干杯

彼得

2个回答

有关策略的一般列表，请参阅：

我在整个答案中提出了一堆问题；请在您的问题中将其作为编辑回复，因为这些问题的某些答案是有助于其他用户回答您的问题的信息。

最好的方法是什么？我正在考虑以下一项或多项：

使用精度更高的求解器。但是，考虑到我的导数仅以约 8 位有效数字进行，这会起作用吗？

我相信在求解器中使用更高精度的主要后果是在求解器中保留更多有效数字。如果有限差分是在求解器之外完成的，这将（显然）不受影响；如果以更高的精度计算导数，则它在选择数值差分参数方面为您提供了更多选项。您是否尝试过改变该差异参数并检查其对系统调节的影响？

使用解析雅可比矩阵（或其相当好的解析近似）是一种选择吗？当数值差分参数与问题中的相关数值参数（或比例）具有相同的顺序时，可能会出现问题。这个问题是我在调试反应扩散模拟时遇到的一个问题：雅可比矩阵的无矩阵近似不够准确，并导致了病态问题。替换雅可比矩阵的无矩阵近似解决了我的问题。（另请参阅：Newton-Krylov 何时不是合适的求解器？）

分解矩阵是一种选择吗？如果是这样，什么分解会起作用？

不，这不对。使用稀疏 LU 分解之类的东西可能需要太多内存（您已经有大约 1000 万个非零，而填充只会增加这个数字）。

其他求解器是否可以选择？

大概。你试过什么？对分析也有帮助的是对系统状况的估计。GMRES（没有重新启动）提供了一种这样的估计，了解您正在求解的线性系统是否病态会很有帮助。

具体来说，了解以下内容会很有帮助：

  • 时间步进器
  • 预处理器
  • 线性求解器
  • 图书馆

您正在使用。（提到 ILU(0) 和 BiCGStab 是一个好的开始。）

查看方程式并更多地了解您的问题也会有所帮助。是抛物线吗？在多个尺度上是否可能存在相同类型的错误？多重网格方法可能有助于预处理，并且比 ILU(0) 更好地扩展。

1e8 是一个非常大的时间步长，没有任何其他上下文。为什么会有这么大的数字？有规模问题吗？问题可以重新缩放或无量纲化吗？

分解矩阵是一种选择吗？如果是这样，什么分解会起作用？

是的，我认为这是一个非常合理的选择。500000 并不是大量的方程，而且您的系统非常稀疏。对于您当前的不对称情况，您需要执行分解。如果您可以按照您的建议对其进行对称化，则可以通过执行或 分解来节省近 1/2 的内存。LULLTLDLT

我建议尝试 MUMPS 稀疏求解器：

http://mumps.enseeiht.fr/

它具有所有三种分解类型的选项。它被设计为在分布式内存计算机上并行运行，但在单 CPU 或多核机器上也能非常高效地运行。它有一个核外选项，您可以使用它来减少内存使用量。

其它你可能感兴趣的问题
上一篇Python 中的自适应 ODE 算法 下一篇如何在 paraview 中进行几何清理？