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 = (AB);







 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):

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 it iterators, so 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 intenionally):

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.