Matlab/CUDA：海浪模拟
Posted
技术标签:
【中文标题】Matlab/CUDA：海浪模拟【英文标题】：Matlab/CUDA: ocean wave simulation 【发布时间】：2015-04-01 11:50:27 【问题描述】：我研究了 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 海洋模拟源代码可以帮助我理解我的错误？快速谷歌搜索没有得到任何结果。
这是我从“CUDA FFT 海洋模拟”中得到的“正确”结果。我还没有在 Matlab 中实现这种行为，但我已经使用来自“CUDA FFT Ocean Simulation”的数据在 matlab 中绘制了“冲浪”。这是它的样子：
我做了一个实验，得到了一个有趣的结果：
我从“CUDA FFT 海洋模拟”中生成了h0。所以我必须做 ifft 从频域转换到空间域来绘制图形。我已经使用 matlab ifft 并使用来自 CUDA 库的 cufftExecC2C 为相同的 h0 完成了它。结果如下：
CUDA IFFT：
Matlab IFFT：
要么我不明白cufftExecC2C或cufftExecC2C的实现的某些方面，matlab ifft是不同的算法，结果不同。
顺便说一下，生成这种表面的参数是：
meshSize = 32
A = 1e-7
patchSize = 80
windSpeed = 10
【问题讨论】：
你能解释一下是什么让你得到“正确”的结果吗？ 我认为它看起来像一个真正的海洋表面。类似于this 但当然没有纹理。 也许您需要重新设计您的体重秤？在我看来，波浪有点“伸展”到更大的平面上。 符合runDOSrun的评论，可能是采样不足。尝试增加meshSize
你的意思是Y方向？ （如果 XZ 是水平面）
【参考方案1】：
这绝对是一个有趣的练习。这是一个完全重写的答案，因为您自己找到了您要问的问题。
与其删除我的答案，不如发帖帮助您矢量化和/或解释一些代码。
我完全重写了我在以前的答案中给出的 GUI，以便合并您的更改并添加几个选项。它开始长出胳膊和腿，所以我不会把清单放在这里，但你可以在那里找到完整的文件：
ocean_simulator.m.
这是完全独立的，它包括我矢量化的所有计算函数，并在下面单独列出。
图形用户界面将允许您使用参数、为波浪设置动画、导出 GIF 文件（以及一些其他选项，如“预设”，但它们还没有完全解决）。您可以实现的几个示例：
基本
这是您通过快速默认设置和几个渲染选项获得的结果。这使用了较小的网格尺寸和较快的时间步长，因此它在任何机器上都运行得非常快。
我在家里很有限（Pentium E2200 32bit），所以我只能在有限的设置下练习。即使设置为最大值，该 gui 也会运行，但真正享受它会变得很慢。
但是，在工作中快速运行ocean_simulator（I7 64 位，8 核，16GB 内存，Raid 中的 2xSSD），它让它变得更加有趣！以下是几个例子：
虽然在一台更好的机器上完成，但我没有使用任何并行功能或任何 GPU 计算，所以 Matlab 只使用了这些规范的一部分，这意味着它可能在任何 64 位上运行得一样好具有良好 RAM 的系统
风湖
这是一个相当平坦的水面，就像一个湖。即使是大风也不会产生高振幅的波（但仍然有很多小小波）。如果您是一名风帆冲浪者，从山顶上的窗户往外看，您的心脏会漏跳一拍，您的下一步行动就是喊戴夫“伙计！加油。 五分钟后在水上见！"
膨胀
这是你在与暴风雨搏斗了一夜之后，早上从船桥上看的照片。暴风雨已经消散，长长的巨浪是绝对摇摇欲坠的夜晚的最后见证（有航海经验的人会知道......）。
T-风暴
这就是你前一天晚上所做的…… 第二张 gif 在家里完成，因此缺少细节……抱歉
到底部：
最后，gui 将允许您在水域周围添加一个补丁。在 gui 中它是透明的，因此您可以在水下或漂亮的海底添加对象。不幸的是，GIF 格式不能包含 alpha 通道，所以这里没有透明度（但如果你导出到视频中，那么你应该没问题）。
此外，导出到 GIF 会降低图像质量，如果在 Matlab 中运行，域边界和水面之间的连接是完美的。在某些情况下，它还会使 Matlab 降低光照渲染，因此这绝对不是导出的最佳选择，但它允许在 matlab 中播放更多内容。
现在进入代码：
我不会列出完整的 GUI，这会很长（这篇文章已经够长了），我将在这里列出 您的 代码的重写版本，并解释更改。
由于剩余的矢量化，您应该注意到执行速度大幅提高（数量级），但主要有两个原因：
(i) 重复了很多计算。缓存值并重用它们比在循环中重新计算完整矩阵（在动画部分）要快得多。
(ii) 注意我是如何定义表面图形对象的。它只定义一次（甚至为空），然后所有进一步的调用（在循环中）仅更新表面对象的底层ZData（而不是在每次迭代时重新创建表面对象。
这里是：
%% // 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])
您会注意到我更改了一些内容，但我可以向您保证，计算结果完全相同。这段代码调用了一些子函数，但它们都是矢量化的，所以如果你愿意，你可以在这里复制/粘贴它们并内联运行所有内容。
第一个调用的函数是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 参数现在用一个标量值表示，该值表示风的“方位角”（以度数）（从 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 在循环开头for time = (1:nSteps)*timeStep 初始化时，很容易为每行获得正确的time 值。
另请注意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。
我更喜欢对称解决方案（而且它也稍微快一点，但它只计算一次，所以没什么大不了的），所以我使用了我的矢量化方式，但这纯粹是一个选择问题，你绝对可以保持你的方式，它只会稍微“偏移”整个结果矩阵，但它不会干扰计算本身。
第一个答案的最后剩余部分
侧面编程注意事项：
我发现您来自 C/C++ 世界或家庭。在 Matlab 中，您不需要用逗号定义十进制数（例如 2.0，您在大多数数字中都使用了它）。除非另有明确定义，否则 Matlab 默认将任意数字转换为 double，这是一个 64 位浮点类型。所以写2 * pi 就足以获得最大精度（Matlab 不会将 pi 转换为整数 ;-)），你不需要写2.0 * pi。 虽然如果你不想改变你的习惯，它仍然有效。
另外，（Matlab 的一大优点），在运算符前添加. 通常意味着“元素级”操作。您可以通过这种方式明智地添加 (.+)、减 (.-)、乘 (.*)、除 (./) 完整矩阵元素。这就是我摆脱代码中所有循环的方法。这也适用于幂运算符：A.^2 将返回一个与A 大小相同且每个元素均平方的矩阵。
【讨论】：
非常感谢您的回答！我真的很感激你所做的巨大工作。但也许你可以再解释一个问题：在执行你的代码时，我注意到可以创建的波非常小，但在文章中它们的大小相当大。如果你对这个问题感兴趣并且知道 CUDA 编程语言，那么在“CUDA FFT Ocean Simulation”中有这个算法的实现。 @EugeneKolesnikov，是的，我注意到 matlab 代码的平均波长明显小于论文示例中的波长。我最初想在它上面添加几个更长的波长分量只是为了模拟某种更大的“膨胀”，但我在你的编辑中看到实际上这是由 iFFT 结果的差异引起的。我不太精通 CUDA 编码（而且我目前的 AMD 显卡也不支持它），但我确信有办法让它们匹配。 @EugeneKolesnikov，您能否在线分享一个在 iFFT 之前使用的矩阵数据H0 的示例，以及 cufftExecC2C 变换的结果？ （我可以尝试让 Matlab 得到同样的结果）。
嘿@AnderBiguri。是的，我对那个有点忘乎所以……但这对我来说真的很有趣。你试过代码吗？ （模拟您的花园池塘）
仅供参考：.+ 和 .- 是 Octave 扩展，在 MATLAB 中不起作用。除此之外......哇，答案需要付出很多努力。也非常令人印象深刻的结果。太棒了！【参考方案2】：
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 只是一个比例因子。
为了“冷静”patchSize / meshSize <= 0.5。
对于“海啸”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 和使用P2 使用phillips2 计算它们的区别之后：
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');
它完美地说明了这两个频谱在实际部分之间存在巨大差异。
【讨论】：
很好地整理了你的计算。您仍然可以获得很多速度改进。查看我修改后的答案。 （注意：你的P 不是 complex 所以检查虚部是没有意义的）以上是关于Matlab/CUDA：海浪模拟的主要内容，如果未能解决你的问题，请参考以下文章