An efficient way to translate a triple for-loop from Matlab to Mathematica

Question

This is my first post on StackExchange, so if I violate any etiquettes, I apologize in advance. I am a senior engineering student who has been using Mathematica for a year now. For one of my classes, called reservoir simulation, I am trying to obtain the pressure distribution matrix. My instructor exclusively uses Matlab to model the problems we do in class. I am not very good at Matlab, so I decided to use Mathematica.

The instructor gave us his code in Matlab, and I want to translate it into Mathematica. There are a lot of constants (All of them are a 12 x 9 matrix) in the problem I managed to get all of them into my Mathematica notebook, so they are all defined. What I will not share are the constant but assume they are already in the notebook. Below is the Matlab code my instructor used which I am explicitly trying to translate.

First, he initialized many matrices like this:

P_old = ones(9,12) * 7750;

Pwf = 5600 ; % psia

Pe = zeros(9,12);

E = zeros(9,12);

Ge = zeros(9,12);

SGe = zeros(9,12);

Pw = zeros(9,12);

W = zeros(9,12);

Gw = zeros(9,12);

SGw = zeros(9,12);

Pn = zeros(9,12);

N = zeros(9,12);

Gn = zeros(9,12);

SGn = zeros(9,12);

Ps = zeros(9,12);

S = zeros(9,12);

Gs = zeros(9,12);

SGs = zeros(9,12);

omega = zeros(9,12); % productivity index

omega_history = zeros(9,12,61);

Store = zeros(9,12,61);

q = zeros(9,12);

q(3,9) = -650;

dt = 1;

A = zeros(108,108);

B = zeros(108,1);

gamma = zeros(9,12);

nx = 12;

ny = 9;

Then he made a very long triple for-loop. The indices for this for-loop are problem specific which is what I am trying to replicate.

for t = 1:61

Store(:,:,t) = P_old;

for j = 1:9

     for i = 1:11


        Pe(j,i) = (Vb(j,i+1) * P_old(j,i+1) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j,i+1));

        E(j,i) = (2 * kx1(j,i+1) * Ax(j,i+1) * kx1(j,i) * Ax(j,i) ) / ((kx1(j,i+1) * Ax(j,i+1) * x(j,i)) + (kx1(j,i) * Ax(j,i) * x(j,i+1))) * (1 + 9.0 * 10^(-6) * (Pe(j,i) - 14.7)) * 1/(  a * (Pe(j,i))^(3) + b * (Pe(j,i))^(2) + c * Pe(j,i) + d  );

        Ge(j,i) = G(j,i+1) - G(j,i);

        SGe(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pe(j,i) - 14.7))) );


     end 
end 

Pe(isnan(Pe)) = 0;

E(isnan(E)) = 0;

SGe(isnan(SGe)) = 0;
for j = 1:9

  for i = 2:12

        Pw(j,i) = (Vb(j,i-1) * P_old(j,i-1) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j,i-1));

        W(j,i) = (2 * kx1(j,i-1) * Ax(j,i-1) * kx1(j,i) * Ax(j,i) ) / ((kx1(j,i-1) * Ax(j,i-1) * x(j,i)) + (kx1(j,i) * Ax(j,i) * x(j,i-1))) *  (1 + 9.0 * 10^(-6) * (Pw(j,i) - 14.7)) * 1/(  a * (Pw(j,i))^(3) + b * (Pw(j,i))^(2) + c * Pw(j,i) + d  ); 

        Gw(j,i) = G(j,i-1) - G(j,i);

        SGw(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pw(j,i) - 14.7))) );
  end 
end 
Pw(isnan(Pw)) = 0; 

W(isnan(W)) = 0;

SGw(isnan(SGe)) = 0;
for j = 2:9

  for i = 1:12

       Pn(j,i) = (Vb(j-1,i) * P_old(j-1,i) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j-1,i));

       N(j,i) = (2 * ky1(j-1,i) * Ay(j-1,i) * ky1(j,i) * Ay(j,i) ) / ((ky1(j-1,i) * Ay(j-1,i) * y(j,i)) + (ky1(j,i) * Ay(j,i) * y(j-1,i))) * (1 + 9.0 * 10^(-6) * (Pn(j,i) - 14.7)) * 1/(  a * (Pn(j,i))^(3) + b * (Pn(j,i))^(2) + c * Pn(j,i) + d  );

       Gn(j,i) = G(j-1,i) - G(j,i);

      SGn(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Pn(j,i) - 14.7))) );

  end
end

Pn(isnan(Pn)) = 0;
N(isnan(N)) = 0;
SGn(isnan(SGn)) = 0;

for j = 1:8

           for i = 1:12

       Ps(j,i) = (Vb(j+1,i) * P_old(j+1,i) + Vb(j,i) * P_old(j,i)) / ( Vb(j,i) + Vb(j+1,i));

       S(j,i) = (2 * ky1(j+1,i) * Ay(j+1,i) * ky1(j,i) * Ay(j,i) ) / ((ky1(j+1,i) * Ay(j+1,i) * y(j,i)) + (ky1(j,i) * Ay(j,i) * y(j+1,i))) * (1 + 9.0 * 10^(-6) * (Ps(j,i) - 14.7)) * 1/(  a * (Ps(j,i))^(3) + b * (Ps(j,i))^(2) + c * Ps(j,i) + d  );                           

       Gs(j,i) = G(j+1,i) - G(j,i);

       SGs(j,i) = (1/144) * (refdensity * ((1 + 9.0 * 10^(-6) .* (Ps(j,i) - 14.7))) );

           end 
     end 


Ps(isnan(Ps)) = 0;   
S(isnan(S)) = 0;
SGs(isnan(SGs)) = 0;


if P_old(4,4) <= Pwf 

omega(4,4) = 0; 
else 

omega(4,4) = ( (2 * pi .* sqrt(kx1(4,4) * ky1(4,4)) .* thickness(4,4) * (1 + 9.0 * 10^(-6) * (P_old(4,4) - 14.7))) / ( (  a * (P_old(4,4))^(3) + b * (P_old(4,4))^(2) + c * P_old(4,4) + d  ) * ( log(re/0.25) ) ) ) ;

end 


omega_history(:, :, t) = omega;


for j = 1:9

   for i = 1:12

      C(j,i) = - (E(j,i) + W(j,i) + N(j,i) + S(j,i) + gamma(j,i) + omega(j,i));

       C(C==0) = 1;

    end 
      end 

    for j = 1:9

      for i = 1:12

        Q(j,i) = (-1 .* omega(j,i) .* Pwf) - (gamma(j,i) .* P_old(j,i)) + ( E(j,i) .* SGe(j,i) .* Ge(j,i) )  + ( W(j,i) .* SGw(j,i) .* Gw(j,i) )  + ( N(j,i) .* SGn(j,i) .* Gn(j,i) )  + ( S(j,i) .* SGs(j,i) .* Gs(j,i) )  - q(j,i);

        end 
         end

        Q(isnan(Q)) = 0;

   for j = 1:9   

     for i = 1:12


     r = i + nx * (j-1);



         A(r,r) = C(j,i);


        if r >1 


            A(r,r-1) = W(j,i);

        end 

        if r < nx * ny 

            A(r,r+1) = E(j,i);
        end 

        if r + nx <= nx * ny 

            A(r,r+nx) = S(j,i);

        end 

        if r > nx

            A(r,r-nx) = N(j,i);

        end 

        B(r) = Q(j,i);

 end
 end 


 P_new = (A\B);



 P_new = reshape(P_new,[12,9]);
 P_new = P_new';



 P_old = P_new;                



 end

As you can see above, it is a very involved messy code. I attempted to do this in Mathematica using the "Module" function and the "Do" function, but I cannot replicate the results. My futile attempt is the following:

Module[{i, j}, Do[pe = (volume[[j, i + 1]] pi[[j, i + 1]] + 
  volume[[j, i]]  pi[[j, i]])/(volume[[j, i]] + 
   volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet, {i, 1, 8, 1}, {j, 1, 10, 1}]]

This is just one small for-loop of the bigger for-loop, and I do not get what I want.

The end (P_new) result should be something like this (The "*" represents an indeterminate):

Answer 1

Okay, let's focus on this snippet:

Module[{i, j}, 
  Do[
   pe = (volume[[j, i + 1]] pi[[j, i + 1]] + 
    volume[[j, i]]  pi[[j, i]])/(volume[[j, i]] + 
    volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet, 
  {i, 1, 8, 1}, {j, 1, 10, 1}]]

Creating some dummy data:

volume = RandomReal[{-1, 1}, {10, 9}];
pi = RandomReal[{-1, 1}, {10, 9}];

You use Do here and not For. That's indeed great for many reasons and shows that you already know about the problems with For. Moreover, Do scopes its iterators, hence we don't need Module.

So let's strip it off. Moreover, by default, Do increases its iterators by 1

Do[
 pe = (volume[[j, i + 1]] pi[[j, i + 1]] + 
       volume[[j, i]] pi[[j, i]])/(volume[[j, i]] + 
       volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
 {i, 1, 8}, {j, 1, 10}]

This reassigns pe over and over again. But I got the impression that you mean pe to be a matrix. The "Matlab-way" of constructing it would be to allocate a 0-matrix and to fill it afterwards (actually, a good Matlab programmer would do it differently):

pe = ConstantArray[0, {10, 8}];
Do[pe[[j, i]] = (volume[[j, i + 1]] pi[[j, i + 1]] + 
       volume[[j, i]] pi[[j, i]])/(volume[[j, i]] + 
       volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
 {i, 1, 8, 1}, {j, 1, 10, 1}]

In Mathematica, we can use Table instead.

pe = Table[(volume[[j, i + 1]] pi[[j, i + 1]] + 
       volume[[j, i]] pi[[j, i]])/(volume[[j, i]] + 
       volume[[j, i + 1]]) //. Indeterminate -> "*" // Quiet,
  {j, 1, 10, 1}, {i, 1, 8, 1}]

Notice that I had to reorder the iterators in order to produce the same matrix as before and not its transpose.

Both Matlab and Mathematica are interpreted languages (in contrast to compiled languages like C or fortran). And in interpreted languages, it is usually way more performant to avoid loop constructions if possible and to replace them by vectorized code. Your instructor must have aimed at showing you how not to program in Matlab. Indeed, the code is overly convoluted and almost all interior loops could be replaced by vectorized operations (at the first glance).

So, in this example, the following produces the same result, and is one to two orders of magnitude faster:

pe = With[{a = volume pi},
  (a[[;; , ;; -2]] + a[[;; , 2 ;;]])/(volume[[;; , ;; -2]] + volume[[;; , 2 ;;]])
  ]/. Indeterminate -> "*" // Quiet

This should also lead to much more concise and more legible code.

Another example (it is so striking that I believe that the instructor did that intentionally):

for j = 1:9

   for i = 1:12

      C(j,i) = - (E(j,i) + W(j,i) + N(j,i) + S(j,i) + gamma(j,i) + omega(j,i));

       C(C==0) = 1;

    end 
      end 

would be simply (in Matlab syntax):

C = - (E + W + N + S + gamma + omega);
C(C==0) = 1;

Moreover, putting C(C==0) = 1; into the body of the loop is just plain crazy; it has to be performed just once.

Remarks

For the understanding of ;; see Span; it is essentially the analogue to Matlab's :. (Notice the three-argument versions of these functions differ in the ordering of the arguments.) See also this post for more analogies between Matlab and Mathematica.