波浪模拟过程中大规模计算的优化与加速

Question

有一个任务：模拟海面3D情况下的波浪。我遇到的问题是，当模拟大面积海洋（1 km x 1 km）中的波浪时，程序需要很长时间才能执行。事实上它根本没有被执行。因为在等待结果的那半个小时里，从来没有计算过。

我的代码示例：

clear,clc
%% Data 
x = 0:1:100; % coordinates х [m]
y = 0:1:100; % coordinates y [m]
g = 9.81; % gravitational constant [m/s^2]
speed = 5; % wind velocity [m/s]
w0 = (g / speed); % norm.frequency [Hz]
dw = 0.1; % frequency step [rad/s]
w = 0.8:dw:11.1; % frequency [rad/s] 
dtt = pi / 18; % angular step [rad]
theta = 0:dtt:pi; % direction angles, angles between the wavevector & coordintae axis [rad]

%% P-M spectrum, Frequency-Angular spectrum & Amplitude
Psi = 8.1e-3 .* ((w/w0).^(-5)) .* exp((-0.74) ./ ((w/w0).^(4))); % P-M spectrum [none]
Phi = ((speed)^(5)/g^(3)) * Psi; % self-similar spectrum [s*m^2]
Sw = Phi / 2; % frequency spectrum [s*m^2]

St = cos(theta).^(4); % angular spectrum [none]
Norm = trapz(dtt, St); % norm.coefficient [none]
Swt = Sw .* St'; % frequency-angular spectrum [s*m^2]

eta0 = sqrt((Swt * dw * dtt) ./ Norm); % amplitude [m]

figure(1);
subplot(2,1,1)
plot(w, Psi);
title('$$\Psi$$($$\omega$$) - P-M spectrum', 'Interpreter', 'LaTex');
xlabel('\omega [rad/s]');
ylabel('\Psi [none]');
grid on;
subplot(2,1,2)
plot(w, Swt); 
title('$$S(\omega , \theta)$$($$\omega$$) - frequency-angular spectrum', 'Interpreter', 'LaTex');
xlabel('\omega [rad/s]');
ylabel('S(\omega,\theta) [s*m^2]');
grid on;

%% Setting the initial phase parameter
phase = 2*pi*rand(length(theta),length(w)); %% random initial phase ranging from 0 to 2pi [rad]

%% Surface Waving [Linear, 3D (eta & x,y)] at different harmonics & random phase (at one moment in time), different directions of the wavevector (multiple angles)
t = 0; % time moment [s]
Kabs = (w.^2) / g; % wavevector module [rad/m]
Kx = Kabs .* cos(theta)'; % projection of the wavevector onto the x-axis [rad/m]
Ky = Kabs .* sin(theta)'; % projection of the wavevector onto the y-axis [rad/m]

eta = zeros(length(x),length(y),length(theta),length(w)); % reserving space for calculation results
tic
for i = 1:length(x)
    for j = 1:length(y)
        eta(i,j,:,:) = eta0 .* cos(w * t - Kx .* i - Ky .* j + phase);
    end
end
toc
% sum(sum(eta,4),3) - double sum of eta over all harmonics (frequencies) and wavevector directions (angles theta),
% where '4' и '3' summation indicator for variable frequency and angle
etaW = sum(eta,4);
etaWA = sum(etaW,3);

figure(2)
surf(x,y,etaWA);
title('\eta(x,y) - surface waving');
xlabel('x [m]');
ylabel('y [m]');
zlabel('\eta [m]');
cbar = colorbar;
cbar.Label.String = '\eta [m]';
grid on
shading flat

我能够使用的代码优化方法之一是创建一个“空”4D 数组（4D 零数组）

eta = zeros(length(x),length(y),length(theta),length(w));

，在执行循环后，计算结果将被填充到其中：

eta = zeros(length(x),length(y),length(theta),length(w)); % reserving space for calculation results
tic
for i = 1:length(x)
    for j = 1:length(y)
        eta(i,j,:,:) = eta0 .* cos(w * t - Kx .* i - Ky .* j + phase);
    end
end
toc

然后我通过频率和角度变量总结结果：

etaW = sum(eta,4);
etaWA = sum(etaW,3);

从而提前为结果做好准备。这有点帮助。例如，使用此方法的区域

x = 0:1:100; y = 0:1:100;

(100 m x 100 m) 的代码执行时间为 0.7 秒（没有它，则为 3.9 秒）。在区域

x = 0:1:500; y = 0:1:500;

（500 m x 500 m）的情况下，执行时间约为19秒（没有它...我不知道，我没有等待代码执行，但事实证明很长）。然而，对于面积

x = 0:1:1000; y = 0:1:1000;

（1000 m x 1000 m），我很长一段时间都没有得到想要的结果（感觉就像我根本得不到）。

在我的情况下，有没有更多的方法来达到预期的结果并优化我的代码，使其能够应对如此规模的数据（同时，不改变数组中的步骤）？

注意：我的电脑有 16 GB RAM。我最初执行计算的第二台计算机在程序执行过程中完全挂起，因此我必须“手动”重新启动它。所以我想它的内存甚至更少。

Answer 1

您想要一个 2D 输出数组，不要计算完整的 4D 数组然后投影它。相反，计算循环中的投影。

在 1000x1000 输出的情况下，中间 4D 数组超过 15 GB，这意味着您的 16 GB 系统无法将其全部保存在内存中（因为 MATLAB 和操作系统也需要一些空间），并且会将部分交换到磁盘。

如果循环将数组索引为

eta(:,:,i,j)

那么它不会那么糟糕，因为你正在写入内存的连续部分。但由于您索引了

eta(i,j,:,:)

，因此每次循环迭代都会写入 15 GB 内存块中的所有字节。因此，每次循环迭代都会让系统从磁盘读取大部分阵列并将其写回磁盘。这显然会减慢速度，被称为“颠簸”，因为在旋转硬盘的时代，这会产生很大的噪音。

一般来说，在编写循环时，应始终考虑数组元素在内存中的顺序，并尽可能按照内存顺序访问数组元素。这优化了系统缓存的使用。在 MATLAB 中，数组的元素沿第一维连续存储在内存中。因此，您不希望第一个维度索引成为外循环，并且您不希望将最后一个维度上的元素作为单个块进行访问。

在循环内部进行投影，而不是

eta = zeros(length(x),length(y),length(theta),length(w));
for i = 1:length(x)
    for j = 1:length(y)
        eta(i,j,:,:) = eta0 .* cos(w * t - Kx .* i - Ky .* j + phase);
    end
end
etaW = sum(eta,4);
etaWA = sum(etaW,3);

写

etaWA = zeros(length(x),length(y));
for j = 1:length(y)
    for i = 1:length(x)
        etaWA(i,j) = sum(eta0 .* cos(w * t - Kx .* i - Ky .* j + phase), 'all');
    end
end

请注意，我不仅在循环内的所有维度上而不是在末尾添加了

sum

。我还交换了循环顺序，以便内部循环越过第一个维度。这是 MATLAB 中最有效的循环方式。

顺便说一句，“为计算结果预留空间”至关重要。这在 MATLAB 中称为预分配，避免了重复扩大数组时所做的大量工作。