7
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It seems faster still to memoize alpha,A,B, and also to make the first arg to fugacitymix a machine double instead of integer. $\endgroup$ – Daniel Lichtblau Apr 26 '18 at 14:23
  • $\begingroup$ @DanielLichtblau Thanks for the suggestions. I added them in a new edit. $\endgroup$ – Henrik Schumacher Apr 26 '18 at 15:18
  • $\begingroup$ Definitely a strong improvement over the original. However when I test the command for the plot, I receive a lot of error messages such as "NRoots::nnumeq: -5.72408*10^-8 x^2+(-0.000370503 x-2.37987*10^-7 x^2) z-z^2+z^3==0 is expected to be a polynomial equation in the variable z with numeric coefficients." The task is eventually performed but it looks like there's some issue with evaluation. Is there a way around that ? $\endgroup$ – Whelp Apr 27 '18 at 10:10
  • $\begingroup$ I added pres_?NumericQ to the argument pattern of fugacitymix0. This should inhibit the symbolic computation that causes the error. $\endgroup$ – Henrik Schumacher Apr 27 '18 at 10:22
  • 1
    $\begingroup$ You might be able to get away with compiling the Cardano-Tartaglia formula for the roots. Care must be taken to allow that the result can be complex even though input is real. Also it is not obvious whether this will offer any speed improvement. $\endgroup$ – Daniel Lichtblau Apr 27 '18 at 14:22

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.