applications/test/Polynomial/Test-Polynomial.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     PolynomialTest
  26
  27 Description
  28     Test application for the templated Polynomial class
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "IStringStream.H"
  33 #include "Polynomial.H"
  34 #include "polynomialFunction.H"
  35 #include "Random.H"
  36 #include "cpuTime.H"
  37
  38 using namespace Foam;
  39
  40 const int PolySize = 8;
  41 const scalar coeff[] = { 0.11, 0.45, -0.94, 1.58, -2.58, 0.08, 3.15, -4.78 };
  42 const char* polyDef = "(0.11 0.45 -0.94 1.58 -2.58 0.08 3.15 -4.78)";
  43
  44
  45 scalar polyValue(const scalar x)
  46 {
  47     // Hard-coded polynomial 8 coeff (7th order)
  48     return
  49         coeff[0]
  50       + coeff[1]*x
  51       + coeff[2]*sqr(x)
  52       + coeff[3]*pow3(x)
  53       + coeff[4]*pow4(x)
  54       + coeff[5]*pow5(x)
  55       + coeff[6]*pow6(x)
  56       + coeff[7]*x*pow6(x);
  57 }
  58
  59
  60 scalar intPolyValue(const scalar x)
  61 {
  62     // Hard-coded integral form of above polynomial
  63     return
  64         coeff[0]*x
  65       + coeff[1]/2.0*sqr(x)
  66       + coeff[2]/3.0*pow3(x)
  67       + coeff[3]/4.0*pow4(x)
  68       + coeff[4]/5.0*pow5(x)
  69       + coeff[5]/6.0*pow6(x)
  70       + coeff[6]/7.0*x*pow6(x)
  71       + coeff[7]/8.0*x*x*pow6(x);
  72 }
  73
  74
  75 scalar polyValue1(const scalar x)
  76 {
  77     // "normal" evaluation using pow()
  78     scalar value = coeff[0];
  79
  80     for (int i=1; i < PolySize; ++i)
  81     {
  82         value += coeff[i]*pow(x, i);
  83     }
  84
  85     return value;
  86 }
  87
  88
  89 scalar polyValue2(const scalar x)
  90 {
  91     // calculation avoiding pow()
  92     scalar value = coeff[0];
  93
  94     scalar powX = x;
  95     for (int i=1; i < PolySize; ++i)
  96     {
  97         value += coeff[i] * powX;
  98         powX *= x;
  99     }
 100
 101     return value;
 102 }
 103
 104
 105 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 106
 107 int main(int argc, char *argv[])
 108 {
 109     const label n = 10000;
 110     const label nIters = 1000;
 111     scalar sum = 0.0;
 112
 113     Info<< "null poly = " << (Polynomial<8>()) << nl
 114         << "null poly = " << (polynomialFunction(8)) << nl
 115         << endl;
 116
 117     Polynomial<8> poly(coeff);
 118     Polynomial<9> intPoly(poly.integral(0.0));
 119
 120     IStringStream is(polyDef);
 121     polynomialFunction polyfunc(is);
 122
 123     Info<< "poly = " << poly << nl
 124         << "intPoly = " << intPoly << nl
 125         << endl;
 126
 127     Info<< "2*poly = " << 2*poly << nl
 128         << "poly+poly = " << poly + poly << nl
 129         << "3*poly = " << 3*poly << nl
 130         << "poly+poly+poly = " << poly + poly + poly << nl
 131         << "3*poly - 2*poly = " << 3*poly - 2*poly << nl
 132         << endl;
 133
 134     Info<< nl << "--- as polynomialFunction" << nl << endl;
 135     Info<< "polyf = " << polyfunc << nl
 136         << "intPoly = " << poly.integral(0.0) << nl
 137         << endl;
 138
 139     Info<< "2*polyf = " << 2*polyfunc << nl
 140         << "polyf+polyf = " << polyfunc + polyfunc << nl
 141         << "3*polyf = " << 3*polyfunc << nl
 142         << "polyf+polyf+polyf = " << polyfunc + polyfunc + polyfunc << nl
 143         << "3*polyf - 2*polyf = " << 3*polyfunc - 2*polyfunc << nl
 144         << endl;
 145
 146
 147     Polynomial<8> polyCopy = poly;
 148     Info<< "poly, polyCopy = " << poly << ", " << polyCopy << nl << endl;
 149     polyCopy = 2.5*poly;
 150     Info<< "2.5*poly = " << polyCopy << nl << endl;
 151
 152     Random rnd(123456);
 153     for (int i=0; i<10; i++)
 154     {
 155         scalar x = rnd.scalar01()*100;
 156
 157         scalar px = polyValue(x);
 158         scalar ipx = intPolyValue(x);
 159
 160         scalar pxTest = poly.value(x);
 161         scalar ipxTest = intPoly.value(x);
 162
 163         Info<<"\nx = " << x << endl;
 164         Info<< "    px, pxTest = " << px << ", " << pxTest << endl;
 165         Info<< "    ipx, ipxTest = " << ipx << ", " << ipxTest << endl;
 166
 167         if (mag(px - pxTest) > SMALL)
 168         {
 169             Info<< "    *** WARNING: px != pxTest: " << px - pxTest << endl;
 170         }
 171
 172         if (mag(ipx - ipxTest) > SMALL)
 173         {
 174             Info<< "    *** WARNING: ipx != ipxTest: " << ipx - ipxTest << endl;
 175         }
 176
 177         Info<< endl;
 178     }
 179
 180
 181     //
 182     // test speed of Polynomial:
 183     //
 184     Info<< "start timing loops" << nl
 185         << "~~~~~~~~~~~~~~~~~~" << endl;
 186
 187     cpuTime timer;
 188
 189     for (int loop = 0; loop < n; ++loop)
 190     {
 191         sum = 0.0;
 192         for (label iter = 0; iter < nIters; ++iter)
 193         {
 194             sum += poly.value(loop+iter);
 195         }
 196     }
 197     Info<< "value:        " << sum
 198         << " in " << timer.cpuTimeIncrement() << " s\n";
 199
 200     for (int loop = 0; loop < n; ++loop)
 201     {
 202         sum = 0.0;
 203         for (label iter = 0; iter < nIters; ++iter)
 204         {
 205             sum += polyfunc.value(loop+iter);
 206         }
 207     }
 208     Info<< "via function: " << sum
 209         << " in " << timer.cpuTimeIncrement() << " s\n";
 210
 211
 212     for (int loop = 0; loop < n; ++loop)
 213     {
 214         sum = 0.0;
 215         for (label iter = 0; iter < nIters; ++iter)
 216         {
 217             sum += polyValue(loop+iter);
 218         }
 219     }
 220     Info<< "hard-coded 0: " << sum
 221         << " in " << timer.cpuTimeIncrement() << " s\n";
 222
 223
 224     for (int loop = 0; loop < n; ++loop)
 225     {
 226         sum = 0.0;
 227         for (label iter = 0; iter < nIters; ++iter)
 228         {
 229             sum += polyValue1(loop+iter);
 230         }
 231     }
 232     Info<< "hard-coded 1: " << sum
 233         << " in " << timer.cpuTimeIncrement() << " s\n";
 234
 235     for (int loop = 0; loop < n; ++loop)
 236     {
 237         sum = 0.0;
 238         for (label iter = 0; iter < nIters; ++iter)
 239         {
 240             sum += polyValue2(loop+iter);
 241         }
 242     }
 243     Info<< "hard-coded 2: " << sum
 244         << " in " << timer.cpuTimeIncrement() << " s\n";
 245
 246
 247     Info<< nl << "Done." << endl;
 248
 249     return 0;
 250 }
 251
 252
 253 // ************************************************************************* //