Я пытаюсь усреднить свои данные CFD (которые представлены в виде скалярного массива N x M x P; N соответствует Y, M к x и P к z) за подмножество временных шагов.Я попытался упростить описание моего желаемого процесса усреднения ниже.
- Поворот сетки на каждом временном шаге на заданный угол ( это потому, что поток имеет когерентную структуру, которая вращаетсяи изменяет форму / размер на каждом временном шаге, и я хочу перекрыть их и найти усредненную по времени форму структуры, которая учитывает изменение формы / размера во времени )
- Рисование сферы по центруна исходной не повернутой сетке
- Определение точек сетки из всех повернутых сеток, которые лежат в сфере
- Определение индексов точек сетки в каждой повернутой сетке
- Использованиеиндексы для поиска скалярных данных в повернутых точках сетки внутри сферы
- Возьмите среднее значение в сфере
- Поместите это новое усредненное значение в местоположение на неповоротной сетке
У меня есть код, который, кажется, правильно делает то, что я хочу, но это занимает слишком много времени, чтобы закончить вычисления.Я хотел бы, чтобы он работал быстрее, и я открыт для изменения кода в случае необходимости.Ниже приведена версия моего кода, которая работает с уменьшенной версией данных.
x = -5:5; % x - position data
y = -2:.5:5; % y - position data
z = -5:5; % z - position data
% my grid is much bigger actually
[X,Y,Z] = meshgrid(x,y,z); % mesh for plotting data
dX = diff(x)'; dX(end+1) = dX(end); % x grid intervals
dY = diff(y)'; dY(end+1) = dY(end); % y grid intervals
dZ = diff(z)'; dZ(end+1) = dZ(end); % z grid intervals
TestPoints = combvec(x,y,z)'; % you need the Matlab Neural Network Toolbox to run this
dXYZ = combvec(dX',dY',dZ')';
% TestPoints is the unrotated grid
M = length(x); % size of grid x - direction
N = length(y); % size of grid y - direction
P = length(z); % size of grid z - direction
D = randi([-10,10],N,M,P,3); % placeholder for data for 3 time steps (I have more than 3, this is a subset)
D2{3,M*N*P} = [];
PosAll{3,2} = [];
[xSph,ySph,zSph] = sphere(50);
c = 0.01; % 1 cm
nu = 8e-6; % 8 cSt
s = 3*c; % span for Aspect Ratio 3
r_g = s/sqrt(3);
U_g = 110*nu/c; % velocity for Reynolds number 110
Omega = U_g/r_g; % angular velocity
T = (2*pi)/Omega; % period
dt = 40*T/1920; % time interval
DeltaRotAngle = ((2*pi)/T)*dt; % angle interval
timesteps = 121:123; % time steps 121, 122, and 123
for ti=timesteps
tj = find(ti==timesteps);
Theta = ti*DeltaRotAngle;
Rotate = [cos(Theta),0,sin(Theta);...
0,1,0;...
-sin(Theta),0,cos(Theta)];
PosAll{tj,1} = (Rotate*TestPoints')';
end
for i=1:M*N*P
aa = TestPoints(i,1);
bb = TestPoints(i,2);
cc = TestPoints(i,3);
rs = 0.8*sqrt(dXYZ(i,1)^2 + dXYZ(i,2)^2 + dXYZ(i,3)^2);
handles.H = figure;
hs = surf(xSph*rs+aa,ySph*rs+bb,zSph*rs+cc);
[Fs,Vs,~] = surf2patch(hs,'triangle');
close(handles.H)
for ti=timesteps
tj = find(timesteps==ti);
f = inpolyhedron(Fs,Vs,PosAll{tj,1},'FlipNormals',false);
TestPointsR_ti = PosAll{tj,1};
PointsInSphere = TestPointsR_ti(f,:);
p1 = [aa,bb,cc];
p2 = [PointsInSphere(:,1),...
PointsInSphere(:,2),...
PointsInSphere(:,3)];
w = 1./sqrt(sum(...
(p2-repmat(p1,size(PointsInSphere,1),1))...
.^2,2));
D_ti = reshape(D(:,:,:,tj),M*N*P,1);
D2{tj,i} = [D_ti(f),w];
end
end
D3{M*N*P,1} = [];
for i=1:M*N*P
D3{i} = vertcat(D2{:,i});
end
D4 = zeros(M*N*P,1);
for i=1:M*N*P
D4(i) = sum(D3{i}(:,1).*D3{i}(:,2))/...
sum(D3{i}(:,2));
end
D_ta = reshape(D4,N,M,P);
Я ожидаю получить массив N x M x P, где каждый индекс является средневзвешенным значением всех точек, охватывающих всевременные шаги в этой конкретной позиции в не повернутой сетке.Как видите, это именно то, что я получаю.Основная проблема, однако, заключается в том, сколько времени требуется для этого, когда я использую больший набор моих «реальных» данных.Код выше занимает всего пару минут, но когда M = 120, N = 24 и P = 120, а количество временных шагов равно 24, это может занять гораздо больше времени.По моим оценкам, для завершения всего расчета потребуется около 25+ дней.