ELECTRIC POTENTIAL - How to speed up the code? (COMPLETE CODE INCLUDED)

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

Answer 1

[Walter Roberson]

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.

[matico]

Thank you for your response, Walter.

I've tried both methods which you suggested and both of them gave very similar results to my initial method. I stopped at "difference == 0.01" and the end result is almost the same compared to my method.

I think I should use the most recently calculated values when calculating new V(i,j) values. Thus, I think your first suggestion should be appropriate. But nevertheless, both of your suggestions give the same result, so probably it's not relevant for this case.

Both of your methods are much faster, too. Both ended in a few minutes, while my method takes few hours.

Thanks again, I truly appreciate your time and effort.