matlab - Matlab/CUDA : ocean wave simulation

Question

我研究了杰里·特森多夫(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海洋模拟”中获得的“正确”结果。我尚未在Matlab中实现此行为，但已使用来自“CUDA FFT海洋模拟”的数据在matlab中绘制了“冲浪”。看起来是这样的:



我做了一个实验，得到了一个有趣的结果:

我已经从“C​​UDA FFT海洋模拟”中提取了生成的 h0。因此，我必须进行ifft从频域转换到空间域才能绘制图形。我已经使用matlab h0和CUDA库中的 ifft对相同的 cufftExecC2C完成了此操作。结果如下:

CUDA ifft:

Matlab ifft:

我要么不了解 cufftExecC2C的实现的某些方面，要么不了解 cufftExecC2C和matlab ifft是具有不同结果的不同算法。

顺便说一下，生成这种表面的参数是:

meshSize = 32

A = 1e-7

patchSize = 80

windSpeed = 10

Answer 1

好吧，那绝对是一个有趣的练习。这是一个完全重写的答案，因为您发现了自己要问的问题。

除了删除我的答案外，发布还可以帮助您向量化和/或解释一些代码。

我完全重写了我以前回答中给出的GUI，以便合并您的更改并添加几个选项。它开始变得越来越长，所以我不会在这里列出 list ，但是您可以在此处找到完整的文件:

ocean_simulator.m 。

这是完全独立的，包括我矢量化的所有计算功能，并在下面分别列出。

GUI将允许您播放参数，设置波形动画，导出GIF文件(以及其他一些选项(例如“预设”)，但还没有解决)。您可以实现的一些示例:

基本的

这就是使用快速默认设置和几个渲染选项所获得的。这使用较小的网格大小和快速的时间步长，因此它可以在任何计算机上快速运行。



我在家里很受限制(奔腾E2200 32位)，所以我只能在有限的设置下练习。 gui将在最大设置下运行，但真正享受起来会变慢。

但是，通过快速运行 ocean_simulator (I7 64位，8核，16GB内存，Raid中的2xSSD)，它变得更加有趣!这里有一些例子:

尽管是在更好的机器上完成的，但我没有使用任何并行功能，也没有进行GPU计算，因此Matlab仅使用了部分规格，这意味着它可以在具有适当RAM的任何64位系统上运行得同样好

风大湖

这是一个相当平坦的水面，就像一个湖泊。即使是强风也不会产生高振幅的波(但是仍然会产生许多小波)。如果您是一名风帆冲浪者，则从山顶上的窗户看风，您的心脏将略微跳动，下一步是叫戴夫“ Man!装备起来。在水面上与您见面五次!”


胀

整夜都在与 Storm 抗争之后，这是您早晨从船桥上眺望的景色。 Storm 消散了，长长的大浪是绝对黑夜的最后见证者(有帆船经验的人会知道...)。


Storm

这就是你前一天晚上所做的...

第二张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之间的任何值)。后来，由于函数 X，它被翻译为其 Y和 pol2cart组件。

[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

最后三行计算需要一些解释，因为它们收到的变化最大(所有结果相同，但速度要好得多)。

您的原始行:

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。

我倾向于对称解决方案(而且它也稍快一些，但只计算一次，所以没什么大不了的)，所以我使用了矢量化方法，但这纯粹是一个选择问题，您绝对可以保持自己的方式，它只会稍微“偏移”整个结果矩阵，但不会干扰计算本身。

第一个答案的最后剩余
侧面编程说明:
我检测到您来自C/C++世界或家庭。在Matlab中，您不需要用逗号定义十进制数字(例如 2.0，您将其用于大多数数字)。除非另有明确定义，否则Matlab默认将任何数字强制转换为 double，这是64位浮点类型。因此，编写 2 * pi足以获得最高的精度(Matlab不会将pi转换为整数;-))，而无需编写 2.0 * pi。尽管如果您不想改变自己的习惯，仍然可以使用。

同样，(Matlab的一大优势)在操作符之前添加 .通常意味着“按元素”操作。您可以通过这种方式明智地添加( .+)，减( .-)，乘以( .*)，除( ./)完整矩阵元素。这就是我摆脱了代码中所有循环的方式。这对于幂运算符也适用: A.^2将返回与 A相同大小的矩阵，每个元素均平方。