# Suggestions to increase the speed of a fugacity calculation

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.

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 seems faster still to memoize alpha,A,B, and also to make the first arg to fugacitymix a machine double instead of integer. – Daniel Lichtblau Apr 26 '18 at 14:23
• I added pres_?NumericQ to the argument pattern of fugacitymix0. This should inhibit the symbolic computation that causes the error. – Henrik Schumacher Apr 27 '18 at 10:22