在时间向前推进后使用 FFTW 和 NaN 的非三角函数、周期函数的精度损失 (Fortran)

计算科学正则纳维斯托克斯精确傅里叶变换 fftw

2021-12-16 00:24:41

我使用 FFTW 开发了 Navier-Stokes 方程的伪谱求解器。我根据标准三角函数（正弦、余弦及其组合）测试了我对 NS 方程右手边 (RHS) 的表述。例如，我设置

density = sin 5x
x_velocity = 5cos 5y + 6sin7z
y_velocity = 4sin4y + cos x
z_velocity = 1
pressure = cos z

将这些值提供给求解器，它计算了 NS 方程的 RHS。我手动做了同样的事情，并将结果与求解器获得的结果进行了比较。结果非常一致。准确答案与求解器计算的结果之间的最大误差E-13约为128*128*128 grid.

接下来，我使用了以下形式的不同功能：

density = constant1+constant2*(tanh(x-constant3)-tanh(x-constant4))
x_velocity = 0
y_velocity = 0
z_velocity = 0
temperature = 1
pressure -> from ideal gas equation connecting density, temperature and pressure

密度根据常数适当调整，周期为2*pi。根据给定的这些值计算 x 动量 Navier Stokes 的 RHS 并将其与我的答案（手动计算）进行比较，我得到了E-03.

此外，使用这些值作为变量的初始值并通过龙格-库塔 4 方案及时向前移动，我得到的密度值似乎很快就发散了。大约 30 个时间步后，我得到了 NaN。

当使用非三角周期函数时，我注意到精度下降有什么具体原因吗？
1.与为什么我的代码在及时前进时似乎会产生不稳定的结果有关吗？

我不介意在此处粘贴代码，但它非常大。

我想我会绘制初始密度及其变化。但事实证明我不能，因为我没有足够的声誉这样做。

初始图 (@t = 0.0s) 是一个密度图，看起来像一个矩形波，其中 tanh 函数用于平滑各个角处的波。

在 t = 0.10 秒左右（时间步长为 0.01 秒，因此，经过 10 次迭代），它会出现尖峰并变得不可微（仍然是连续的）。

1个回答

在伪光谱方法中，有三个问题可能会导致此类问题：

吉布斯振荡
别名
时间步长太大

在任何情况下，您都可能在解决方案中产生振荡，直到某个点以负密度结束，从而在计算压力或声速或其他项时产生 NaN。3 的解很明显，减小时间步长直到时间积分稳定。其他两个更微妙。

吉布斯振荡

计算不连续函数的傅里叶级数时会出现吉布斯振荡。如果函数不光滑，则在导数中会出现吉布斯振荡。如果您的初始条件有较大的跳跃，则傅立叶级数将与网格点处的值完全匹配，但导数会有很大的振荡，导致导数（右手边）计算的精度损失。请参见下图进行演示，值匹配但导数不匹配。根据经验，跳跃必须在大约 10 个网格上平滑，以防止这种行为。

吉布斯振荡的图像

即使您的初始条件在网格尺度上是平滑的，状态变量也可能会迅速变陡。在可压缩的 Navier-Stokes 中，粘性项的作用是防止形成冲击，但如果您的模拟没有得到充分解决，您仍然会在模拟中产生跳跃。充分解决意味着网格间距足够小以捕获粘性耗散，这可以通过查看 Kolmogorov 标度来估计，请参阅此PDF。这很快会导致大的、非物理的振荡和不同的解决方案。

别名

由于演化方程中存在非线性项（例如计算光谱空间中的导数假设您正在解析一定数量的波长。然而，非线性项不断产生越来越高的波数。在离散问题中，这些较高的波数被“混叠”回来以影响实际可以在所选分辨率下表示的较低波数。这会破坏较低的波数值，并可能迅速导致振荡、非物理结果和模拟崩溃。 $u u_x$

以下场景显示了非线性项如何生成更高波数以及这些波数如何混叠为更低波数的简单演示（见下图）：
在区间 [0,1) 上取 7 个点的网格。设和。这些项的（角）频率均为。项的频率为。该频率无法在所选网格上解析，但网格上的值等于网格上的值。所以不是 $a = \cos(6\pi x )$ $b = \sin(6\pi x)$ $6\pi$

a b = \cos (6 π x) \sin (6 π x) = \frac{1}{2} \sin (12 π x)

$ab = \cos(6\pi x )\sin(6\pi x ) = \frac{1}{2}\sin(12\pi x)$

12 π

$12\pi$

a b

$ab$

- \sin (2 π x) / 2

$-\sin(2\pi x)/2$

12 π

$12\pi$ 频率分量（未在离散空间中捕获），的结果显示为频率分量。

a b

$ab$

2 π

$2\pi$

演示锯齿的图像

有多个选项可用于防止混叠破坏您的结果，通常称为去混叠。最常见的方法是零填充（3/2 规则，在非线性乘法之前有效地增加网格分辨率，然后丢弃较高的频率），2/3 规则截断（在非线性乘法之前将最高波数归零，原始论文），或过滤程序（例如这种方法）。

此外，有证据表明，对于解析良好的模拟，去混叠并不重要，因为最高波数分量（可能产生混叠）由于粘性耗散而很小。

此演示文稿 (PDF)还提供了一个很好的去锯齿概述，包括一些历史。

其它你可能感兴趣的问题

上一篇如何有效计算网格边缘中点？下一篇与 GPU 相比，为什么 CPU 中的 SIMD 编程水平如此之低？