Улучшение расчета 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. Есть ли способ получить сюжет (kabs, beta) и сюжет (kabs, Ek0)?

Пожалуйста, будьте терпеливы со мной, я все еще новичок в мире Matlab.

Спасибо заранее и наилучшими пожеланиями, Francesco

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 секунды. Vecotrized версия:

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);