src/thermophysicalModels/chemistryModel/chemPoint/chemPoint.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 #include "IOstream.H"
  32 #include "dictionary.H"
  33 #include "Switch.H"
  34 #include "scalarField.H"
  35 #include "chemPoint.H"
  36 #include "binaryNode.H"
  37
  38 // #define debugEOA
  39 // #define debugGrow
  40 // #define debugCheckSolution
  41
  42 namespace Foam
  43 {
  44
  45 chemPoint::chemPoint
  46 (
  47     const scalarField& v0,
  48     const scalarField& r,
  49     const scalarField& tolerances,
  50     const scalarField& tolerancesSolutions,
  51     const scalar& absErr,
  52     const Switch& logT,
  53     const scalar& deltaT
  54 )
  55 :
  56     node_(NULL),
  57     nUsed_(0),
  58     v0_(v0),
  59     r_(r),
  60     EOA_(tolerances.size()),
  61     solutionsEOA_(tolerancesSolutions.size()),
  62     absErr_(absErr),
  63     logT_(logT),
  64     deltaT_(deltaT),
  65     timeEOA_(0)
  66
  67 {
  68
  69
  70     forAll(tolerances, i)
  71     {
  72         if(tolerances[i] == 0)
  73         {
  74             EOA_[i] = GREAT;
  75         }
  76         else
  77         {
  78             EOA_[i] = absErr * tolerances[i];
  79         }
  80     }
  81
  82     forAll(tolerancesSolutions, i)
  83     {
  84         if(tolerances[i] == 0)
  85         {
  86             solutionsEOA_[i] = GREAT;
  87         }
  88         else
  89         {
  90             solutionsEOA_[i] = absErr * tolerancesSolutions[i];
  91         }
  92     }
  93
  94 }
  95
  96
  97 chemPoint::chemPoint
  98 (
  99     const scalarField& v0,
 100     const scalarField& r,
 101     const scalarField& tolerances,
 102     const scalarField& tolerancesSolutions,
 103     const scalar& absErr,
 104     const Switch& logT,
 105     const scalar& deltaT,
 106     binaryNode* node
 107 )
 108 :
 109     node_(node),
 110     nUsed_(0),
 111     v0_(v0),
 112     r_(r),
 113     EOA_(tolerances.size()),
 114     solutionsEOA_(tolerancesSolutions.size()),
 115     absErr_(absErr),
 116     logT_(logT),
 117     deltaT_(deltaT),
 118     timeEOA_(0)
 119
 120 {
 121     forAll(tolerances, i)
 122     {
 123         if(tolerances[i] == 0)
 124         {
 125             EOA_[i] = 0.0;
 126         }
 127         else
 128         {
 129             EOA_[i] = absErr * tolerances[i];
 130         }
 131     }
 132
 133     forAll(tolerancesSolutions, i)
 134     {
 135         if(tolerancesSolutions[i] == 0)
 136         {
 137             solutionsEOA_[i] = 0.0;
 138         }
 139         else
 140         {
 141             solutionsEOA_[i] = absErr * tolerancesSolutions[i];
 142         }
 143     }
 144
 145
 146 }
 147
 148 chemPoint::chemPoint
 149 (
 150     const chemPoint& p,
 151     binaryNode* node
 152 )
 153 :
 154     node_(node),
 155     nUsed_(0),
 156     v0_(p.v0()),
 157     r_(p.r()),
 158     EOA_(p.EOA()),
 159     solutionsEOA_(p.solutionsEOA()),
 160     absErr_(p.absErr()),
 161     logT_(p.logT()),
 162     deltaT_(p.deltaT()),
 163     timeEOA_(p.timeEOA())
 164 {}
 165
 166 chemPoint::chemPoint
 167 (
 168     const chemPoint& p
 169 )
 170 :
 171     node_(NULL),
 172     nUsed_(0),
 173     v0_(p.v0()),
 174     r_(p.r()),
 175     EOA_(p.EOA()),
 176     solutionsEOA_(p.solutionsEOA()),
 177     absErr_(p.absErr()),
 178     logT_(p.logT()),
 179     deltaT_(p.deltaT()),
 180     timeEOA_(p.timeEOA())
 181 {}
 182
 183 bool chemPoint::inEOA
 184 (
 185     const scalarField& point,
 186     const scalarField& Wi,
 187     const scalar& rhoi,
 188     const scalar& deltaT,
 189     const scalarField& scaleFactor
 190 )
 191 {
 192
 193     if(nUsed_ < INT_MAX)
 194     {
 195         nUsed_++;
 196     }
 197
 198     label size = point.size();
 199
 200     scalar eps2 = 0.0;
 201
 202     for(label i=0; i < size; i++)
 203     {
 204
 205 //      Temperature
 206         if(i == size-2 && logT_)
 207         {
 208             scalar diffAxis = sqr(Foam::scalar(mag(log(point[i])-log(v0_[i]))))/sqr(EOA_[i]) ;
 209             eps2 += diffAxis;
 210         }
 211         else if(i == size-2 && !logT_)
 212         {
 213             scalar diffAxis = sqr(Foam::scalar(mag((point[i])-(v0_[i])))) / sqr(EOA_[i]) ;
 214             eps2 += diffAxis;
 215         }
 216         if(i == size-1)
 217         {
 218             scalar diffAxis = sqr(Foam::scalar(mag((point[i]-v0_[i])))) / sqr(EOA_[i]) ;
 219             eps2 += diffAxis;
 220         }
 221 //      Species
 222         else
 223         {
 224             scalar diffAxis = sqr(Foam::scalar(mag((point[i]-v0_[i])))) / sqr(EOA_[i]) ;
 225             eps2 += diffAxis;
 226         }
 227     }
 228
 229
 230     if(eps2 > 1.0)
 231     {
 232         return false;
 233     }
 234     else
 235     {
 236         return true;
 237     }
 238
 239 }
 240
 241 void chemPoint::grow
 242 (
 243     const scalarField& v0,
 244     const scalarField& Wi,
 245     const scalar& rhoi,
 246     const scalarField& v,
 247     const scalar& deltaT
 248 )
 249 {
 250
 251     label size = v.size();
 252
 253     for(label i=0; i < size-2; i++)
 254     {
 255         if(i == size-2 && logT_)
 256         {
 257             scalar diffAxis = Foam::scalar(mag(log(v0[i])-log(v0_[i]))) - EOA_[i] ;
 258
 259             if(diffAxis > Foam::VSMALL)
 260             {
 261                 EOA_[i] = mag(log(v0[i])-log(v0_[i]));
 262             }
 263         }
 264         else if(i == size-2 && !logT_)
 265         {
 266             scalar diffAxis = Foam::scalar(mag((v0[i])-(v0_[i]))) - EOA_[i] ;
 267
 268             if(diffAxis > Foam::VSMALL)
 269             {
 270                 EOA_[i] = mag((v0[i])-(v0_[i]));
 271             }
 272         }
 273         else
 274         {
 275             scalar diffAxis = Foam::scalar(mag((v0[i]-v0_[i]))) - EOA_[i] ;
 276
 277             if(diffAxis > Foam::VSMALL)
 278             {
 279                 EOA_[i] = mag((v0[i]-v0_[i]));
 280             }
 281
 282         }
 283     }
 284
 285     timeEOA_ = Foam::scalar(mag(deltaT_ - deltaT));
 286
 287 }
 288
 289 bool chemPoint::checkSolution
 290 (
 291     const scalarField& v0,
 292     const scalarField& v,
 293     const scalarField& Wi,
 294     const scalar& rhoi,
 295     const scalar& T,
 296     const scalar& p,
 297     const scalar& deltaT,
 298     const scalarField& tolerances
 299 )
 300 {
 301
 302     label size = v.size();
 303
 304     scalarField s = v;
 305
 306     s.setSize(size+2);
 307
 308 //    s[size-2] = T;
 309 //    s[size-1] = p;
 310
 311     s[s.size()-2] = T;
 312     s[s.size()-1] = p;
 313     size = s.size();
 314
 315     scalar eps2 = 0.0;
 316
 317
 318     for(label i=0; i < size; i++)
 319     {
 320         if(i == size-2 && logT_)
 321         {
 322             scalar diffAxis = (Foam::scalar(mag(log(s[i])-log(r_[i])))) / (solutionsEOA_[i]) ;
 323             eps2 += sqr(diffAxis);
 324         }
 325         else if(i == size-2 && !logT_)
 326         {
 327             scalar diffAxis = (Foam::scalar(mag(s[i]-r_[i]))) / (solutionsEOA_[i]) ;
 328             eps2 += sqr(diffAxis);
 329         }
 330         else if(i == size-1)
 331         {
 332             scalar diffAxis = (Foam::scalar(mag(s[i]-r_[i]))) / (solutionsEOA_[i]) ;
 333             eps2 += sqr(diffAxis);
 334         }
 335         else
 336         {
 337             scalar diffAxis = (Foam::scalar(mag((s[i]*Wi[i]/rhoi-r_[i])))) / (solutionsEOA_[i]) ;
 338             eps2 += sqr(diffAxis);
 339         }
 340     }
 341
 342     if(eps2 > 1.0)
 343     {
 344         return false;
 345     }
 346     else
 347     {
 348         // if the solution is in the ellipsoid of accuracy, grow it!
 349         grow(v0, Wi, rhoi, s, deltaT);
 350         return true;
 351     }
 352
 353
 354 }
 355
 356 void chemPoint::setFree()
 357 {
 358     node_ = NULL;
 359     nUsed_ = 0;
 360     /*
 361         v0_.clear();
 362         r_.clear();
 363         EOA_.clear();
 364         solutionsEOA_.clear();
 365         absErr_ = 0;
 366         logT_ = 0;
 367         deltaT_ = 0;
 368         timeEOA_ = 0;
 369     */
 370 }
 371
 372 void chemPoint::clearData()
 373 {
 374         nUsed_ = 0;
 375 //        Info << "nUsed" << endl;
 376         v0_.clear();
 377 //        Info << "v0.clear" << endl;
 378         r_.clear();
 379 //        Info << "r_.clear()" << endl;
 380         EOA_.clear();
 381 //        Info << "EOA_.clear()" << endl;
 382         solutionsEOA_.clear();
 383 //        Info << "solutionsEOA_.clear()" << endl;
 384         absErr_ = 0;
 385 //        Info << "absErr_ = 0" << endl;
 386 //        logT_ = 0;
 387 //        Info << "logT_ = 0" << endl;
 388         deltaT_ = 0;
 389 //        Info << "deltaT_ = 0;" << endl;
 390         timeEOA_ = 0;
 391 //        node_ = NULL;
 392 //        Info << "chemPoint::clearData():: ENDDDDDD!" << endl;
 393
 394 }
 395
 396
 397 }