Matlab/CUDA: 海浪模拟
Matlab/CUDA: ocean wave simulation
我研究了 Jerry Tessendorf 的 "Simulating Ocean Water" 文章并尝试编写统计波模型，但我没有得到正确的结果，我不明白为什么。
在我的程序中，我只尝试在 t = 0 时间创建一个波高场，而没有任何进一步的时间变化。执行我的程序后，我得到的不是我所期望的：
这是我的源代码：
clear all; close all; clc;
rng(11); % setting seed for random numbers
meshSize = 64; % field size
windDir = [1, 0]; % ||windDir|| = 1
patchSize = 64;
A = 1e+4;
g = 9.81; % gravitational constant
windSpeed = 1e+2;
x1 = linspace(-10, 10, meshSize+1); x = x1(1:meshSize);
y1 = linspace(-10, 10, meshSize+1); y = y1(1:meshSize);
[X,Y] = meshgrid(x, y);
H0 = zeros(size(X)); % height field at time t = 0
for i = 1:meshSize
for j = 1:meshSize
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + x(i)); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + y(j)); % = 2*pi*m / Ly
P = phillips(kx, ky, windDir, windSpeed, A, g); % phillips spectrum
H0(i,j) = 1/sqrt(2) * (randn(1) + 1i * randn(1)) * sqrt(P);
end
end
H0 = H0 + conj(H0);
surf(X,Y,abs(ifft(H0)));
axis([-10 10 -10 10 -10 10]);
和phillips函数：
function P = phillips(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
end
有没有matlab海洋模拟源代码可以帮助我理解我的错误？快速 google 搜索没有得到任何结果。
这是我从 "CUDA FFT Ocean Simulation" 得到的 "correct" 结果。我还没有在 Matlab 中实现这种行为，但我已经使用来自 "CUDA FFT Ocean Simulation" 的数据在 Matlab 中绘制了 "surf"。这是它的样子：
我做了一个实验，得到了一个有趣的结果：
我从 "CUDA FFT Ocean Simulation" 中生成了 h0。所以我必须做 ifft 从频域转换到空间域来绘制图形。我已经使用 matlab ifft 并使用 CUDA 库中的 cufftExecC2C 完成了相同的 h0。结果如下：
CUDA ifft:
Matlab ifft:
要么我不明白cufftExecC2C的某些实现，要么cufftExecC2C和matlab ifft是不同的算法，结果不同。
顺便说一句，生成这种表面的参数是：
meshSize = 32
A = 1e-7
补丁大小=80
windSpeed = 10
那绝对是一个有趣的练习。这是一个完全重写的答案，因为你自己找到了你要问的问题。
与其删除我的答案，posting 仍然有帮助您向量化的优点and/or 解释一些代码。
我完全重写了我在之前的回答中提供的 GUI，以便合并您的更改并添加几个选项。它开始长出胳膊和腿，所以我不会把清单放在这里，但你可以在那里找到完整的文件：
这是完全独立的，它包括我矢量化并在下面单独列出的所有计算函数。
GUI 将允许您使用参数、动画波、导出 GIF 文件（以及一些其他选项，如 "preset"，但它们还没有完全解决）。您可以实现的一些示例：
基本
这是您通过快速默认设置和几个渲染选项获得的结果。这使用较小的网格大小和较快的时间步长，因此 运行 在任何机器上都非常快。
我在家里比较有限（Pentium E2200 32bit），所以我只能在有限的设置下练习。 gui 将 运行 即使设置已最大化，但真正享受它会变得很慢。
然而，在 Raid 中使用 ocean_simulator 的快速 运行（I7 64 位，8 核，16GB 内存，2xSSD），它让它变得更有趣！这里有几个例子：
虽然在更好的机器上完成，但我没有使用任何并行功能或任何 GPU 计算，所以 Matlab 只使用了这些规范的一部分，这意味着它可能 运行 在任何具有良好 RAM 的 64 位系统上都一样好
风湖
这是一个像湖一样平坦的水面。即使大风也不会产生高振幅波（但仍然有很多小波）。如果你是一名风帆冲浪者，在山顶上从你的 window 看，你的心会跳动，你的下一步是打电话给戴夫“伙计！准备好。 五点后在水面上见！"
膨胀
这是你在与暴风雨搏斗了一夜之后，早上从你的船桥上看到的。风暴已经消散，长长的巨浪是这绝对是一个摇摇欲坠的夜晚的最后见证（有航海经验的人都知道......）。
T-Storm
这就是你前一天晚上所做的...
第二张gif是在家做的，所以细节不够...抱歉
到底部：
最后，gui 会让你在水域周围添加一个补丁。在 gui 中，它是透明的，因此您可以在水下或漂亮的海底添加 objects。不幸的是，GIF 格式不能包含 alpha 通道，因此这里没有透明度（但如果您导出视频则应该没问题）。
此外，导出为 GIF 会降低图像质量，如果您在 Matlab 中 运行 域边界和水面之间的接合是完美的。在某些情况下，它还会使 Matlab 降低照明的渲染，因此这绝对不是导出的最佳选择，但它允许在 matlab 中发挥更多的作用。
现在进入代码：
而不是列出完整的 GUI，这将是超长的（这个 post 已经足够长了），我将在这里列出 re-written 版本的 你的 代码，并解释变化。
由于剩余的矢量化，您应该注意到执行速度大幅提高（数量级），但主要有两个原因：
(i) 重复了很多计算。缓存值并重新使用它们比在循环中（在动画部分期间）重新计算完整矩阵要快得多。
(ii) 注意我是如何定义表面图形的object。它只定义一次（即使是空的），然后所有进一步的调用（在循环中）只更新表面 object 的底层 ZData （而不是 re-creating 表面 object 在每次迭代中。
这里是：
%% // clear workspace
clear all; close all; clc;
%% // Default parameters
param.meshsize = 128 ; %// main grid size
param.patchsize = 200 ;
param.windSpeed = 100 ; %// what unit ? [m/s] ??
param.winddir = 90 ; %// Azimuth
param.rng = 13 ; %// setting seed for random numbers
param.A = 1e-7 ; %// Scaling factor
param.g = 9.81 ; %// gravitational constant
param.xLim = [-10 10] ; %// domain limits X
param.yLim = [-10 10] ; %// domain limits Y
param.zLim = [-1e-4 1e-4]*2 ;
gridSize = param.meshsize * [1 1] ;
%% // Define the grid X-Y domain
x = linspace( param.xLim(1) , param.xLim(2) , param.meshsize ) ;
y = linspace( param.yLim(1) , param.yLim(2) , param.meshsize ) ;
[X,Y] = meshgrid(x, y);
%% // get the grid parameters which remain constants (not time dependent)
[H0, W, Grid_Sign] = initialize_wave( param ) ;
%% // calculate wave at t0
t0 = 0 ;
Z = calc_wave( H0 , W , t0 , Grid_Sign ) ;
%% // populate the display panel
h.fig = figure('Color','w') ;
h.ax = handle(axes) ; %// create an empty axes that fills the figure
h.surf = handle( surf( NaN(2) ) ) ; %// create an empty "surface" object
%% // Display the initial wave surface
set( h.surf , 'XData',X , 'YData',Y , 'ZData',Z )
set( h.ax , 'XLim',param.xLim , 'YLim',param.yLim , 'ZLim',param.zLim )
%% // Change some rendering options
axis off %// make the axis grid and border invisible
shading interp %// improve shading (remove "faceted" effect)
blue = linspace(0.4, 1.0, 25).' ; cmap = [blue*0, blue*0, blue]; %'// create blue colormap
colormap(cmap)
%// configure lighting
h.light_handle = lightangle(-45,30) ; %// add a light source
set(h.surf,'FaceLighting','phong','AmbientStrength',.3,'DiffuseStrength',.8,'SpecularStrength',.9,'SpecularExponent',25,'BackFaceLighting','unlit')
%% // Animate
view(75,55) %// no need to reset the view inside the loop ;)
timeStep = 1./25 ;
nSteps = 2000 ;
for time = (1:nSteps)*timeStep
%// update wave surface
Z = calc_wave( H0,W,time,Grid_Sign ) ;
h.surf.ZData = Z ;
pause(0.001);
end
%% // This block of code is only if you want to generate a GIF file
%// be carefull on how many frames you put there, the size of the GIF can
%// quickly grow out of proportion ;)
nFrame = 55 ;
gifFileName = 'MyDancingWaves.gif' ;
view(-70,40)
clear im
f = getframe;
[im,map] = rgb2ind(f.cdata,256,'nodither');
im(1,1,1,20) = 0;
iframe = 0 ;
for time = (1:nFrame)*.5
%// update wave surface
Z = calc_wave( H0,W,time,Grid_Sign ) ;
h.surf.ZData = Z ;
pause(0.001);
f = getframe;
iframe= iframe+1 ;
im(:,:,1,iframe) = rgb2ind(f.cdata,map,'nodither');
end
imwrite(im,map,gifFileName,'DelayTime',0,'LoopCount',inf)
disp([num2str(nFrame) ' frames written in file: ' gifFileName])
您会注意到我更改了一些内容，但我可以向您保证计算结果完全相同。这段代码调用了一些子函数，但它们都是矢量化的，所以如果你愿意，你可以在这里 copy/paste 它们，运行 一切都是内联的。
调用的第一个函数是initialize_wave.m
此处计算的所有内容以后都会保持不变（当您稍后为波浪设置动画时，它不会随时间变化），因此将其单独放入一个块中是有意义的。
function [H0, W, Grid_Sign] = initialize_wave( param )
% function [H0, W, Grid_Sign] = initialize_wave( param )
%
% This function return the wave height coefficients H0 and W for the
% parameters given in input. These coefficients are constants for a given
% set of input parameters.
% Third output parameter is optional (easy to recalculate anyway)
rng(param.rng); %// setting seed for random numbers
gridSize = param.meshsize * [1 1] ;
meshLim = pi * param.meshsize / param.patchsize ;
N = linspace(-meshLim , meshLim , param.meshsize ) ;
M = linspace(-meshLim , meshLim , param.meshsize ) ;
[Kx,Ky] = meshgrid(N,M) ;
K = sqrt(Kx.^2 + Ky.^2); %// ||K||
W = sqrt(K .* param.g); %// deep water frequencies (empirical parameter)
[windx , windy] = pol2cart( deg2rad(param.winddir) , 1) ;
P = phillips(Kx, Ky, [windx , windy], param.windSpeed, param.A, param.g) ;
H0 = 1/sqrt(2) .* (randn(gridSize) + 1i .* randn(gridSize)) .* sqrt(P); % height field at time t = 0
if nargout == 3
Grid_Sign = signGrid( param.meshsize ) ;
end
请注意，初始 winDir 参数现在用表示 "azimuth"（ 度的单个标量值表示) 的风（从 0 到 360 的任何值）。由于函数 pol2cart，它后来被翻译成它的 X 和 Y 组件。
[windx , windy] = pol2cart( deg2rad(param.winddir) , 1) ;
这确保范数始终为 1。
该函数单独调用有问题的 phillips.m，但如前所述，它甚至可以完全矢量化工作，因此您可以根据需要将其内联复制回去。 （别担心，我根据您的版本检查了结果 => 完全相同）。请注意，此函数不输出复数，因此无需比较虚部。
function P = phillips(Kx, Ky, windDir, windSpeed, A, g)
%// The function now accept scalar, vector or full 2D grid matrix as input
K_sq = Kx.^2 + Ky.^2;
L = windSpeed.^2 ./ g;
k_norm = sqrt(K_sq) ;
WK = Kx./k_norm * windDir(1) + Ky./k_norm * windDir(2);
P = A ./ K_sq.^2 .* exp(-1.0 ./ (K_sq * L^2)) .* WK.^2 ;
P( K_sq==0 | WK<0 ) = 0 ;
end
主程序调用的下一个函数是calc_wave.m。该函数完成给定时间的波场计算。单独使用它绝对值得，因为这是一组最简单的计算，当您想要为波浪设置动画时，必须在每个给定时间重复这些计算。
function Z = calc_wave( H0,W,time,Grid_Sign )
% Z = calc_wave( H0,W,time,Grid_Sign )
%
% This function calculate the wave height based on the wave coefficients H0
% and W, for a given "time". Default time=0 if not supplied.
% Fourth output parameter is optional (easy to recalculate anyway)
% recalculate the grid sign if not supplied in input
if nargin < 4
Grid_Sign = signGrid( param.meshsize ) ;
end
% Assign time=0 if not specified in input
if nargin < 3 ; time = 0 ; end
wt = exp(1i .* W .* time ) ;
Ht = H0 .* wt + conj(rot90(H0,2)) .* conj(wt) ;
Z = real( ifft2(Ht) .* Grid_Sign ) ;
end
最后 3 行计算需要一些解释，因为它们收到了最大的变化（所有结果相同但速度更快）。
你原来的线路：
Ht = H0 .* exp(1i .* W .* (t * timeStep)) + conj(flip(flip(H0,1),2)) .* exp(-1i .* W .* (t * timeStep));
同一件事重新计算太多次效率不高：
(t * timeStep) 在每一行上计算两次，在每个循环中，而当 time 在开始时初始化时，很容易为每一行获得正确的 time 值循环 for time = (1:nSteps)*timeStep.
另请注意 exp(-1i .* W .* time) 与 conj(exp(1i .* W .* time)) 相同。与其进行 2*m*n 次乘法来分别计算它们，不如先计算一次，然后使用更快的 conj() 运算。
所以你的单行会变成：
wt = exp(1i .* W .* time ) ;
Ht = H0 .* wt + conj(flip(flip(H0,1),2)) .* conj(wt) ;
最后一点点，flip(flip(H0,1),2)) 可以替换为 rot90(H0,2)（也稍微快一点）。
请注意，由于函数 calc_wave 将被广泛重复，因此绝对值得减少计算次数（如我们上面所做的），但也可以通过向其发送 Grid_Sign 参数（而不是让函数在每次迭代时重新计算它）。这就是为什么：
你的神秘函数signCor(ifft2(Ht),meshSize))，简单地反转Ht的所有其他元素的符号。有一种更快的方法可以实现这一点：只需将 Ht 乘以一个相同大小的矩阵 (Grid_Sign)，这是一个交替 +1 -1 ... 的矩阵，依此类推。
所以 signCor(ifft2(Ht),meshSize) 变成 ifft2(Ht) .* Grid_Sign.
由于Grid_Sign只依赖于矩阵的大小，它不会随着循环中的每个time而改变，你只需计算一次（在循环之前）然后按原样使用它每隔一次迭代。它的计算方式如下（矢量化，因此您也可以将其内联到您的代码中）：
function sgn = signGrid(n)
% return a matrix the size of n with alternate sign for every indice
% ex: sgn = signGrid(3) ;
% sgn =
% -1 1 -1
% 1 -1 1
% -1 1 -1
[x,y] = meshgrid(1:n,1:n) ;
sgn = ones( n ) ;
sgn(mod(x+y,2)==0) = -1 ;
end
最后，您会注意到您的版本与此版本之间网格 [Kx,Ky] 的定义方式有所不同。它们确实会产生略有不同的结果，这只是一个选择问题。
为了用一个简单的例子来解释，让我们考虑一个小的meshsize=5。你的做事方式会将其分成 5 个值，等距，如下所示：
Kx(first line)=[-1.5 -0.5 0.5 1.5 2.5] * 2 * pi / patchSize
虽然我生成网格的方式会生成等距值，但也会以域限制为中心，如下所示：
Kx(first line)=[-2.50 -1.25 0.0 1.25 2.50] * 2 * pi / patchSize
似乎更尊重你在定义它的那一行的评论% = 2*pi*n / Lx, -N/2 <= n < N/2。
我倾向于偏向于对称的解决方案（加上它也稍微快一点但是只计算一次所以没什么大不了的），所以我使用了我的矢量化方式，但这纯粹是一个选择问题，你绝对可以按照您的方式行事，它只是稍微 "offset" 整个结果矩阵，但它不会扰乱计算本身。
第一个答案的最后遗留
侧面编程说明：
我检测到你来自 C/C++ 世界或家庭。在 Matlab 中，您不需要用逗号定义十进制数（例如 2.0，您将其用于大多数数字）。除非另有特别定义，否则 Matlab 默认将任何数字转换为 double，这是一个 64 位浮点类型。所以写 2 * pi 足以获得最大精度（Matlab 不会将 pi 转换为整数 ;-)），你不需要写 2.0 * pi . 虽然如果你不想改变你的习惯，它仍然有效。
此外，（Matlab 的一大好处），在运算符前添加 . 通常意味着 "element-wise" 运算。您可以通过这种方式明智地添加 (.+)、减去 (.-)、乘法 (.*)、除法 (./) 全矩阵元素。这就是我摆脱代码中所有循环的方法。这也适用于幂运算符：A.^2 将 return 一个与 A 大小相同的矩阵，每个元素都是平方的。
Here 是工作程序。
首先-源代码：
clear all; close all; clc;
rng(13); % setting seed for random numbers
meshSize = 128; % field size
windDir = [0.1,1];
patchSize = 200;
A = 1e-7;
g = 9.81; % gravitational constant
windSpeed = 100;
timeStep = 1/25;
x1 = linspace(-10, 10, meshSize+1); x = x1(1:meshSize);
y1 = linspace(-10, 10, meshSize+1); y = y1(1:meshSize);
[X,Y] = meshgrid(x,y); % wave field
i = 1:meshSize; j = 1:meshSize; % indecies
[I,J] = meshgrid(i,j); % field of indecies
Kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + I); % = 2*pi*n / Lx, -N/2 <= n < N/2
Ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + J); % = 2*pi*m / Ly, -M/2 <= m < M/2
K = sqrt(Kx.^2 + Ky.^2); % ||K||
W = sqrt(K .* g); % deep water frequencies (empirical parameter)
P = zeros(size(X)); % Cant compute P without loops
for i = 1:meshSize
for j = 1:meshSize
P(i,j) = phillips(Kx(i,j), Ky(i,j), windDir, windSpeed, A, g); % phillips spectrum
end
end
H0 = 1/sqrt(2) .* (randn(size(X)) + 1i .* randn(size(X))) .* sqrt(P); % height field at time t = 0
rotate3d on;
for t = 1:10000 % 10000 * timeStep (sec)
Ht = H0 .* exp(1i .* W .* (t * timeStep)) + ...
conj(flip(flip(H0,1),2)) .* exp(-1i .* W .* (t * timeStep));
[az,el] = view;
surf(X,Y,real(signCor(ifft2(Ht),meshSize)));
axis([-10 10 -10 10 -1e-4 1e-4]); view(az,el);
blue = linspace(0.4, 1.0, 25)'; map = [blue*0, blue*0, blue];
%shading interp; % improve shading (remove "faceted" effect)
colormap(map);
pause(1/60);
end
phillips.m：（我试图向量化菲利普斯谱的计算，但我遇到了一个困难，我将进一步展示）
function P = phillips(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
if k_sq == 0
P = 0;
else
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
if wk < 0
P = 0;
end
end
end
signCor.m：（这个函数对我来说绝对是个谜...我从 "CUDA FFT Ocean Simulation" 实现中复制了它。没有它，模拟效果会更差。而且我也不知道如何向量化这个函数。）
function H = signCor(H1, meshSize)
H = H1;
for i = 1:meshSize
for j = 1:meshSize
if mod(i+j,2) == 0
sign = -1; % works fine if we change signs vice versa
else
sign = 1;
end
H(i,j) = H1(i,j) * sign;
end
end
end
我犯的最大错误是我使用了 ifft 而不是 ifft2，这就是为什么 CUDA ifft 和 Matlab ifft 不匹配的原因。
我的第二个错误是在这行代码中：
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + x(i)); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + y(j)); % = 2*pi*m / Ly
我应该写：
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + i); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + j); % = 2*pi*m / Ly
我对参数 A、meshSize、patchSize 进行了一些研究，得出的结论是：
波幅的合理参数是 A * (patchSize / meshSize)，其中 A 只是一个比例因子。
为'calm' patchSize / meshSize <= 0.5.
为'tsunami' patchSize / meshSize >= 3.0.
菲利普斯谱矢量化困难：
我有两个功能：
% non-vectorized spectrum
function P = phillips1(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
if k_sq == 0
P = 0;
else
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
if wk < 0
P = 0;
end
end
end
% vectorized spectrum
function P = phillips2(Kx, Ky, windDir, windSpeed, A, g)
K_sq = Kx .^ 2 + Ky .^ 2;
L = -g^2 / windSpeed^4;
WK = (Kx ./ K_sq) .* windDir(1) + (Ky ./ K_sq) .* windDir(2);
P = (A ./ (K_sq .^ 2)) .* ( exp(L ./ K_sq) .* (WK .^ 2) );
P(K_sq == 0) = 0;
P(WK < 0) = 0;
P(isinf(P)) = 0;
end
在我使用 phillips1 计算 P1 和使用 phillips2 计算 P2 之后，我绘制了它们的差异：
subplot(2,1,1); surf(X,Y,real(P2-P1)); title('Difference in real part');
subplot(2,1,2); surf(X,Y,imag(P2-P1)); title('Difference in imaginary part');
完美地说明了这2个光谱在实部之间存在着巨大的差异。
我研究了 Jerry Tessendorf 的 "Simulating Ocean Water" 文章并尝试编写统计波模型，但我没有得到正确的结果，我不明白为什么。
在我的程序中，我只尝试在 t = 0 时间创建一个波高场，而没有任何进一步的时间变化。执行我的程序后，我得到的不是我所期望的：
这是我的源代码：
clear all; close all; clc;
rng(11); % setting seed for random numbers
meshSize = 64; % field size
windDir = [1, 0]; % ||windDir|| = 1
patchSize = 64;
A = 1e+4;
g = 9.81; % gravitational constant
windSpeed = 1e+2;
x1 = linspace(-10, 10, meshSize+1); x = x1(1:meshSize);
y1 = linspace(-10, 10, meshSize+1); y = y1(1:meshSize);
[X,Y] = meshgrid(x, y);
H0 = zeros(size(X)); % height field at time t = 0
for i = 1:meshSize
for j = 1:meshSize
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + x(i)); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + y(j)); % = 2*pi*m / Ly
P = phillips(kx, ky, windDir, windSpeed, A, g); % phillips spectrum
H0(i,j) = 1/sqrt(2) * (randn(1) + 1i * randn(1)) * sqrt(P);
end
end
H0 = H0 + conj(H0);
surf(X,Y,abs(ifft(H0)));
axis([-10 10 -10 10 -10 10]);
和phillips函数：
function P = phillips(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
end
有没有matlab海洋模拟源代码可以帮助我理解我的错误？快速 google 搜索没有得到任何结果。
这是我从 "CUDA FFT Ocean Simulation" 得到的 "correct" 结果。我还没有在 Matlab 中实现这种行为，但我已经使用来自 "CUDA FFT Ocean Simulation" 的数据在 Matlab 中绘制了 "surf"。这是它的样子：
我做了一个实验，得到了一个有趣的结果：
我从 "CUDA FFT Ocean Simulation" 中生成了 h0。所以我必须做 ifft 从频域转换到空间域来绘制图形。我已经使用 matlab ifft 并使用 CUDA 库中的 cufftExecC2C 完成了相同的 h0。结果如下：
CUDA ifft:
Matlab ifft:
要么我不明白cufftExecC2C的某些实现，要么cufftExecC2C和matlab ifft是不同的算法，结果不同。
顺便说一句，生成这种表面的参数是：
meshSize = 32
A = 1e-7
补丁大小=80
windSpeed = 10
那绝对是一个有趣的练习。这是一个完全重写的答案，因为你自己找到了你要问的问题。
与其删除我的答案，posting 仍然有帮助您向量化的优点and/or 解释一些代码。
我完全重写了我在之前的回答中提供的 GUI，以便合并您的更改并添加几个选项。它开始长出胳膊和腿，所以我不会把清单放在这里，但你可以在那里找到完整的文件：
这是完全独立的，它包括我矢量化并在下面单独列出的所有计算函数。
GUI 将允许您使用参数、动画波、导出 GIF 文件（以及一些其他选项，如 "preset"，但它们还没有完全解决）。您可以实现的一些示例：
基本
这是您通过快速默认设置和几个渲染选项获得的结果。这使用较小的网格大小和较快的时间步长，因此 运行 在任何机器上都非常快。
我在家里比较有限（Pentium E2200 32bit），所以我只能在有限的设置下练习。 gui 将 运行 即使设置已最大化，但真正享受它会变得很慢。
然而，在 Raid 中使用 ocean_simulator 的快速 运行（I7 64 位，8 核，16GB 内存，2xSSD），它让它变得更有趣！这里有几个例子：
虽然在更好的机器上完成，但我没有使用任何并行功能或任何 GPU 计算，所以 Matlab 只使用了这些规范的一部分，这意味着它可能 运行 在任何具有良好 RAM 的 64 位系统上都一样好
风湖
这是一个像湖一样平坦的水面。即使大风也不会产生高振幅波（但仍然有很多小波）。如果你是一名风帆冲浪者，在山顶上从你的 window 看，你的心会跳动，你的下一步是打电话给戴夫“伙计！准备好。 五点后在水面上见！"
膨胀
这是你在与暴风雨搏斗了一夜之后，早上从你的船桥上看到的。风暴已经消散，长长的巨浪是这绝对是一个摇摇欲坠的夜晚的最后见证（有航海经验的人都知道......）。
T-Storm
这就是你前一天晚上所做的...
到底部：
最后，gui 会让你在水域周围添加一个补丁。在 gui 中，它是透明的，因此您可以在水下或漂亮的海底添加 objects。不幸的是，GIF 格式不能包含 alpha 通道，因此这里没有透明度（但如果您导出视频则应该没问题）。
此外，导出为 GIF 会降低图像质量，如果您在 Matlab 中 运行 域边界和水面之间的接合是完美的。在某些情况下，它还会使 Matlab 降低照明的渲染，因此这绝对不是导出的最佳选择，但它允许在 matlab 中发挥更多的作用。
现在进入代码：
而不是列出完整的 GUI，这将是超长的（这个 post 已经足够长了），我将在这里列出 re-written 版本的 你的 代码，并解释变化。
由于剩余的矢量化，您应该注意到执行速度大幅提高（数量级），但主要有两个原因：
(i) 重复了很多计算。缓存值并重新使用它们比在循环中（在动画部分期间）重新计算完整矩阵要快得多。
(ii) 注意我是如何定义表面图形的object。它只定义一次（即使是空的），然后所有进一步的调用（在循环中）只更新表面 object 的底层 ZData （而不是 re-creating 表面 object 在每次迭代中。
这里是：
%% // clear workspace
clear all; close all; clc;
%% // Default parameters
param.meshsize = 128 ; %// main grid size
param.patchsize = 200 ;
param.windSpeed = 100 ; %// what unit ? [m/s] ??
param.winddir = 90 ; %// Azimuth
param.rng = 13 ; %// setting seed for random numbers
param.A = 1e-7 ; %// Scaling factor
param.g = 9.81 ; %// gravitational constant
param.xLim = [-10 10] ; %// domain limits X
param.yLim = [-10 10] ; %// domain limits Y
param.zLim = [-1e-4 1e-4]*2 ;
gridSize = param.meshsize * [1 1] ;
%% // Define the grid X-Y domain
x = linspace( param.xLim(1) , param.xLim(2) , param.meshsize ) ;
y = linspace( param.yLim(1) , param.yLim(2) , param.meshsize ) ;
[X,Y] = meshgrid(x, y);
%% // get the grid parameters which remain constants (not time dependent)
[H0, W, Grid_Sign] = initialize_wave( param ) ;
%% // calculate wave at t0
t0 = 0 ;
Z = calc_wave( H0 , W , t0 , Grid_Sign ) ;
%% // populate the display panel
h.fig = figure('Color','w') ;
h.ax = handle(axes) ; %// create an empty axes that fills the figure
h.surf = handle( surf( NaN(2) ) ) ; %// create an empty "surface" object
%% // Display the initial wave surface
set( h.surf , 'XData',X , 'YData',Y , 'ZData',Z )
set( h.ax , 'XLim',param.xLim , 'YLim',param.yLim , 'ZLim',param.zLim )
%% // Change some rendering options
axis off %// make the axis grid and border invisible
shading interp %// improve shading (remove "faceted" effect)
blue = linspace(0.4, 1.0, 25).' ; cmap = [blue*0, blue*0, blue]; %'// create blue colormap
colormap(cmap)
%// configure lighting
h.light_handle = lightangle(-45,30) ; %// add a light source
set(h.surf,'FaceLighting','phong','AmbientStrength',.3,'DiffuseStrength',.8,'SpecularStrength',.9,'SpecularExponent',25,'BackFaceLighting','unlit')
%% // Animate
view(75,55) %// no need to reset the view inside the loop ;)
timeStep = 1./25 ;
nSteps = 2000 ;
for time = (1:nSteps)*timeStep
%// update wave surface
Z = calc_wave( H0,W,time,Grid_Sign ) ;
h.surf.ZData = Z ;
pause(0.001);
end
%% // This block of code is only if you want to generate a GIF file
%// be carefull on how many frames you put there, the size of the GIF can
%// quickly grow out of proportion ;)
nFrame = 55 ;
gifFileName = 'MyDancingWaves.gif' ;
view(-70,40)
clear im
f = getframe;
[im,map] = rgb2ind(f.cdata,256,'nodither');
im(1,1,1,20) = 0;
iframe = 0 ;
for time = (1:nFrame)*.5
%// update wave surface
Z = calc_wave( H0,W,time,Grid_Sign ) ;
h.surf.ZData = Z ;
pause(0.001);
f = getframe;
iframe= iframe+1 ;
im(:,:,1,iframe) = rgb2ind(f.cdata,map,'nodither');
end
imwrite(im,map,gifFileName,'DelayTime',0,'LoopCount',inf)
disp([num2str(nFrame) ' frames written in file: ' gifFileName])
您会注意到我更改了一些内容，但我可以向您保证计算结果完全相同。这段代码调用了一些子函数，但它们都是矢量化的，所以如果你愿意，你可以在这里 copy/paste 它们，运行 一切都是内联的。
调用的第一个函数是initialize_wave.m
此处计算的所有内容以后都会保持不变（当您稍后为波浪设置动画时，它不会随时间变化），因此将其单独放入一个块中是有意义的。
function [H0, W, Grid_Sign] = initialize_wave( param )
% function [H0, W, Grid_Sign] = initialize_wave( param )
%
% This function return the wave height coefficients H0 and W for the
% parameters given in input. These coefficients are constants for a given
% set of input parameters.
% Third output parameter is optional (easy to recalculate anyway)
rng(param.rng); %// setting seed for random numbers
gridSize = param.meshsize * [1 1] ;
meshLim = pi * param.meshsize / param.patchsize ;
N = linspace(-meshLim , meshLim , param.meshsize ) ;
M = linspace(-meshLim , meshLim , param.meshsize ) ;
[Kx,Ky] = meshgrid(N,M) ;
K = sqrt(Kx.^2 + Ky.^2); %// ||K||
W = sqrt(K .* param.g); %// deep water frequencies (empirical parameter)
[windx , windy] = pol2cart( deg2rad(param.winddir) , 1) ;
P = phillips(Kx, Ky, [windx , windy], param.windSpeed, param.A, param.g) ;
H0 = 1/sqrt(2) .* (randn(gridSize) + 1i .* randn(gridSize)) .* sqrt(P); % height field at time t = 0
if nargout == 3
Grid_Sign = signGrid( param.meshsize ) ;
end
请注意，初始 winDir 参数现在用表示 "azimuth"（ 度的单个标量值表示) 的风（从 0 到 360 的任何值）。由于函数 pol2cart，它后来被翻译成它的 X 和 Y 组件。
[windx , windy] = pol2cart( deg2rad(param.winddir) , 1) ;
这确保范数始终为 1。
该函数单独调用有问题的 phillips.m，但如前所述，它甚至可以完全矢量化工作，因此您可以根据需要将其内联复制回去。 （别担心，我根据您的版本检查了结果 => 完全相同）。请注意，此函数不输出复数，因此无需比较虚部。
function P = phillips(Kx, Ky, windDir, windSpeed, A, g)
%// The function now accept scalar, vector or full 2D grid matrix as input
K_sq = Kx.^2 + Ky.^2;
L = windSpeed.^2 ./ g;
k_norm = sqrt(K_sq) ;
WK = Kx./k_norm * windDir(1) + Ky./k_norm * windDir(2);
P = A ./ K_sq.^2 .* exp(-1.0 ./ (K_sq * L^2)) .* WK.^2 ;
P( K_sq==0 | WK<0 ) = 0 ;
end
主程序调用的下一个函数是calc_wave.m。该函数完成给定时间的波场计算。单独使用它绝对值得，因为这是一组最简单的计算，当您想要为波浪设置动画时，必须在每个给定时间重复这些计算。
function Z = calc_wave( H0,W,time,Grid_Sign )
% Z = calc_wave( H0,W,time,Grid_Sign )
%
% This function calculate the wave height based on the wave coefficients H0
% and W, for a given "time". Default time=0 if not supplied.
% Fourth output parameter is optional (easy to recalculate anyway)
% recalculate the grid sign if not supplied in input
if nargin < 4
Grid_Sign = signGrid( param.meshsize ) ;
end
% Assign time=0 if not specified in input
if nargin < 3 ; time = 0 ; end
wt = exp(1i .* W .* time ) ;
Ht = H0 .* wt + conj(rot90(H0,2)) .* conj(wt) ;
Z = real( ifft2(Ht) .* Grid_Sign ) ;
end
最后 3 行计算需要一些解释，因为它们收到了最大的变化（所有结果相同但速度更快）。
你原来的线路：
Ht = H0 .* exp(1i .* W .* (t * timeStep)) + conj(flip(flip(H0,1),2)) .* exp(-1i .* W .* (t * timeStep));
同一件事重新计算太多次效率不高：
(t * timeStep) 在每一行上计算两次，在每个循环中，而当 time 在开始时初始化时，很容易为每一行获得正确的 time 值循环 for time = (1:nSteps)*timeStep.
另请注意 exp(-1i .* W .* time) 与 conj(exp(1i .* W .* time)) 相同。与其进行 2*m*n 次乘法来分别计算它们，不如先计算一次，然后使用更快的 conj() 运算。
所以你的单行会变成：
wt = exp(1i .* W .* time ) ;
Ht = H0 .* wt + conj(flip(flip(H0,1),2)) .* conj(wt) ;
最后一点点，flip(flip(H0,1),2)) 可以替换为 rot90(H0,2)（也稍微快一点）。
请注意，由于函数 calc_wave 将被广泛重复，因此绝对值得减少计算次数（如我们上面所做的），但也可以通过向其发送 Grid_Sign 参数（而不是让函数在每次迭代时重新计算它）。这就是为什么：
你的神秘函数signCor(ifft2(Ht),meshSize))，简单地反转Ht的所有其他元素的符号。有一种更快的方法可以实现这一点：只需将 Ht 乘以一个相同大小的矩阵 (Grid_Sign)，这是一个交替 +1 -1 ... 的矩阵，依此类推。
所以 signCor(ifft2(Ht),meshSize) 变成 ifft2(Ht) .* Grid_Sign.
由于Grid_Sign只依赖于矩阵的大小，它不会随着循环中的每个time而改变，你只需计算一次（在循环之前）然后按原样使用它每隔一次迭代。它的计算方式如下（矢量化，因此您也可以将其内联到您的代码中）：
function sgn = signGrid(n)
% return a matrix the size of n with alternate sign for every indice
% ex: sgn = signGrid(3) ;
% sgn =
% -1 1 -1
% 1 -1 1
% -1 1 -1
[x,y] = meshgrid(1:n,1:n) ;
sgn = ones( n ) ;
sgn(mod(x+y,2)==0) = -1 ;
end
最后，您会注意到您的版本与此版本之间网格 [Kx,Ky] 的定义方式有所不同。它们确实会产生略有不同的结果，这只是一个选择问题。
为了用一个简单的例子来解释，让我们考虑一个小的meshsize=5。你的做事方式会将其分成 5 个值，等距，如下所示：
Kx(first line)=[-1.5 -0.5 0.5 1.5 2.5] * 2 * pi / patchSize
虽然我生成网格的方式会生成等距值，但也会以域限制为中心，如下所示：
Kx(first line)=[-2.50 -1.25 0.0 1.25 2.50] * 2 * pi / patchSize
似乎更尊重你在定义它的那一行的评论% = 2*pi*n / Lx, -N/2 <= n < N/2。
我倾向于偏向于对称的解决方案（加上它也稍微快一点但是只计算一次所以没什么大不了的），所以我使用了我的矢量化方式，但这纯粹是一个选择问题，你绝对可以按照您的方式行事，它只是稍微 "offset" 整个结果矩阵，但它不会扰乱计算本身。
第一个答案的最后遗留
侧面编程说明：
我检测到你来自 C/C++ 世界或家庭。在 Matlab 中，您不需要用逗号定义十进制数（例如 2.0，您将其用于大多数数字）。除非另有特别定义，否则 Matlab 默认将任何数字转换为 double，这是一个 64 位浮点类型。所以写 2 * pi 足以获得最大精度（Matlab 不会将 pi 转换为整数 ;-)），你不需要写 2.0 * pi . 虽然如果你不想改变你的习惯，它仍然有效。
此外，（Matlab 的一大好处），在运算符前添加 . 通常意味着 "element-wise" 运算。您可以通过这种方式明智地添加 (.+)、减去 (.-)、乘法 (.*)、除法 (./) 全矩阵元素。这就是我摆脱代码中所有循环的方法。这也适用于幂运算符：A.^2 将 return 一个与 A 大小相同的矩阵，每个元素都是平方的。
Here 是工作程序。
首先-源代码：
clear all; close all; clc;
rng(13); % setting seed for random numbers
meshSize = 128; % field size
windDir = [0.1,1];
patchSize = 200;
A = 1e-7;
g = 9.81; % gravitational constant
windSpeed = 100;
timeStep = 1/25;
x1 = linspace(-10, 10, meshSize+1); x = x1(1:meshSize);
y1 = linspace(-10, 10, meshSize+1); y = y1(1:meshSize);
[X,Y] = meshgrid(x,y); % wave field
i = 1:meshSize; j = 1:meshSize; % indecies
[I,J] = meshgrid(i,j); % field of indecies
Kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + I); % = 2*pi*n / Lx, -N/2 <= n < N/2
Ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + J); % = 2*pi*m / Ly, -M/2 <= m < M/2
K = sqrt(Kx.^2 + Ky.^2); % ||K||
W = sqrt(K .* g); % deep water frequencies (empirical parameter)
P = zeros(size(X)); % Cant compute P without loops
for i = 1:meshSize
for j = 1:meshSize
P(i,j) = phillips(Kx(i,j), Ky(i,j), windDir, windSpeed, A, g); % phillips spectrum
end
end
H0 = 1/sqrt(2) .* (randn(size(X)) + 1i .* randn(size(X))) .* sqrt(P); % height field at time t = 0
rotate3d on;
for t = 1:10000 % 10000 * timeStep (sec)
Ht = H0 .* exp(1i .* W .* (t * timeStep)) + ...
conj(flip(flip(H0,1),2)) .* exp(-1i .* W .* (t * timeStep));
[az,el] = view;
surf(X,Y,real(signCor(ifft2(Ht),meshSize)));
axis([-10 10 -10 10 -1e-4 1e-4]); view(az,el);
blue = linspace(0.4, 1.0, 25)'; map = [blue*0, blue*0, blue];
%shading interp; % improve shading (remove "faceted" effect)
colormap(map);
pause(1/60);
end
phillips.m：（我试图向量化菲利普斯谱的计算，但我遇到了一个困难，我将进一步展示）
function P = phillips(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
if k_sq == 0
P = 0;
else
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
if wk < 0
P = 0;
end
end
end
signCor.m：（这个函数对我来说绝对是个谜...我从 "CUDA FFT Ocean Simulation" 实现中复制了它。没有它，模拟效果会更差。而且我也不知道如何向量化这个函数。）
function H = signCor(H1, meshSize)
H = H1;
for i = 1:meshSize
for j = 1:meshSize
if mod(i+j,2) == 0
sign = -1; % works fine if we change signs vice versa
else
sign = 1;
end
H(i,j) = H1(i,j) * sign;
end
end
end
我犯的最大错误是我使用了 ifft 而不是 ifft2，这就是为什么 CUDA ifft 和 Matlab ifft 不匹配的原因。
我的第二个错误是在这行代码中：
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + x(i)); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + y(j)); % = 2*pi*m / Ly
我应该写：
kx = 2.0 * pi / patchSize * (-meshSize / 2.0 + i); % = 2*pi*n / Lx
ky = 2.0 * pi / patchSize * (-meshSize / 2.0 + j); % = 2*pi*m / Ly
我对参数 A、meshSize、patchSize 进行了一些研究，得出的结论是：
波幅的合理参数是 A * (patchSize / meshSize)，其中 A 只是一个比例因子。
为'calm'
patchSize / meshSize <= 0.5.为'tsunami'
patchSize / meshSize >= 3.0.
菲利普斯谱矢量化困难：
我有两个功能：
% non-vectorized spectrum
function P = phillips1(kx, ky, windDir, windSpeed, A, g)
k_sq = kx^2 + ky^2;
if k_sq == 0
P = 0;
else
L = windSpeed^2 / g;
k = [kx, ky] / sqrt(k_sq);
wk = k(1) * windDir(1) + k(2) * windDir(2);
P = A / k_sq^2 * exp(-1.0 / (k_sq * L^2)) * wk^2;
if wk < 0
P = 0;
end
end
end
% vectorized spectrum
function P = phillips2(Kx, Ky, windDir, windSpeed, A, g)
K_sq = Kx .^ 2 + Ky .^ 2;
L = -g^2 / windSpeed^4;
WK = (Kx ./ K_sq) .* windDir(1) + (Ky ./ K_sq) .* windDir(2);
P = (A ./ (K_sq .^ 2)) .* ( exp(L ./ K_sq) .* (WK .^ 2) );
P(K_sq == 0) = 0;
P(WK < 0) = 0;
P(isinf(P)) = 0;
end
在我使用 phillips1 计算 P1 和使用 phillips2 计算 P2 之后，我绘制了它们的差异：
subplot(2,1,1); surf(X,Y,real(P2-P1)); title('Difference in real part');
subplot(2,1,2); surf(X,Y,imag(P2-P1)); title('Difference in imaginary part');
完美地说明了这2个光谱在实部之间存在着巨大的差异。