ELECTRIC POTENTIAL - How to speed up the code? (COMPLETE CODE INCLUDED)
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
Hi.
I have a code which calculates an electrical potential around stormy cloud, which is placed above the lighthouse. I compute this using "Finite difference method". A "box" is an 320x320 matrix with known boundary conditions (V=0) and V of a cloud is also known (V1). V of a lighthouse is 0.
Now I have 2 nested "for" loops which go from left to right side of a matrix and from top to bottom. Every round I compare new and old value of matrix. When they differentiate for less than 0.01 I am happy and quit.
That method works OK, except it is very slow. It needs more than 4 hours to complete. Do you have any idea how to speed it up? Is possible not to use nested "for" loops? I read a lot about vectorization, but still don't know how to put this into my code.
I include a complete code, so you are welcome to look at it.
Thanks
clc; clear;
% Cloud Potential
V1 = 6500;
% Cloud Radius
CloudRad = 20;
% Cloud Height
CloudHeight = 40;
% Number of points (vertically * horizontally)
n = 320;
% Height and Radius of a Space (width is 2 * rho)
h = 80;
rho = 40;
% Space Width - Cloud Width
Width = 2*rho - 2*CloudRad;
% Lighthouse measurements
rad1 = 6;
rad2 = 7;
h1 = 22;
h2 = 17;
r = linspace(0,2*rho,n);
z = linspace(h,0,n);
% Mesh
[rr,zz] = meshgrid(r,z);
% Cloud measurements based on mesh
RelCloudHeight = round(((h-CloudHeight)/h) * n);
if RelCloudHeight == 0
RelCloudHeight = 1;
end
RelCloudWidth = round(((2*rho-Width)/(2*rho))*n);
% Matrix of calculated Potential
V = zeros(n);
% Cloud is always on V1
V(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
% Boundary conditions - Potential there is 0
V(1:n,1) = zeros(n,1);
V(1:n,n) = zeros(n,1);
% Lighthouse measurements based on mesh
UpperHight = (round(((h-h1)/h)*n));
MiddleHight = (round((h-h2)/h*n));
ExtremeLeft = round(((rho-rad1)/(2*rho))*n)+1;
ExtremeRight = round(((rho+rad1)/(2*rho))*n);
Left = round(((rho-rad2)/(2*rho))*n);
Right = round(((rho+rad2)/(2*rho))*n);
Finite difference method
difference = 100;
while difference > 0.01 % Desired final accuracy
VV = V;
for i = 2:n-1 % Along the rows of matrix V
for j = 2:n-1 % Along the columns of matrix V
% Lighthouse is always on V=0
V(UpperHight:MiddleHight,ExtremeLeft:ExtremeRight) = 0;
V(MiddleHight:n,Left:Right) = 0;
% Cloud is always on V1
V(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
% Iteration equation
V(i,j) = (1/6)*(2*V(i,j+1)+2*V(i,j-1)+V(i+1,j)+V(i-1,j));
end
end
% Watching the difference between old and new value of V
difference = max(max(abs(V-VV)));
end
% At the end I set the predetermined elements of V to the required value
% for the last time
V(UpperHight:MiddleHight,ExtremeLeft:ExtremeRight) = 0;
V(MiddleHight:n,Left:Right) = 0;
V(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
save LighthouseOUT.mat
PLOTTING
figure;
EkvNum = 150; % Number of equipotential contours
contour(rr,zz,V,EkvNum); axis equal;
% Plotting area
xlim([0 2*rho]); ylim([0 h]);
set(gca,'fontsize',14,'fontname','Times','box','on');
title(['Equipotentials (' num2str(EkvNum) ')']);
xlabel('{\itWidth} / m');
ylabel('{\itHeight} / m');
% Plotting lighthouse:
line([rho-rad2 rho-rad2],[0 h2],'LineWidth',2,'Color','black');
line([rho+rad2 rho+rad2],[0 h2],'LineWidth',2,'Color','black');
line([rho-rad2 rho-rad1],[h2 h2],'LineWidth',2,'Color','black');
line([rho+rad2 rho+rad1],[h2 h2],'LineWidth',2,'Color','black');
line([rho-rad1 rho-rad1],[h2 h1],'LineWidth',2,'Color','black');
line([rho+rad1 rho+rad1],[h2 h1],'LineWidth',2,'Color','black');
line([rho-rad1 rho+rad1],[h1 h1],'LineWidth',2,'Color','black');
% Plotting cloud:
line([rho-CloudRad rho+CloudRad],[CloudHeight CloudHeight],...
'LineWidth',2,'Color','black');
0 Commenti
Risposta accettata
Walter Roberson
il 9 Ago 2015
I am puzzled as to why you are doing
% Lighthouse is always on V=0
V(UpperHight:MiddleHight,ExtremeLeft:ExtremeRight) = 0;
V(MiddleHight:n,Left:Right) = 0;
% Cloud is always on V1
V(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
for every (i,j) pair . You only need to do that if the previous (i,j) was in that region, in which case you would be undoing the effect of the assignment to V(i,j) for that previous (i,j). The only time it is worth doing and then immediately undoing an assignment is the boundary case, the very last (i,j) as that would not loop back up to have that region assigned.
I have to ask you though whether that is intentional, that you are counting on the possibility that V(n,n) is within the regions that get assigned every iteration and want your "difference" variable to reflect that V(n,n) has not been re-assigned the 0 or V1 value? It looks to me as if it would not be expected that V(n,n) is in Lighthouse or Cloud region; if I am correct, then the effect of those loops would be the same as setting the values in those regions once before the nested "for" and then skip the assignments if (i,j) is within those regions.
If you combine this with some logical indexing...
%do once before we start iterations
% Lighthouse is always on V=0
V(UpperHight:MiddleHight,ExtremeLeft:ExtremeRight) = 0;
V(MiddleHight:n,Left:Right) = 0;
% Cloud is always on V1
V(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
%M is true in the locations V is to be updated
M = true(n,n);
% Lighthouse is always on V = 0
M(UpperHight:MiddleHight,ExtremeLeft:ExtremeRight) = false;
M(MiddleHight:n,Left:Right) = false;
% Cloud is always on V1
M(RelCloudHeight,(n/2-RelCloudWidth/2)+1:(n/2+RelCloudWidth/2)) = V1;
%and we never update the first or last row or column
M([1,n],:) = false;
M(:,[1,n]) = false;
%now we have to find the indices in column order because the original code changed column most quickly
[Vrow, Vcol] = find(M);
Vrowcol = sortrows([Vrow, Vcol]);
Vidx = sub2ind(size(V), Vrowcol(:,1), Vrowcol(:,2));
%we use linear indexing. V(i,j+1) is n elements further on than V(i,j), V(i+1,j) is 1 element further on than V(i,j)
while difference > 0
VV = V;
for idx = Vidx
V(idx) = (1/6)*(2*V(idx+n)+2*V(idx-n)+V(idx+1)+V(idx-1));
end
difference = max(abs(V(:)-VV(:)));
end
I have replicated in this code that the value at V(i,j) depends upon the changed value of V(i,j-1) and V(i-1,j) and upon the unchanged (yet) V(i,j+1) and V(i+1,j) . The sortrows() and sub2idx() is there to order the iterations the same way that you had before, proceeding first to last row, first to last column, using the changed values from above-or-left and the unchanged values from below-or-right.
If that was an accident of coding and each iteration should be using the unchanged values, reading from the VV rather from the V, then the calculation becomes completely vectorizable and with simpler setup logic:
%we can find the indices in any order as V is "read-only" for calculating the next iteration
Vidx = find(M);
%we use linear indexing. V(i,j+1) is n elements further on than V(i,j), V(i+1,j) is 1 element further on than V(i,j)
while difference > 0
VV = V;
V(Vidx) = (1/6)*(2*V(Vidx+n)+2*V(Vidx-n)+V(Vidx+1)+V(Vidx-1)); %Everything at once
difference = max(abs(V(:)-VV(:)))
end
I think you would find this second version much faster if it matches what the update equations should be. To re-emphasize, this version is for the case where the reference to V(i,j-1) is intended to be to the V(i,j-1) from before the nested loop, rather than to the V(i,j-1) that was just assigned on the previous "for j" iteration.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Graphics Objects in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!