Suggestions to increase the speed of a fugacity calculation

Question

I'm writing code to simulate a chemical reaction. This requires the calculation of the fugacity of a mixture of gases. The underlying model is the SRK equation of state.

First the physical properties are declared:

*fugacity calculations according to SRK equation*)
(*Critical properties for selected gases*)
r = 8.314;
presc = {35.0, 73.8, 20.5, 220.5, 81.0, 33.9}*100000;
tempc = {132.9, 304.2, 43.6, 647.3, 512.6, 126.2};
volc = {93.1, 94.0, 51.5, 56.0, 118.0}*10^-6;
omegac = {0.049, 0.255, 0.0, 0.344, 0.572, 0.040};

The equation of state is set up:

(*SRK equation*)
a = 0.42748*(r*tempc)^2/presc;
b = 0.08664*(r*tempc)/presc;
alpha[temp_] = (1 + (0.480 + 1.574 omegac - 0.176 omegac^2) (1- 
Sqrt[temp/tempc]))^2;

A[temp_, pres_] = Table[(a[[i]]*alpha[temp][[i]]*pres)/(r*temp)^2, {i, 1, 
6}];

B[temp_, pres_] = Table[(b[[i]]*pres)/(r*temp), {i, 1, 6}];

Amix[temp_, pres_, yfrac_] := Sum[yfrac[[i]]*yfrac[[j]]*Sqrt[A[temp, 
pres][[i]]*A[temp, pres][[j]]], {i, 1, 6}, {j, 1, 6}];

Bmix[temp_, pres_, yfrac_] := Sum[yfrac[[i]]*B[temp, pres][[i]], {i, 1, 6}];

Finally, the fugacity is calculated:

(Fugacity calculation)

fugacitymix[temp_, pres_, i_, yfrac_] := With[{Amix2 = Amix[temp, pres, 
yfrac], Bmix2 = Bmix[temp, pres, yfrac]},
                        Exp[(z - 1) B[temp, pres][[i]]/Bmix2 - Log[z - 
Bmix2] - Amix2/Bmix2 (2 (A[temp, pres][[i]]/Amix2)^0.5 - 
B[temp, pres][[i]]/Bmix2) Log[1 + Bmix2/z] /. 
z -> Last[
  Cases[z /. {ToRules@NRoots[z^3 - 
         z^2 + (Amix2 - Bmix2 - (Bmix2)^2) z - (Amix2) (Bmix2) == 
        0, z]}, _Real]]
                        ]]

The issue is that the code is somewhat slow. It makes running the model costly. For example:

yfra = N@{3, 20, 75, 1, 1, 1}/100;

Plot[Table[fugacitymix[273 + 200, x*100000, i, yfra], {i, 1, 6}], {x, 
1, 300}] // AbsoluteTiming

Takes approximately 40s on my machine.

I tried compiling the function with:

cfugacity = 
Compile[{{temp, _Real}, {pres, _Real}, {i, _Integer}, {yfrac, _Real, 
6}},
With[{Amix2 = Amix[temp, pres, yfrac], 
Bmix2 = Bmix[temp, pres, yfrac]},
                        Exp[
                        (z - 1) B[temp, pres][[i]]/Bmix2 - Log[z - Bmix2] - 
  Amix2/Bmix2 (2 (A[temp, pres][[i]]/Amix2)^0.5 - 
     B[temp, pres][[i]]/Bmix2) Log[1 + Bmix2/z] /. 
 z -> Last[
   Cases[z /. {ToRules@
       NRoots[z^3 - 
          z^2 + (Amix2 - Bmix2 - (Bmix2)^2) z - (Amix2) (Bmix2) ==
          0, z]}, _Real]]
                            ]]
  ]

Plot[Table[cfugacity[273 + 200, x*100000, i, yfra], {i, 1, 6}], {x, 
1, 300}] // AbsoluteTiming

This one takes in fact longer than the uncompiled version, about 47s. I'm looking for suggestions to improve the speed.

Answer 1

The functions Amix and Bmix are called quite often with the same arguments. Memoization (see below) is a easily implementable way to speed the computations up. Moreover, the sum in Bmix can be rewritten as a simple Dot-product and the double sum in Amix can be also simplified tremendously. Try your code with these new definitions; it makes the execution about 80 times faster on my machine:

ClearAll[Amix, Bmix, A, B];

A[temp_, pres_] := (a alpha[temp]) (pres/(r*temp)^2);
B[temp_, pres_] := b (pres/(r*temp));

Bmix[temp_, pres_, yfrac_] := Bmix[temp, pres, yfrac] = yfrac.B[temp, pres];
Amix[temp_, pres_, yfrac_] := Amix[temp, pres, yfrac] = (yfrac.Sqrt[A[temp, pres]])^2;

Edit

I reworked the code a bit in the meantime. Due to Daniel Lichtblau's suggestion, I added memoization also for A, B, and alpha. Moreover, I removed the dependence of A and B on pres (they simply depend linearly on pres) so that fewer memoized values have to be computed. But most importantly, I vectorized fugacitymix along the argument i and I observed that the value of z can be reused for all i. These are the new parts of the code:

ClearAll[Amix0, Bmix0, A0, B0, alpha]

alpha[temp_] := alpha[temp] = (1. + (0.480 + 1.574 omegac - 0.176 omegac^2) (1 - Sqrt[temp/tempc]))^2;

A0[temp_] := A0[temp] = (a alpha[N[temp]]) (1./(r N[temp])^2);
B0[temp_] := B0[temp] = b (1./(r N[temp]));

Amix0[temp_, yfrac_] := Amix0[temp, yfrac] = (yfrac.Sqrt[A0[temp]])^2;
Bmix0[temp_, yfrac_] := Bmix0[temp, yfrac] = yfrac.B0[temp]

fugacitymix0[temp_, pres_?NumericQ, yfrac_] := Block[{Amix2, Bmix2, z, zz},
  Amix2 = Amix0[temp, yfrac] pres;
  Bmix2 = Bmix0[temp, yfrac] pres;
  z = Last[
    Cases[
     z /. {ToRules@NRoots[
           z^3 - z^2 + (Amix2 - Bmix2 - (Bmix2)^2) z - (Amix2) (Bmix2) == 0,
           z]}, _Real]
    ];
  Exp[(z - 1) pres B0[temp]/Bmix2 - Log[z - Bmix2] - 
    Amix2/Bmix2 (2. (pres A0[temp]/Amix2)^0.5 - pres B0[temp]/Bmix2) Log[1 + Bmix2/z]
   ]
  ]

And this is how it is used:

Plot[fugacitymix0[273 + 200, x*100000, yfra], {x, 1, 300}] // AbsoluteTiming

It performs the task in about 0.12 seconds.