src/thermophysicalModels/chemistryModel/chemPoint/chemPoint.C

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