5
$\begingroup$

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): enter image description here

$\endgroup$
  • 5
    $\begingroup$ Can you provide a smaller yet essence-preserving example? $\endgroup$ – Αλέξανδρος Ζεγγ Nov 11 '18 at 6:56
  • $\begingroup$ hey @H.Alanzi welcome! Are you sure "(All of them are a 12 x 9 matrix)" because most of zeros and ones are called with arguments (9, 12) which-if I'm not mistaken-means that they are the transpose of what the quoted claim above states (please check the documentation)... $\endgroup$ – user42582 Nov 11 '18 at 8:22
  • 2
    $\begingroup$ If you just want to make a naive translation, I don't think there'll be any barrier because the MATLAB code only involves basic algebraic calculation, you must have made simple mistake if you don't get what you want. If you want to implement the algorithm in a better way (make the code more elegant, efficient, etc. ) then you should show us the algorithm in traditional math notation, because it's hard to figure it out directly from the code. $\endgroup$ – xzczd Nov 11 '18 at 8:43
17
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ “instructor must have aimed at showing you how not to program in Matlab” – I like your optimism, but going by what atrocities I've seen in quite a few other Matlab scripts, I wouldn't be so sure... $\endgroup$ – leftaroundabout Nov 11 '18 at 23:59
  • $\begingroup$ Is there some packages that we can convert matlab script to mma script, it would be great helpful for many MMA users. $\endgroup$ – ABCDEMMM Nov 17 '18 at 23:01
  • $\begingroup$ @ABCDEMMM Unfortunately, I don't think that such a package exists. It would be really hard to develop. But it is possible to run Matlab code from Mathematica with MatLink. $\endgroup$ – Henrik Schumacher Nov 17 '18 at 23:10
  • $\begingroup$ @HenrikSchumacher, thanks a lot, I have some code about topo, still I want to use it in MMA. okay, thanks a lot again! $\endgroup$ – ABCDEMMM Nov 18 '18 at 0:09
  • $\begingroup$ @ABCDEMMM You're welcome. Glad to be of help. $\endgroup$ – Henrik Schumacher Nov 18 '18 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.