公告:1)网站程序升级:Q2A升级到1.8.6,Wordpress升级到5.7.2
2)修复了头像加载慢与提交问题反应慢等问题
2021-06-16

欢迎来到 Mathematica 问答社区

提问时请贴上文本代码

语法高亮:在编辑器中点击

被禁止的话题:广告破解

请阅读:《提问的智慧》

备用域名:mma.ooo

支持LaTex数学公式
行内公式标识符:\$ 或“$\backslash ($”+“$\backslash )$”,
行间公式标识符:\$\$ 或 “$\backslash [$”+“$\backslash ]$”

社区建议QQ群:365716997

分类

0 投票
20 浏览

想了很久,还是觉得这个“问题”不太适合发在 stackexchange ,姑且在中文社区试试水吧。

大约八年前,我在这个网站(需爬墙)看到了一段有问题的卡门涡街模拟,代码如下:

(* Note: Something is wrong with this Mathematica code. 
   Please tell me if you find out what it is. 
   runtime: 14 seconds *)
n = 64; dt = 0.3; mu = 0.001; v = Table[{0, 0}, {n}, {n}];
Do[
Do[If[i < 5, v[[i, j]] = {0.1, 0}]; 
If[(i - n/4)^2 + (j -n/2)^2 < 4^2, v[[i, j]] = {0, 0}], {i, 1, n}, {j, 1, n}]; 
ui = ListInterpolation[v[[All, All,1]]]; 
vi = ListInterpolation[v[[All, All, 2]]]; 
v = Table[{i2, j2} = {i, j} - n dt v[[i, j]]; 
{ui[i2, j2], vi[i2, j2]}, {i, 1, n}, {j, 1, n}]; 
v = Transpose[Map[Fourier[v[[All, All, #]]] &, {1,2}], {3, 1, 2}]; 
v = Table[x = Mod[i + n/2, n] - n/2; y = Mod[j + n/2, n] - n/2; 
          k = x^2 + y^2; 
          Exp[-k dt mu]If[k > 0, (v[[i, j]].{-y, x}/k){-y, x}, v[[i, j]]], 
          {i, 1, n}, {j, 1, n}]; 
v = Transpose[Map[Re[InverseFourier[v[[All, All, #]]]] &, {1, 2}], {3, 1, 2}]; 
ListDensityPlot[
  Table[((v[[i + 1, j, 2]] - v[[i - 1, j, 2]]) - 
         (v[[i, j + 1, 1]] - v[[i, j - 1, 1]]))/2, 
    {j, 2, n - 1}, {i, 2, n - 1}], 
         Mesh -> False,Frame -> False, 
         ColorFunction -> (Hue[2#/3] &)], {t, 1, 25}];

如注释所说,这段代码有问题。在过去八年时间里,每隔一两年,我就会心血来潮地试图修复这段代码,然后以失败告终。大家有什么主意吗?

照例说说我的失败尝试。按网页上的说法,这段代码参考的是  Jos Stam 的论文《A Simple Fluid Solver based on the FFT》及其中给出的C代码。为方便阅读,论文及代码我已在本站传了一份:

论文《A Simple Fluid Solver based on the FFT》

solver.c

但这篇文章的细节并不十分清楚,我的C语言水平也很有限(我甚至搞不清楚文章中所用的 fftw 库要怎么在 windows 下“安装”),所以看完了照样不知道该怎么修。目前,我能看出来的只有2点:

1. 这段C代码的快速傅立叶变换基于的是 fftw-2.1.5 ,从相应的文档来看,该程序包所含的快速傅立叶变换大概相当于 Fourier[…, FourierParameters -> {1, -1}] ,但把 FourierParameters 加上似乎并不足以修复上述代码。

2. ListInterpolation 应该是得强制搞成周期性的,但这依旧不足以修复上述代码。

顺便,网页上附了一个模拟结果:

卡门涡街

据网页上的解说,此模拟用的是C++,但具体代码未给出。

 

用户: xzczd (2.2k 分)
修改于 用户:xzczd
这个给力了,不过看上去对我来说已经超纲了

登录 或者 注册 后回答这个问题。

...