Matlab/CUDA: моделирование океанских волн

Вопрос:

Я изучил “Моделирование океанской воды” в статье Джерри Тессендорфа и попытался запрограммировать статистическую волновую модель, но я не получил правильный результат, и я не понимаю, почему.

В моей программе я попытался создать поле высоты волны в момент 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” . Я еще не достиг такого поведения в Matlab, но я остановил “серфинг” в Matlab, используя данные “CUDA FFT Ocean Simulation” . Вот как это выглядит:

Я сделал эксперимент и получил интересный результат:

Я взял сгенерированный h0 из “CUDA FFT Ocean Simulation” . Поэтому я должен сделать ifft для преобразования из частотной области в пространственную область для построения графика. Я сделал это для того же h0, используя matlab ifft и используя cufftExecC2C из библиотеки CUDA. Здесь результат:

CUDA ifft:

Matlab ifft:

Либо я не понимаю некоторые аспекты реализации cufftExecC2C или cufftExecC2C, а matlab ifft – это разные алгоритмы с разными результатами.

Кстати, параметры для создания такой поверхности:

meshSize = 32
A = 1e-7
patchSize = 80
windSpeed = 10

Лучший ответ:

Ну, это было забавное упражнение. Это полностью переписанный ответ, так как вы нашли проблемы, о которых вы спрашивали сами.

Вместо того, чтобы удалять мой ответ, по-прежнему есть заслуга в публикации, чтобы помочь вам векторизовать и/или объяснить несколько бит кода.

Я полностью переписал GUI, который я дал в своем прежнем ответе, чтобы включить ваши изменения и добавить несколько вариантов. Он начал расти руками и ногами, поэтому я не стану перечислять здесь, но вы можете найти полный файл там:

ocean_simulator.m.

Это полностью самодостаточно и включает в себя все вычислительные функции, которые я вектурировал и перечислял отдельно ниже.

Графический интерфейс позволит вам играть с параметрами, анимировать волны, экспортировать GIF файл (и несколько других опций, таких как “пресет”, но они еще не сглажены). Несколько примеров того, что вы можете достичь:

Содержание

Основные
Ветренное озеро
набухают
Т-шторм
Внизу:
Теперь на код:

Основные

Это то, что вы получаете с быстрыми настройками по умолчанию и несколькими вариантами рендеринга. Это использует небольшой размер сетки и быстрый временной шаг, поэтому он работает довольно быстро на любой машине.

Я довольно ограничен дома (Pentium E2200 32bit), поэтому я мог заниматься только с ограниченными настройками. Gui будет работать даже с настройками maxed, но он станет медленным, чтобы действительно наслаждаться.

Однако при быстром запуске ocean_simulator на работе (I7 64 бит, 8 ядер, 16 ГБ RAM, 2xSSD в рейде), это делает его намного больше весело! Вот несколько примеров:

_{Несмотря на то, что сделано на гораздо более качественной машине, я не использовал никаких параллельных функций и вычислений графических процессоров, поэтому Matlab использовал только часть этих спецификаций, а это значит, что он мог бы работать так же хорошо на любых 64-битных система с достойным ОЗУ}

Ветренное озеро

Это довольно плоская водная поверхность, как озеро. Даже сильные ветры не создают волны с высокой амплитудой (но все же много мини-всплесков). Если вы являетесь владельцем ветра, смотрящим на это из окна на вершине холма, ваше сердце собирается пропустить удар, и ваш следующий шаг – позвонить Дейву “ Man!”. Встречайте вас в пять на воде! ”

набухают

Это ты смотришь с моста своей лодки утром, после того, как всю ночь сражался со штормом. Буря рассеялась, и длинные волны стали последним свидетелем того, что было определенно шаткой ночью (люди с опытом парусного спорта узнают…).

Т-шторм

И это было то, что вы были до ночи раньше…
_{второй gif, выполненный дома, отсюда отсутствие деталей… извините}

Внизу:

Наконец, gui позволит вам добавить патч вокруг водного домена. В gui это прозрачно, чтобы вы могли добавлять объекты под водой или красивое дно океана. К сожалению, формат GIF не может включать альфа-канал, поэтому здесь нет прозрачности (но если вы экспортируете в видео, тогда вы должны быть в порядке).

Кроме того, экспорт в 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. Эта функция заканчивает расчеты волнового поля за заданное время. Это определенно стоит того, что это само по себе, потому что это набор вычислений mimimun, который нужно будет повторять за каждый данный момент времени, когда вы хотите оживить волны.

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 больше соответствует вашей строке.

Я предпочитаю симметричные решения (плюс это также немного быстрее, но он вычисляется только один раз, так что это не очень важно), поэтому я использовал свой векторный способ, но это чисто вопрос выбора, вы можете определенно держите свой путь, он только когда-либо слегка “смещает” всю матрицу результатов, но он не возмущает вычисления как таковые.

_{последние остатки первого ответа}
Боковые программные заметки:
Я обнаружил, что вы пришли из мира или семьи C/С++. В Matlab вам не нужно определять десятичное число с комой (например, 2.0, вы использовали это для большинства своих номеров). Если не указано иное, Matlab по умолчанию передает любой номер в double, который является 64-битным типом с плавающей запятой. Поэтому писать 2 * pi достаточно, чтобы получить максимальную точность (Matlab не будет отличать pi как целое;-)), вам не нужно писать 2.0 * pi. Хотя он все равно будет работать, если вы не захотите изменить свои привычки.

Кроме того, (одно из больших преимуществ Matlab), добавляя ., прежде чем оператор обычно означает операцию “по элементам”. Вы можете добавить (.+), substract (.-), умножить (.*), разделить (./) полный матричный элемент таким образом. Вот как я избавился от всех циклов в вашем коде. Это также работает для оператора питания: A.^2 вернет матрицу того же размера, что и A, с каждым квадратом элемента.

Ответ №1

Здесь рабочая программа.

Прежде всего – исходный код:

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: (Я попытался векторизовать вычисление спектра Phillips, но я столкнулся с трудностью, которую я покажу дальше)

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.

Сложность с векторизацией спектра Филлипса:

У меня есть 2 функции:

% 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

После вычисления P1 с помощью phillips1 и 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');

Это прекрасно иллюстрирует, что существует огромная разница между этими двумя спектрами в реальной части.