Hi everybody, i have to do a homework which simulates fire spreading and the quest is; convert serial matlab code to parallel and then into CUDA (GPU computing). The serial matlab code is :
_________________________________________________________________
function celluralFire292(fire_spots, time_steps, CAwidth, CAheight, noDynamics)
%2D Cellular Automaton For Fire Simulation
% 'time_steps' are the number of steps to simulate.
%
% 'fire_spots' is the number of initial fire spots.
%
% 'CAwidth' and 'CAheight' are the width and height of automaton,
% respectively.
%
% 'noDynamics' makes everything homogeneous
%----------------------------------------------------------------
% FIRST STEP. CHECK INPUT ARGUMENTS FOR VALIDITY
%----------------------------------------------------------------
% Check number of arguments.
if (nargin > 5)
error('Way too many arguments passed to function!');
end
close all; %close open figures
% Fill in potentially unset optional values with default ones.
if (nargin <1)
fire_spots=3;
end
if (nargin <2)
time_steps=300;
end
if (nargin <3)
CAwidth=120;
end
if (nargin <4)
CAheight=120;
end
if (nargin <5)
noDynamics=false;
end
%Check fire_spots
if isempty(fire_spots)
error('Invalid input parameter fire_spots (empty)!');
end
if (~isscalar(time_steps))
error('Invalid input parameter fire_spots. Must be a scalar value!');
end
if (~isnumeric(fire_spots))
error('Invalid input parameter fire_spots. Must be a numeric value!');
end
if ( fire_spots<1)
error('Input parameter fire_spots out of range. Valid values greater or equal to 1!');
end
%Check time_steps
if isempty(time_steps)
error('Invalid input parameter time_steps (empty)!');
end
if (~isscalar(time_steps))
error('Invalid input parameter time_steps. Must be a scalar value!');
end
if (~isnumeric(time_steps))
error('Invalid input parameter time_steps. Must be a numeric value!');
end
if ( time_steps<1)
error('Input parameter time_steps out of range. Valid values greater or equal to 1!');
end
%Check CAwidth
if isempty(CAwidth)
error('Invalid input parameter CAwidth (empty)!');
end
if (~isscalar(CAwidth))
error('Invalid input parameter CAwidth. Must be a scalar value!');
end
if (~isnumeric(CAwidth))
error('Invalid input parameter CAwidth. Must be a numeric value!');
end
if ( CAwidth<5)
error('Input parameter CAwidth out of range. Valid values greater or equal to 5!');
end
%Check CAheight
if isempty(CAheight)
error('Invalid input parameter CAheight (empty)!');
end
if (~isscalar(CAheight))
error('Invalid input parameter CAheight. Must be a scalar value!');
end
if (~isnumeric(CAheight))
error('Invalid input parameter CAheight. Must be a numeric value!');
end
if ( CAheight<5)
error('Input parameter CAheight out of range. Valid values greater or equal to 5!');
end
%Check noDynamics
if isempty(noDynamics)
error('Invalid input parameter noDynamics (empty)!');
end
if (~isscalar(noDynamics))
error('Invalid input parameter noDynamics. Must be a scalar value!');
end
if (~islogical(noDynamics))
error('Invalid input parameter noDynamics. Must be a logical value!');
end
%----------------------------------------------------------------
%TRUTH AT LAST. START CODING CELLURAL AUTOMATON
%----------------------------------------------------------------
width=CAwidth;
height=CAheight;
C=zeros(width,height); %cell states. Initialized to zero (non-burning)
%initialize fire spots
if (fire_spots==1)
C(floor(width/2),floor(height/2))=1; %a fire spot
else
spots_i=randi([2,width-1],1,fire_spots);
spots_j=randi([2,height-1],1,fire_spots);
for k=1:fire_spots
C(spots_i(k),spots_j(k))=1; %a fire spot
end
end
%height coefficients
H=zeros(width,height);
x=1:width;
y=1:height;
H=10+0.13*(sin(5.5*pi*x/width)')*sin(7.2*pi*y/height);
H=H+1.1*(sin(1.5*pi*x/width)')*sin(1.2*pi*y/height);
H=H+0.2*(sin(1.12*pi*x/width)')*sin(1.3*pi*y/height);
H=H+0.15*(sin(0.2*pi*x/width)')*sin(0.52*pi*y/height);
%Temperatures spatial coefficients.
Fuel=1+0.2*randn(width,height);
Fuel=Fuel- min(0,min(min(Fuel)) );
Fuel=Fuel/max(max(Fuel));
if ~(noDynamics)
jk_greenmap=[0.1, 0.4, 0.1; 0.1,1,0.1]; %maps 0 to black and 1 to green
surf(x,y,H,Fuel);
colormap(jk_greenmap);
title('Topology and Vegetation');
end
%Wind coefficients with respect to each neighbor
W=zeros(3,3);
windStrength=0.4; %0 to 1
windVec=[1,1]; %(x,y) coordinates of end point. Start point assumed (x,y)=(0,0)
%weights of neighbors proportional to projection of
%wind vector to the neighbor's direction
for m=-1:1
for n=-1:1
nVec=[m n]; %vector pointing to neighbor
nVec=nVec/norm(nVec); %normalize length for fairness
W(2+n,2-m)=windStrength*nVec*(windVec')/norm(windVec); %dot product. Depends on wind length and angle
end
end
%distance weights (for circular distribution)
DIS=ones(3,3);
DIS(1,1)=sqrt(2);
DIS(1,3)=sqrt(2);
DIS(3,1)=sqrt(2);
DIS(3,3)=sqrt(2);
% Cellural Automata Output Movie - Intialization to zeros.
movie_out = zeros(width, height, time_steps);
movie_out( :, :, 1) = C; %initial output state
%I assume that the boundary of te used tables are dixed
%boundary conditions. Automaton is defined within this boundary.
disp('PATIENCE MY FRIEND!');
%-------------------------------------------------
% CELLURAL SIMULATION. I HAVE NO CLUE WHETHER THIS IS WHAT
% YOU WERE ASKING FOR, BUT IT IS A CELLURAL AUTOMATON
% ANYWAYS (ALL RULES ARE LOCAL, both for updating Cnew and Rnew)...
%-------------------------------------------------
stop_frame=time_steps;
Cnew=C;
for k = 2:time_steps
stop_simulation=true; %used to stop simulation when all burned or stable...
%simulate everything but boundary
for i=2:width-1
for j=2:height-1
if (C(i,j)==1) %skip if totally burned
continue;
end
%update CA
Cnew(i,j)=0;
%iteratively apply rule using all neighbors
for m=-1:1
for n=-1:1
%update neighbors weight
if (m==0 && n==0)
if (noDynamics)
R=1/9;
else
R=1/9;
end
else
if (noDynamics)
R=4/(9*DIS(2+m,2+n));
else
Hdiff=1+( H(i,j)-H(i+m,j+n) );
R=1.4*(1+Hdiff*Fuel(i,j)*W(2+m,2+n))/(9*DIS(2+m,2+n));
end
end
Cnew(i,j)=Cnew(i,j)+R*C(i+m,j+n);
end
end
%make sure you update
if ( ( Cnew(i,j)- C(i,j) )< 0.1 )
Cnew(i,j)=C(i,j)*1.1;
end
%make sure you do not exceed 1.0
Cnew(i,j)=min(1,Cnew(i,j));
if ~(Cnew(i,j)==1)
stop_simulation=false;
end
end
end
C=Cnew;
movie_out
,:,k)=C; %update current state array for next iteration
if (stop_simulation )
stop_frame=k;
break
end
end
%-------------------------------------------------------------------------
% FINAL STEP. PLOT THE OUTPUT AS MOVIE.
%-------------------------------------------------------------------------
figure;
%plot and captute every frame of the simulation output
for n=1:stop_frame
jk_colormap=[1,1,1;0.2,0.2,0.2]; %maps 0 to white and 1 to black (dark gray)
%pcolor(cellural_out);
%imshow(cellural_out, 'InitialMagnification', 100);
imagesc(movie_out
,:,n)); %image values scaled automatically
%colormap(jk_colormap);
axis ij; %places the coordinate system origin in the upper left corner
% The cell array ( {...} ) in title is used to write separate lines...
title(['\...}']);
text(3,3,['\color[rgb]{0 .3 .6}Time Step: ', num2str(n) ] );
set(gca,'XAxisLocation','top');
axis image; %keeps relative dimensions of axes proportional to image size
%capture and store current plot frame
M(n)=getframe;
end
figure;
imagesc(movie_out
,:,2));
figure;
imagesc(movie_out
,:,50));
figure;
imagesc(movie_out
,:,100));
% %%CREATE A MOVIE FROM FRAMES
% desired_playback_time=7.0; %desired playback time in seconds (approximately)
% fps=(round(1.0*stop_frame/desired_playback_time));
% fps=max(1,fps); %make sure its higher than 1
% fps=min(25,fps); %make sure its less than 50
% fps=min(fps,stop_frame); %make sure its less than stop_frame
% movie(M,1,fps);
%Save output to a movie file. Quite problematic...
%movie2avi(M, [imageName , '.avi'], 'compression', 'None', 'fps', 2 );
end
%PLEASE DONT USE THE CODE FOR COMMERCIAL PURPOSES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
_________________________________________________________________
I have done something with parallel matlab conversion but i am not sure is completely correct. Any suggestions (for both parallel matlab and CUDA conversion)?
->You can also send me email.
_________________________________________________________________
function celluralFire292(fire_spots, time_steps, CAwidth, CAheight, noDynamics)
%2D Cellular Automaton For Fire Simulation
% 'time_steps' are the number of steps to simulate.
%
% 'fire_spots' is the number of initial fire spots.
%
% 'CAwidth' and 'CAheight' are the width and height of automaton,
% respectively.
%
% 'noDynamics' makes everything homogeneous
%----------------------------------------------------------------
% FIRST STEP. CHECK INPUT ARGUMENTS FOR VALIDITY
%----------------------------------------------------------------
% Check number of arguments.
if (nargin > 5)
error('Way too many arguments passed to function!');
end
close all; %close open figures
% Fill in potentially unset optional values with default ones.
if (nargin <1)
fire_spots=3;
end
if (nargin <2)
time_steps=300;
end
if (nargin <3)
CAwidth=120;
end
if (nargin <4)
CAheight=120;
end
if (nargin <5)
noDynamics=false;
end
%Check fire_spots
if isempty(fire_spots)
error('Invalid input parameter fire_spots (empty)!');
end
if (~isscalar(time_steps))
error('Invalid input parameter fire_spots. Must be a scalar value!');
end
if (~isnumeric(fire_spots))
error('Invalid input parameter fire_spots. Must be a numeric value!');
end
if ( fire_spots<1)
error('Input parameter fire_spots out of range. Valid values greater or equal to 1!');
end
%Check time_steps
if isempty(time_steps)
error('Invalid input parameter time_steps (empty)!');
end
if (~isscalar(time_steps))
error('Invalid input parameter time_steps. Must be a scalar value!');
end
if (~isnumeric(time_steps))
error('Invalid input parameter time_steps. Must be a numeric value!');
end
if ( time_steps<1)
error('Input parameter time_steps out of range. Valid values greater or equal to 1!');
end
%Check CAwidth
if isempty(CAwidth)
error('Invalid input parameter CAwidth (empty)!');
end
if (~isscalar(CAwidth))
error('Invalid input parameter CAwidth. Must be a scalar value!');
end
if (~isnumeric(CAwidth))
error('Invalid input parameter CAwidth. Must be a numeric value!');
end
if ( CAwidth<5)
error('Input parameter CAwidth out of range. Valid values greater or equal to 5!');
end
%Check CAheight
if isempty(CAheight)
error('Invalid input parameter CAheight (empty)!');
end
if (~isscalar(CAheight))
error('Invalid input parameter CAheight. Must be a scalar value!');
end
if (~isnumeric(CAheight))
error('Invalid input parameter CAheight. Must be a numeric value!');
end
if ( CAheight<5)
error('Input parameter CAheight out of range. Valid values greater or equal to 5!');
end
%Check noDynamics
if isempty(noDynamics)
error('Invalid input parameter noDynamics (empty)!');
end
if (~isscalar(noDynamics))
error('Invalid input parameter noDynamics. Must be a scalar value!');
end
if (~islogical(noDynamics))
error('Invalid input parameter noDynamics. Must be a logical value!');
end
%----------------------------------------------------------------
%TRUTH AT LAST. START CODING CELLURAL AUTOMATON
%----------------------------------------------------------------
width=CAwidth;
height=CAheight;
C=zeros(width,height); %cell states. Initialized to zero (non-burning)
%initialize fire spots
if (fire_spots==1)
C(floor(width/2),floor(height/2))=1; %a fire spot
else
spots_i=randi([2,width-1],1,fire_spots);
spots_j=randi([2,height-1],1,fire_spots);
for k=1:fire_spots
C(spots_i(k),spots_j(k))=1; %a fire spot
end
end
%height coefficients
H=zeros(width,height);
x=1:width;
y=1:height;
H=10+0.13*(sin(5.5*pi*x/width)')*sin(7.2*pi*y/height);
H=H+1.1*(sin(1.5*pi*x/width)')*sin(1.2*pi*y/height);
H=H+0.2*(sin(1.12*pi*x/width)')*sin(1.3*pi*y/height);
H=H+0.15*(sin(0.2*pi*x/width)')*sin(0.52*pi*y/height);
%Temperatures spatial coefficients.
Fuel=1+0.2*randn(width,height);
Fuel=Fuel- min(0,min(min(Fuel)) );
Fuel=Fuel/max(max(Fuel));
if ~(noDynamics)
jk_greenmap=[0.1, 0.4, 0.1; 0.1,1,0.1]; %maps 0 to black and 1 to green
surf(x,y,H,Fuel);
colormap(jk_greenmap);
title('Topology and Vegetation');
end
%Wind coefficients with respect to each neighbor
W=zeros(3,3);
windStrength=0.4; %0 to 1
windVec=[1,1]; %(x,y) coordinates of end point. Start point assumed (x,y)=(0,0)
%weights of neighbors proportional to projection of
%wind vector to the neighbor's direction
for m=-1:1
for n=-1:1
nVec=[m n]; %vector pointing to neighbor
nVec=nVec/norm(nVec); %normalize length for fairness
W(2+n,2-m)=windStrength*nVec*(windVec')/norm(windVec); %dot product. Depends on wind length and angle
end
end
%distance weights (for circular distribution)
DIS=ones(3,3);
DIS(1,1)=sqrt(2);
DIS(1,3)=sqrt(2);
DIS(3,1)=sqrt(2);
DIS(3,3)=sqrt(2);
% Cellural Automata Output Movie - Intialization to zeros.
movie_out = zeros(width, height, time_steps);
movie_out( :, :, 1) = C; %initial output state
%I assume that the boundary of te used tables are dixed
%boundary conditions. Automaton is defined within this boundary.
disp('PATIENCE MY FRIEND!');
%-------------------------------------------------
% CELLURAL SIMULATION. I HAVE NO CLUE WHETHER THIS IS WHAT
% YOU WERE ASKING FOR, BUT IT IS A CELLURAL AUTOMATON
% ANYWAYS (ALL RULES ARE LOCAL, both for updating Cnew and Rnew)...
%-------------------------------------------------
stop_frame=time_steps;
Cnew=C;
for k = 2:time_steps
stop_simulation=true; %used to stop simulation when all burned or stable...
%simulate everything but boundary
for i=2:width-1
for j=2:height-1
if (C(i,j)==1) %skip if totally burned
continue;
end
%update CA
Cnew(i,j)=0;
%iteratively apply rule using all neighbors
for m=-1:1
for n=-1:1
%update neighbors weight
if (m==0 && n==0)
if (noDynamics)
R=1/9;
else
R=1/9;
end
else
if (noDynamics)
R=4/(9*DIS(2+m,2+n));
else
Hdiff=1+( H(i,j)-H(i+m,j+n) );
R=1.4*(1+Hdiff*Fuel(i,j)*W(2+m,2+n))/(9*DIS(2+m,2+n));
end
end
Cnew(i,j)=Cnew(i,j)+R*C(i+m,j+n);
end
end
%make sure you update
if ( ( Cnew(i,j)- C(i,j) )< 0.1 )
Cnew(i,j)=C(i,j)*1.1;
end
%make sure you do not exceed 1.0
Cnew(i,j)=min(1,Cnew(i,j));
if ~(Cnew(i,j)==1)
stop_simulation=false;
end
end
end
C=Cnew;
movie_out
if (stop_simulation )
stop_frame=k;
break
end
end
%-------------------------------------------------------------------------
% FINAL STEP. PLOT THE OUTPUT AS MOVIE.
%-------------------------------------------------------------------------
figure;
%plot and captute every frame of the simulation output
for n=1:stop_frame
jk_colormap=[1,1,1;0.2,0.2,0.2]; %maps 0 to white and 1 to black (dark gray)
%pcolor(cellural_out);
%imshow(cellural_out, 'InitialMagnification', 100);
imagesc(movie_out
%colormap(jk_colormap);
axis ij; %places the coordinate system origin in the upper left corner
% The cell array ( {...} ) in title is used to write separate lines...
title(['\...}']);
text(3,3,['\color[rgb]{0 .3 .6}Time Step: ', num2str(n) ] );
set(gca,'XAxisLocation','top');
axis image; %keeps relative dimensions of axes proportional to image size
%capture and store current plot frame
M(n)=getframe;
end
figure;
imagesc(movie_out
figure;
imagesc(movie_out
figure;
imagesc(movie_out
% %%CREATE A MOVIE FROM FRAMES
% desired_playback_time=7.0; %desired playback time in seconds (approximately)
% fps=(round(1.0*stop_frame/desired_playback_time));
% fps=max(1,fps); %make sure its higher than 1
% fps=min(25,fps); %make sure its less than 50
% fps=min(fps,stop_frame); %make sure its less than stop_frame
% movie(M,1,fps);
%Save output to a movie file. Quite problematic...
%movie2avi(M, [imageName , '.avi'], 'compression', 'None', 'fps', 2 );
end
%PLEASE DONT USE THE CODE FOR COMMERCIAL PURPOSES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
_________________________________________________________________
I have done something with parallel matlab conversion but i am not sure is completely correct. Any suggestions (for both parallel matlab and CUDA conversion)?
->You can also send me email.