Parfor loop for finite difference scheme

Hey, i'm trying to convert the for loop(s) in a finite difference code into parfor, but I can't get the right result.

Here's the code:

clc; clear all; close all;

Dy = 1.0204;
Nphi = 150;
Ny = 50;
lambda = 1.8;
Dphi = 0.0422;
phi = linspace(0,2*pi,Nphi);

p(Nphi,Ny)  = 0;
iter = 0;

while (1)
    p_old = p;
    parfor j = 2:Ny-1
        for i = 2:Nphi-1
            B = 1/2;
            D = 1/2;
            E = (eps*sin(phi(i)))/(2*(1+eps*cos(phi(i)))^3)*Dy^2;
            p(i,j) = (1-lambda)*p_old(i,j) + lambda*(B*p_old(i,j-1) + D*p_old(i,j+1) + E)
        end
    end
    ea = max(max(abs(p-p_old)))/max(max(p));
    iter = iter + 1;
    if ea <= 1e-10, break, end
end

I'm not sure if it's appropriate at all to convert it into parfor since the data size isn't large.

[Brendan Hamm]

p_old never changes, so you have the same calculation in every iteration of the while loop. What is the actual mathematical form of the function you are applying the Finite Difference to?

[Andreas]

Okay, where should it be ? I've tried different things but I can't make it work. I'm solving Reynolds equation for a narrow bearing

[Brendan Hamm]

Either I missed the p_old at the beginning of the while loop, or it has been added back in. Obviously this is not the problem. I had thought the update had to do with a difference in time. I will try and look back at this later if no one else gets to it prior. Apologies that I will not have time for the next several days.

Clarifying what you mean by "does not produce the right result will likely be useful information.

At first glance, this may have to do with using eps which is the smallest distance the next distinguishable value is from 1:

That is:

1 + eps == 1

returns false. So this is very small indeed.

As you are considering parallelizing code, have you got this to work in serial prior? Why would you parallelize prior to vectorizing this code? The calculation of E can be done outside of the for loops. Same with B and D.

Looking at the problem it seems it is a series of linear problems. That is if you fix i = 1 then you can represent this problem in a linear form (looking at the first equation):

A*p = b

where A is a tridiagonal matrix, and b is your equation of phi. This can be solved in MATLAB with:

p = A\b;

There is no need to try and implement SOR yourself.

[Tim Berk]

@Brendan

I believe p_old does change every iteration of the while loop.

@Andreas

" I cant get the right result ", Can you elaborate on the result you get vs what you expect? Do you get the right result without the parfor? Why do you want to use parfor? Are you planning to increase the size, or is it a mandatory part of a homework assignment?

In this example it is indeed not appropriate at all. It takes my computer longer to start the parfor loop than to run the cope without parfor.

Answer 1

[Andreas]

Hey, sorry for not answering. I decided to use the fully implicit method, which @Brendan was referring to. This method is WAY faster than the previous one, so it's not really a problem anymore. However, I've read that the grid numbering can be done in way that makes parallel computing beneficial (I think it's called Red-black ordering). Since my grid size is fairly small (200x200) I don't think it's necessary though.

Thanks for the help anyway :)