改进5D矩阵计算

Question

大家早上好，

我写了这段代码：

clc
clear all
close all

%Box size
Nx=4096;
Ny=15;
Nz=15;

%Spatial gird resolution
delta=6;

%WT / turbulence condition
UHub=11.4;
HubHt=90;
z0=0.03;
IECturbC='B';

%%INITIALISATION

% definition of constants
twopi=2*pi;
fourpi=4*pi;
sqrt2=sqrt(2);

%constants and derived parameters from IEC
gamma = 3.9; %IEC, (B.12)
alpha = 0.2; %IEC, sect. 6.3.1.2

%set delta1 according to guidelines (chap.6)
if HubHt<=60,
    delta1=0.7*HubHt;
else
    delta1=42;
end;

%IEC, Table 1, p.22
if IECturbC == 'A',
    Iref=0.16;
elseif IECturbC == 'B',
    Iref=0.14; 
elseif IECturbC == 'C',
    Iref=0.12;
else
    error('IECturbC can be equal to A,B or C;adjust the input value')
end;

%IEC, sect. 6.3.1.3
b=6.5;
sigma1=Iref*(0.75*UHub+b);
%derived constants
l=0.8*delta1; %IEC, (B.12)
sigmaiso=0.55*sigma1; %IEC, (B.12)

%%MAIN PROGRAM
Cij=zeros(3,3,Nx,Ny,Nz);
k = zeros(3,1); %current vector k

for ikx=1:(Nx),
    m = -1.*Nx/2+ikx;
    k(1)=m*l/(Nx*delta)*twopi;
    for iky=1:(Ny),
        m= -1.*Ny/2+iky;
        k(2)=m*l/(Ny*delta)*twopi;
        for ikz=1:(Nz),
        m= -1.*Nz/2+ikz;
        k(3)=m*l/(Nz*delta)*twopi;


           if k(1)==0,
            Cij(:,:,ikx,iky,ikz)=0;
             else
            kabs=sqrt(k(1)^2+k(2)^2+k(3)^2);
            beta= gamma./(kabs.^(2/3));
            k0(3)=k(3)+beta.*k(1);
            k0abs=sqrt(k(1)^2+k(2)^2+k0(3)^2);
            Ek0=1.453*k0abs^4/(1.+k0abs.^2)^(17/6);
            C1=beta.*k(1)^2*( k0abs.^2 - 2*k0(3)^2 + beta.*k(1)*k0(3) )/( kabs.^2*( k(1)^2 + k(2)^2 ));
            C2=k(2).*k0abs.^2./ (exp( (3/2).*log( k(1).^2 + k(2).^2 ) )) .* atan2( beta.*k(1).* sqrt( k(1)^2 + k(2)^2 ) ,( k0abs.^2 - k0(3).*k(1).*beta));
            xhsi1=C1 - k(2).*C2./k(1);
            xhsi2=k(2).*C1./k(1) + C2;

            Cij(1,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( k(2).*xhsi1);
            Cij(1,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( k(3) - k(1).*xhsi1 + beta.*k(1));
            Cij(1,3,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( -k(2));
            Cij(2,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( k(2).*xhsi2 - k(3) - beta.*k(1));
            Cij(2,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( -k(1).*xhsi2);
            Cij(2,3,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( k(1));
            Cij(3,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( k0abs.^2.*k(2) ./ (kabs.^2));
            Cij(3,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*( -k0abs.^2*k(1) ./ (kabs.^2));
            Cij(3,3,ikx,iky,ikz)= 0;
           end;       
        end;
    end;
end;

我想问你: 1.有没有更快的方法来获得Cij矩阵？当Nx，Ny，Nz增加时，Cij的计算相当慢；2.有没有办法得到plot(kabs，beta)和plot(kabs，Ek0)？

请对我耐心点，在matlab的世界里，我还是个新手。

提前感谢并致以最良好的问候，弗朗西斯科

Answer 1

如果你想在matlab中有好的性能，你应该尝试尽可能多地对你的代码进行向量化。

例如，不是执行以下操作：

for x=1:n
    A(x)=x^2
end

做

x=1:n;
A=x.^2;

当您有多个索引时，可以使用ndgrid。因此，与其这样做，不如：

for x=1:nx
  for y=1:ny
    for z=1:nz
      A(x,y,z)=x^2+y-2*z;
    end
  end
end

做

[x y z]=ndgrid(1:nx,1:ny,1:nz)
A=x.^2+y-2*z

由于您看起来很努力，我已经为您更改了代码。执行时间现在是0.33秒。矢量化的版本是：

clc
clear all
close all
tic

%Box size
Nx=1024;
Ny=15;
Nz=15;

%Spatial gird resolution
delta=6;

%WT / turbulence condition
UHub=11.4;
HubHt=90;
z0=0.03;
IECturbC='B';

%%INITIALISATION

% definition of constants
twopi=2*pi;
fourpi=4*pi;
sqrt2=sqrt(2);

%constants and derived parameters from IEC
gamma = 3.9; %IEC, (B.12)
alpha = 0.2; %IEC, sect. 6.3.1.2

%set delta1 according to guidelines (chap.6)
if HubHt<=60,
    delta1=0.7*HubHt;
else
    delta1=42;
end;

%IEC, Table 1, p.22
if IECturbC == 'A',
    Iref=0.16;
elseif IECturbC == 'B',
    Iref=0.14; 
elseif IECturbC == 'C',
    Iref=0.12;
else
    error('IECturbC can be equal to A,B or C;adjust the input value')
end;

%IEC, sect. 6.3.1.3
b=6.5;
sigma1=Iref*(0.75*UHub+b);
%derived constants
l=0.8*delta1; %IEC, (B.12)
sigmaiso=0.55*sigma1; %IEC, (B.12)

Cij2=zeros(3,3,Nx,Ny,Nz);
[x y z]=ndgrid(1:Nx,1:Ny,1:Nz);
k1=(x-Nx/2)*l/(Nx*delta)*twopi;
k2=(y-Ny/2)*l/(Ny*delta)*twopi;
k3=(z-Nz/2)*l/(Nz*delta)*twopi;
kabs=sqrt(k1.^2+k2.^2+k3.^2);
beta= gamma./(kabs.^(2/3));
k03=k3+beta.*k1;
k0abs=sqrt(k1.^2+k2.^2+k03.^2);
Ek0=1.453*k0abs.^4./(1+k0abs.^2).^(17/6);
C1=beta.*k1.^2.*( k0abs.^2 - 2*k03.^2 + beta.*k1.*k03 )./( kabs.^2.*( k1.^2 + k2.^2 ));
C2=k2.*k0abs.^2./ (exp( (3/2).*log( k1.^2 + k2.^2 ) )) .* atan2( beta.*k1.* sqrt( k1.^2 + k2.^2 ) ,( k0abs.^2 - k03.*k1.*beta));
xhsi1=C1 - k2.*C2./k1;
xhsi2=k2.*C1./k1 + C2;
CC=sigmaiso*sqrt(twopi*pi*l^3.*Ek0./(Nx*Ny*Nz*delta^3.*k0abs.^4));
Cij2(1,1,:,:,:)= CC.*( k2.*xhsi1);
Cij2(1,2,:,:,:)= CC.*( k3 - k1.*xhsi1 + beta.*k1);
Cij2(1,3,:,:,:)= CC.*( -k2);
Cij2(2,1,:,:,:)= CC.*( k2.*xhsi2 - k3 - beta.*k1);
Cij2(2,2,:,:,:)= CC.*( -k1.*xhsi2);
Cij2(2,3,:,:,:)= CC.*( k1);
Cij2(3,1,:,:,:)= CC.*( k0abs.^2.*k2 ./ (kabs.^2));
Cij2(3,2,:,:,:)= CC.*( -k0abs.^2.*k1 ./ (kabs.^2));

Answer 2

我将尝试更全面地回答您的问题，以便其他人也能从中受益。

您应该对for循环进行矢量化，以便更快地编写代码。而不是像这样的东西：

for i=1:n
    for j=1:m
        M(i,j)=sqrt(i) + sqrt(j);
    end
end

根据以下代码向量化循环：

[xi,xj] = ndgrid(1:n,1:m);
M = sqrt(xi)+sqrt(xj);