How can use parallel programming in the below code?

Hello. The below code calculates a PDE equation. I'd like to use parallel to speed up in finer discritization.

How is it possible to use parallel in this structure ?

Thanks.

clc
close all

%Time
delt=0.2;
totalTime = 1e3;
timestepn=round(totalTime/delt);%size(delt,2);
%  parameters
L=0.29;
gamma=3;
kappa=2;
% geometry settings; % phase field numbers
mboxsize=64; % system grid numbers in x direction
nboxsize=64; % system grid numbers in y direction
delx=2;      % Delta x
Nx=mboxsize;
Ny=Nx;
dx=delx;
dy=delx;
% calculation-loading data
load etas
ngrain=25;
eta2 = zeros(Nx,Ny,ngrain);
eta = etas;
eta2 = eta;
p=ngrain;
%evolution
for tn=1:10
   
    for i=1:mboxsize
        for j=1:nboxsize
            % calculation of nabla square eta(Laplacian)
            del2=1/delx^2*(0.5*(eta(indg(i+1,nboxsize),j,:)-2*eta(i,j,:)+eta(indg(i-1,nboxsize),j,:))...
                +0.25*(eta(indg(i+2,nboxsize),j,:)-2*eta(i,j,:)+eta(indg(i-2,nboxsize),j,:)))...
                +1/delx^2*(0.5*(eta(i,indg(j+1,mboxsize),:)-2*eta(i,j,:)+eta(i,indg(j-1,mboxsize),:))...
                +0.25*(eta(i,indg(j+2,mboxsize),:)-2*eta(i,j,:)+eta(i,indg(j-2,mboxsize),:)));
    
           
            % summation
            sumterm=eta(i,j,:)*sum(eta(i,j,:).^2)-eta(i,j,:).^3;
            detadtM=(-eta(i,j,:)+eta(i,j,:).^3-kappa*del2);
              
            detadt=-L*(detadtM+2*gamma*(sumterm));
            eta2(i,j,:)=eta(i,j,:)+delt*detadt;
            % for making sure eta is not outside the equilibrium values
            % actually it is unnecessary
            for pind=1:p
                if eta2(i,j,pind)>1
                    eta2(i,j,pind)=1;
                end
                if eta2(i,j,pind)<0
                    eta2(i,j,pind)=0;
                end
            end

        end
    end
   
    eta=eta2;
    phi=sum(eta(:,:,1:p).^2,3);
    imagesc(phi)
    colorbar
    title(tn)

    waitbar(tn/10)
end

 

Answer 1

[Thiago Henrique Gomes Lobato]

You can indeed try to do the code in parallel with a parfor loop in the "i" loop and some tweak in the eta2 variable bound like this:

clc
close all
%Time
delt=0.2;
totalTime = 1e3;
timestepn=round(totalTime/delt);%size(delt,2);
%  parameters
L=0.29;
gamma=3;
kappa=2;
% geometry settings; % phase field numbers
mboxsize=64; % system grid numbers in x direction
nboxsize=64; % system grid numbers in y direction
delx=2;      % Delta x
Nx=mboxsize;
Ny=Nx;
dx=delx;
dy=delx;
% calculation-loading data
load etas
ngrain=25;
eta2 = zeros(Nx,Ny,ngrain);
eta = etas;
eta2 = eta;
p=ngrain;
%evolution
for tn=1:10
   
    parfor i=1:mboxsize
        for j=1:nboxsize
            % calculation of nabla square eta(Laplacian)
            del2=1/delx^2*(0.5*(eta(indg(i+1,nboxsize),j,:)-2*eta(i,j,:)+eta(indg(i-1,nboxsize),j,:))...
                +0.25*(eta(indg(i+2,nboxsize),j,:)-2*eta(i,j,:)+eta(indg(i-2,nboxsize),j,:)))...
                +1/delx^2*(0.5*(eta(i,indg(j+1,mboxsize),:)-2*eta(i,j,:)+eta(i,indg(j-1,mboxsize),:))...
                +0.25*(eta(i,indg(j+2,mboxsize),:)-2*eta(i,j,:)+eta(i,indg(j-2,mboxsize),:)));
   
            % summation
            sumterm=eta(i,j,:)*sum(eta(i,j,:).^2)-eta(i,j,:).^3;
            detadtM=(-eta(i,j,:)+eta(i,j,:).^3-kappa*del2);
              
            detadt=-L*(detadtM+2*gamma*(sumterm));
            eta2(i,j,:)=eta(i,j,:)+delt*detadt;

        end
    end
    % for making sure eta is not outside the equilibrium values
    % actually it is unnecessary
    eta2(eta2>1) = 1;
    eta2(eta2<0) = 0;
   
    eta=eta2;
    phi=sum(eta(:,:,1:p).^2,3);
    imagesc(phi)
    colorbar
    title(tn)
    waitbar(tn/10)
end

I however don't think this is the best solution, since the only loop you actually need to calculate is the time on. If you substitute your i and j loops for vector operations (calculate all matrix values at once) it will probably be way faster than the parallelization and your problem will have only one loop.

[Amir Torabi]

I'd like to solve the probelm for specific values, but, it has error in parfor usage. Could you please take a look at this? The %code is similar the previouse codes.

Thanks.

for tn = 1:10 
[yii,xjj]=find(...
         imerode((phi>0.999),se)==0);
    % space discretization loop
    parfor ii = 1:mboxsize
         i=yii(ii);j=xjj(ii);
         ....
             %code
         
             
    end
    %code
end

[Thiago Henrique Gomes Lobato]

In your new example you're looping for "i" and then replacing it inside the loop

[Amir Torabi]

I modified it.