From bd7d1a438af3d0ca92b605786cbbc9cd79dc93e5 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 10 May 2023 11:02:18 +0100
Subject: [PATCH] scalar: Added optimised integer powers and roots

---
 applications/test/integerPow/Make/files       |  3 +
 applications/test/integerPow/Make/options     |  0
 .../test/integerPow/Test-integerPow.C         | 87 +++++++++++++++++++
 .../fields/Fields/scalarField/scalarField.C   |  8 +-
 .../fields/Fields/scalarField/scalarField.H   |  9 +-
 .../primitives/Scalar/scalar/scalar.H         | 11 +++
 .../primitives/Scalar/scalar/scalarI.H        | 74 ++++++++++++++++
 .../reaction/specieExponent/specieExponentI.H | 25 +-----
 8 files changed, 192 insertions(+), 25 deletions(-)
 create mode 100644 applications/test/integerPow/Make/files
 create mode 100644 applications/test/integerPow/Make/options
 create mode 100644 applications/test/integerPow/Test-integerPow.C
 create mode 100644 src/OpenFOAM/primitives/Scalar/scalar/scalarI.H

diff --git a/applications/test/integerPow/Make/files b/applications/test/integerPow/Make/files
new file mode 100644
index 0000000000..8e07339c39
--- /dev/null
+++ b/applications/test/integerPow/Make/files
@@ -0,0 +1,3 @@
+Test-integerPow.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-integerPow
diff --git a/applications/test/integerPow/Make/options b/applications/test/integerPow/Make/options
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/applications/test/integerPow/Test-integerPow.C b/applications/test/integerPow/Test-integerPow.C
new file mode 100644
index 0000000000..b6d58654d6
--- /dev/null
+++ b/applications/test/integerPow/Test-integerPow.C
@@ -0,0 +1,87 @@
+#include "cpuTime.H"
+#include "Pair.H"
+#include "Random.H"
+
+using namespace Foam;
+
+//- Number of calculations to do. Need 10^8-ish to get enough time (~2 seconds)
+//  to meaningfully compare.
+static const label n = 100000000;
+
+//- Range of exponents to calculate. From -E to +E. Zero will be omitted so
+//  that we can do roots, too. At small values (<8) integer powers win. At
+//  larger values (>16), scalar powers do better.
+static const label E = 4;
+
+//- Run a given power method
+template<class Pow>
+Pair<scalar> run(Pow pow)
+{
+    Random rndGen(0);
+
+    cpuTime time;
+
+    scalar y = 0;
+
+    for (label i = 0; i < n; ++ i)
+    {
+        const scalar x = rndGen.sampleAB<scalar>(rootSmall, rootGreat);
+        const label e =
+            rndGen.sample01<label>()
+          ? rndGen.sampleAB<label>(- E, 0)
+          : rndGen.sampleAB<label>(1, E + 1);
+
+        y = y*scalar(i)/(i + 1) + pow(x, e)/(i + 1);
+    }
+
+    return Pair<scalar>(y, time.cpuTimeIncrement());
+}
+
+//- Scalar power function
+scalar scalarPow(const scalar x, const scalar e)
+{
+    return Foam::pow(x, e);
+}
+
+//- Scalar root function
+scalar scalarRoot(const scalar x, const scalar e)
+{
+    return Foam::pow(x, 1/e);
+}
+
+//- Run tests
+int main()
+{
+    Info<< "Averaging " << n << " power/root calculations with exponents "
+        << "ranging from " << - E << " to " << + E << ":" << endl;
+
+    Info<< incrIndent;
+
+    Info<< indent << " Scalar power: " << flush;
+    const Pair<scalar> scalarPowYTime =
+        run(static_cast<scalar (*)(scalar, scalar)>(&scalarPow));
+    Info<< "Calculated value " << scalarPowYTime.first()
+        << " in " << scalarPowYTime.second() << " s" << endl;
+
+    Info<< indent << " Integer power: " << flush;
+    const Pair<scalar> integerPowYTime =
+        run(static_cast<scalar (*)(scalar, label)>(&integerPow));
+    Info<< "Calculated value " << integerPowYTime.first()
+        << " in " << integerPowYTime.second() << " s" << endl;
+
+    Info<< indent << " Scalar root: " << flush;
+    const Pair<scalar> scalarRootYTime =
+        run(static_cast<scalar (*)(scalar, scalar)>(&scalarRoot));
+    Info<< "Calculated value " << scalarRootYTime.first()
+        << " in " << scalarRootYTime.second() << " s" << endl;
+
+    Info<< indent << " Integer root: " << flush;
+    const Pair<scalar> integerRootYTime =
+        run(static_cast<scalar (*)(scalar, label)>(&integerRoot));
+    Info<< "Calculated value " << integerRootYTime.first()
+        << " in " << integerRootYTime.second() << " s" << endl;
+
+    Info<< decrIndent;
+
+    return 0;
+}
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
index 2fe51a697b..8da5ffac40 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,6 +50,7 @@ void component(scalarField& sf, const UList<scalar>& f, const direction)
     sf = f;
 }
 
+
 template<>
 void scalarField::replace(const direction, const UList<scalar>& sf)
 {
@@ -117,10 +118,15 @@ BINARY_TYPE_OPERATOR_SF(scalar, scalar, scalar, /, divide)
 
 BINARY_FUNCTION(scalar, scalar, scalar, pow)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, pow)
+BINARY_FUNCTION(scalar, scalar, label, integerPow)
+BINARY_TYPE_FUNCTION_FS(scalar, scalar, label, integerPow)
+BINARY_FUNCTION(scalar, scalar, label, integerRoot)
+BINARY_TYPE_FUNCTION_FS(scalar, scalar, label, integerRoot)
 
 BINARY_FUNCTION(scalar, scalar, scalar, atan2)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, atan2)
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 UNARY_FUNCTION(scalar, scalar, pow3)
diff --git a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
index 4e25f35f48..0d6cace18d 100644
--- a/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
+++ b/src/OpenFOAM/fields/Fields/scalarField/scalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,6 +60,7 @@ void component
     const direction
 );
 
+
 template<>
 void scalarField::replace(const direction, const UList<scalar>& sf);
 
@@ -68,7 +69,9 @@ void scalarField::replace(const direction, const scalar& s);
 
 
 void stabilise(scalarField& Res, const UList<scalar>& sf, const scalar s);
+
 tmp<scalarField> stabilise(const UList<scalar>&, const scalar s);
+
 tmp<scalarField> stabilise(const tmp<scalarField>&, const scalar s);
 
 
@@ -90,6 +93,10 @@ BINARY_TYPE_OPERATOR_SF(scalar, scalar, scalar, /, divide)
 
 BINARY_FUNCTION(scalar, scalar, scalar, pow)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, pow)
+BINARY_FUNCTION(scalar, scalar, label, integerPow)
+BINARY_TYPE_FUNCTION_FS(scalar, scalar, label, integerPow)
+BINARY_FUNCTION(scalar, scalar, label, integerRoot)
+BINARY_TYPE_FUNCTION_FS(scalar, scalar, label, integerRoot)
 
 BINARY_FUNCTION(scalar, scalar, scalar, atan2)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, atan2)
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H b/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
index fdf91ddd20..fdb6cf3037 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
+++ b/src/OpenFOAM/primitives/Scalar/scalar/scalar.H
@@ -124,6 +124,7 @@ namespace Foam
 
 namespace Foam
 {
+    //- Read a single scalar from an input stream
     scalar readScalar(Istream& is);
 
     //- Normalised upper incomplete gamma function
@@ -140,9 +141,19 @@ namespace Foam
 
     //- Inverse normalised lower incomplete gamma function
     scalar invIncGammaRatio_P(const scalar a, const scalar P);
+
+    //- Compute the power of the number x to the integer e
+    inline scalar integerPow(const scalar x, const label e);
+
+    //- Compute the power of the number x to the reciprocal integer 1/e
+    inline scalar integerRoot(const scalar x, const label e);
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "scalarI.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/scalarI.H b/src/OpenFOAM/primitives/Scalar/scalar/scalarI.H
new file mode 100644
index 0000000000..55bff35b61
--- /dev/null
+++ b/src/OpenFOAM/primitives/Scalar/scalar/scalarI.H
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::integerPow(const scalar x, const label e)
+{
+    if (e == 0)
+    {
+        return 1;
+    }
+
+    scalar xx = e > 0 ? x : 1/x;
+
+    scalar y = 1;
+
+    for (label i = e; i != 0; i /= 2)
+    {
+        if (i % 2)
+        {
+            y *= xx;
+        }
+
+        xx *= xx;
+    }
+
+    return y;
+}
+
+
+inline Foam::scalar Foam::integerRoot(const scalar x, const label e)
+{
+    switch (e)
+    {
+        case -4: return 1/pow025(x);
+        case -3: return 1/cbrt(x);
+        case -2: return 1/sqrt(x);
+        case -1: return 1/x;
+        case 0: return NaN;
+        case 1: return x;
+        case 2: return sqrt(x);
+        case 3: return cbrt(x);
+        case 4: return pow025(x);
+    }
+
+    return pow(x, scalar(1)/e);
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H
index eb9068e982..c2c2870717 100644
--- a/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H
+++ b/src/thermophysicalModels/specie/reaction/specieExponent/specieExponentI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,28 +103,7 @@ inline Foam::scalar Foam::pow(const scalar x, const specieExponent& e)
 {
     if (e.hasIntegerExponent())
     {
-        if (e.integerExponent_ == 0)
-        {
-            return 1;
-        }
-        else
-        {
-            scalar xx = e.integerExponent_ > 0 ? x : 1/x;
-
-            scalar y = 1;
-
-            for (label i = e.integerExponent_; i != 0; i /= 2)
-            {
-                if (i % 2)
-                {
-                    y *= xx;
-                }
-
-                xx *= xx;
-            }
-
-            return y;
-        }
+        return integerPow(x, e.integerExponent_);
     }
     else
     {