src/dbns/timeStepping/EulerLocalDdtScheme/EulerLocalDdtScheme.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 Class
  25     Foam::fv::EulerLocalDdtScheme
  26
  27 Author
  28     Oliver Borm
  29     Aleksandar Jemcov
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "EulerLocalDdtScheme.H"
  34 #include "surfaceInterpolate.H"
  35 #include "fvcDiv.H"
  36 #include "fvMatrices.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 namespace fv
  46 {
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 template<class Type>
  51 tmp<GeometricField<Type, fvPatchField, volMesh> >
  52 EulerLocalDdtScheme<Type>::fvcDdt
  53 (
  54     const dimensioned<Type>& dt
  55 )
  56 {
  57     const objectRegistry& registry = this->mesh();
  58
  59     // get access to the scalar beta[i]
  60     const scalarField& beta =
  61         registry.lookupObject<scalarField>(deltaTName_);
  62
  63     volScalarField rDeltaT =
  64         1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
  65
  66     IOobject ddtIOobject
  67     (
  68         "ddt("+dt.name()+')',
  69         mesh().time().timeName(),
  70         mesh()
  71     );
  72
  73     if (mesh().moving())
  74     {
  75         tmp<GeometricField<Type, fvPatchField, volMesh> > tdtdt
  76         (
  77             new GeometricField<Type, fvPatchField, volMesh>
  78             (
  79                 ddtIOobject,
  80                 mesh(),
  81                 dimensioned<Type>
  82                 (
  83                     "0",
  84                     dt.dimensions()/dimTime,
  85                     pTraits<Type>::zero
  86                 )
  87             )
  88         );
  89
  90         tdtdt().internalField() =
  91             rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
  92
  93         return tdtdt;
  94     }
  95     else
  96     {
  97         return tmp<GeometricField<Type, fvPatchField, volMesh> >
  98         (
  99             new GeometricField<Type, fvPatchField, volMesh>
 100             (
 101                 ddtIOobject,
 102                 mesh(),
 103                 dimensioned<Type>
 104                 (
 105                     "0",
 106                     dt.dimensions()/dimTime,
 107                     pTraits<Type>::zero
 108                 ),
 109                 calculatedFvPatchField<Type>::typeName
 110             )
 111         );
 112     }
 113 }
 114
 115
 116 template<class Type>
 117 tmp<GeometricField<Type, fvPatchField, volMesh> >
 118 EulerLocalDdtScheme<Type>::fvcDdt
 119 (
 120     const GeometricField<Type, fvPatchField, volMesh>& vf
 121 )
 122 {
 123     const objectRegistry& registry = this->mesh();
 124
 125     // get access to the scalar beta[i]
 126     const scalarField& beta =
 127         registry.lookupObject<scalarField>(deltaTName_);
 128
 129     volScalarField rDeltaT =
 130         1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
 131
 132     IOobject ddtIOobject
 133     (
 134         "ddt("+vf.name()+')',
 135         mesh().time().timeName(),
 136         mesh()
 137     );
 138
 139     if (mesh().moving())
 140     {
 141         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 142         (
 143             new GeometricField<Type, fvPatchField, volMesh>
 144             (
 145                 ddtIOobject,
 146                 mesh(),
 147                 rDeltaT.dimensions()*vf.dimensions(),
 148                 rDeltaT.internalField()*
 149                 (
 150                     vf.internalField()
 151                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 152                 ),
 153                 rDeltaT.boundaryField()*
 154                 (
 155                     vf.boundaryField() - vf.oldTime().boundaryField()
 156                 )
 157             )
 158         );
 159     }
 160     else
 161     {
 162         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 163         (
 164             new GeometricField<Type, fvPatchField, volMesh>
 165             (
 166                 ddtIOobject,
 167                 rDeltaT*(vf - vf.oldTime())
 168             )
 169         );
 170     }
 171 }
 172
 173
 174 template<class Type>
 175 tmp<GeometricField<Type, fvPatchField, volMesh> >
 176 EulerLocalDdtScheme<Type>::fvcDdt
 177 (
 178     const dimensionedScalar& rho,
 179     const GeometricField<Type, fvPatchField, volMesh>& vf
 180 )
 181 {
 182     const objectRegistry& registry = this->mesh();
 183
 184     // get access to the scalar beta[i]
 185     const scalarField& beta =
 186         registry.lookupObject<scalarField>(deltaTName_);
 187
 188     volScalarField rDeltaT =
 189         1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
 190
 191     IOobject ddtIOobject
 192     (
 193         "ddt("+rho.name()+','+vf.name()+')',
 194         mesh().time().timeName(),
 195         mesh()
 196     );
 197
 198     if (mesh().moving())
 199     {
 200         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 201         (
 202             new GeometricField<Type, fvPatchField, volMesh>
 203             (
 204                 ddtIOobject,
 205                 mesh(),
 206                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 207                 rDeltaT.internalField()*rho.value()*
 208                 (
 209                     vf.internalField()
 210                   - vf.oldTime().internalField()*mesh().V0()/mesh().V()
 211                 ),
 212                 rDeltaT.boundaryField()*rho.value()*
 213                 (
 214                     vf.boundaryField() - vf.oldTime().boundaryField()
 215                 )
 216             )
 217         );
 218     }
 219     else
 220     {
 221         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 222         (
 223             new GeometricField<Type, fvPatchField, volMesh>
 224             (
 225                 ddtIOobject,
 226                 rDeltaT*rho*(vf - vf.oldTime())
 227             )
 228         );
 229     }
 230 }
 231
 232
 233 template<class Type>
 234 tmp<GeometricField<Type, fvPatchField, volMesh> >
 235 EulerLocalDdtScheme<Type>::fvcDdt
 236 (
 237     const volScalarField& rho,
 238     const GeometricField<Type, fvPatchField, volMesh>& vf
 239 )
 240 {
 241     const objectRegistry& registry = this->mesh();
 242
 243     // get access to the scalar beta[i]
 244     const scalarField& beta =
 245         registry.lookupObject<scalarField>(deltaTName_);
 246
 247     volScalarField rDeltaT =
 248         1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
 249
 250     IOobject ddtIOobject
 251     (
 252         "ddt("+rho.name()+','+vf.name()+')',
 253         mesh().time().timeName(),
 254         mesh()
 255     );
 256
 257     if (mesh().moving())
 258     {
 259         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 260         (
 261             new GeometricField<Type, fvPatchField, volMesh>
 262             (
 263                 ddtIOobject,
 264                 mesh(),
 265                 rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
 266                 rDeltaT.internalField()*
 267                 (
 268                     rho.internalField()*vf.internalField()
 269                   - rho.oldTime().internalField()
 270                    *vf.oldTime().internalField()*mesh().V0()/mesh().V()
 271                 ),
 272                 rDeltaT.boundaryField()*
 273                 (
 274                     rho.boundaryField()*vf.boundaryField()
 275                   - rho.oldTime().boundaryField()
 276                    *vf.oldTime().boundaryField()
 277                 )
 278             )
 279         );
 280     }
 281     else
 282     {
 283         return tmp<GeometricField<Type, fvPatchField, volMesh> >
 284         (
 285             new GeometricField<Type, fvPatchField, volMesh>
 286             (
 287                 ddtIOobject,
 288                 rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
 289             )
 290         );
 291     }
 292 }
 293
 294
 295 template<class Type>
 296 tmp<fvMatrix<Type> >
 297 EulerLocalDdtScheme<Type>::fvmDdt
 298 (
 299     GeometricField<Type, fvPatchField, volMesh>& vf
 300 )
 301 {
 302     const objectRegistry& registry = this->mesh();
 303
 304     // get access to the scalar beta[i]
 305     const scalarField& beta =
 306         registry.lookupObject<scalarField>(deltaTName_);
 307
 308     tmp<fvMatrix<Type> > tfvm
 309     (
 310         new fvMatrix<Type>
 311         (
 312             vf,
 313             vf.dimensions()*dimVol/dimTime
 314         )
 315     );
 316
 317     fvMatrix<Type>& fvm = tfvm();
 318
 319     scalarField rDeltaT =
 320         1.0/(beta[0]*registry.lookupObject<volScalarField>
 321         (
 322             deltaTauName_
 323         ).internalField());
 324
 325     fvm.diag() = rDeltaT*mesh().V();
 326
 327     if (mesh().moving())
 328     {
 329         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
 330     }
 331     else
 332     {
 333         fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
 334     }
 335
 336     return tfvm;
 337 }
 338
 339
 340 template<class Type>
 341 tmp<fvMatrix<Type> >
 342 EulerLocalDdtScheme<Type>::fvmDdt
 343 (
 344     const dimensionedScalar& rho,
 345     GeometricField<Type, fvPatchField, volMesh>& vf
 346 )
 347 {
 348     const objectRegistry& registry = this->mesh();
 349
 350     // get access to the scalar beta[i]
 351     const scalarField& beta =
 352         registry.lookupObject<scalarField>(deltaTName_);
 353
 354     tmp<fvMatrix<Type> > tfvm
 355     (
 356         new fvMatrix<Type>
 357         (
 358             vf,
 359             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 360         )
 361     );
 362     fvMatrix<Type>& fvm = tfvm();
 363
 364     scalarField rDeltaT =
 365         1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_).internalField());
 366
 367     fvm.diag() = rDeltaT*rho.value()*mesh().V();
 368
 369     if (mesh().moving())
 370     {
 371         fvm.source() = rDeltaT
 372             *rho.value()*vf.oldTime().internalField()*mesh().V0();
 373     }
 374     else
 375     {
 376         fvm.source() = rDeltaT
 377             *rho.value()*vf.oldTime().internalField()*mesh().V();
 378     }
 379
 380     return tfvm;
 381 }
 382
 383
 384 template<class Type>
 385 tmp<fvMatrix<Type> >
 386 EulerLocalDdtScheme<Type>::fvmDdt
 387 (
 388     const volScalarField& rho,
 389     GeometricField<Type, fvPatchField, volMesh>& vf
 390 )
 391 {
 392     const objectRegistry& registry = this->mesh();
 393
 394     // get access to the scalar beta[i]
 395     const scalarField& beta =
 396         registry.lookupObject<scalarField>(deltaTName_);
 397
 398     tmp<fvMatrix<Type> > tfvm
 399     (
 400         new fvMatrix<Type>
 401         (
 402             vf,
 403             rho.dimensions()*vf.dimensions()*dimVol/dimTime
 404         )
 405     );
 406     fvMatrix<Type>& fvm = tfvm();
 407
 408     scalarField rDeltaT =
 409         1.0/(beta[0]*registry.lookupObject<volScalarField>
 410         (
 411             deltaTauName_
 412         ).internalField());
 413
 414     fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
 415
 416     if (mesh().moving())
 417     {
 418         fvm.source() = rDeltaT
 419             *rho.oldTime().internalField()
 420             *vf.oldTime().internalField()*mesh().V0();
 421     }
 422     else
 423     {
 424         fvm.source() = rDeltaT
 425             *rho.oldTime().internalField()
 426             *vf.oldTime().internalField()*mesh().V();
 427     }
 428
 429     return tfvm;
 430 }
 431
 432
 433 template<class Type>
 434 tmp<typename EulerLocalDdtScheme<Type>::fluxFieldType>
 435 EulerLocalDdtScheme<Type>::fvcDdtPhiCorr
 436 (
 437     const volScalarField& rA,
 438     const GeometricField<Type, fvPatchField, volMesh>& U,
 439     const fluxFieldType& phi
 440 )
 441 {
 442     IOobject ddtIOobject
 443     (
 444         "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
 445         mesh().time().timeName(),
 446         mesh()
 447     );
 448
 449     if (mesh().moving())
 450     {
 451         return tmp<fluxFieldType>
 452         (
 453             new fluxFieldType
 454             (
 455                 ddtIOobject,
 456                 mesh(),
 457                 dimensioned<typename flux<Type>::type>
 458                 (
 459                     "0",
 460                     rA.dimensions()*phi.dimensions()/dimTime,
 461                     pTraits<typename flux<Type>::type>::zero
 462                 )
 463             )
 464         );
 465     }
 466     else
 467     {
 468         const objectRegistry& registry = this->mesh();
 469
 470         // get access to the scalar beta[i]
 471         const scalarField& beta =
 472             registry.lookupObject<scalarField>(deltaTName_);
 473
 474         volScalarField rDeltaT =
 475             1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
 476
 477         return tmp<fluxFieldType>
 478         (
 479             new fluxFieldType
 480             (
 481                 ddtIOobject,
 482                 this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
 483                 (
 484                     fvc::interpolate(rDeltaT*rA)*phi.oldTime()
 485                   - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
 486                 )
 487             )
 488         );
 489     }
 490 }
 491
 492
 493 template<class Type>
 494 tmp<typename EulerLocalDdtScheme<Type>::fluxFieldType>
 495 EulerLocalDdtScheme<Type>::fvcDdtPhiCorr
 496 (
 497     const volScalarField& rA,
 498     const volScalarField& rho,
 499     const GeometricField<Type, fvPatchField, volMesh>& U,
 500     const fluxFieldType& phi
 501 )
 502 {
 503     IOobject ddtIOobject
 504     (
 505         "ddtPhiCorr("
 506       + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
 507         mesh().time().timeName(),
 508         mesh()
 509     );
 510
 511     if (mesh().moving())
 512     {
 513         return tmp<fluxFieldType>
 514         (
 515             new fluxFieldType
 516             (
 517                 ddtIOobject,
 518                 mesh(),
 519                 dimensioned<typename flux<Type>::type>
 520                 (
 521                     "0",
 522                     rA.dimensions()*phi.dimensions()/dimTime,
 523                     pTraits<typename flux<Type>::type>::zero
 524                 )
 525             )
 526         );
 527     }
 528     else
 529     {
 530         const objectRegistry& registry = this->mesh();
 531
 532         // get access to the scalar beta[i]
 533         const scalarField& beta =
 534             registry.lookupObject<scalarField>(deltaTName_);
 535
 536         volScalarField rDeltaT =
 537             1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));
 538
 539         if
 540         (
 541             U.dimensions() == dimVelocity
 542          && phi.dimensions() == dimVelocity*dimArea
 543         )
 544         {
 545             return tmp<fluxFieldType>
 546             (
 547                 new fluxFieldType
 548                 (
 549                     ddtIOobject,
 550                     this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
 551                     (
 552                         fvc::interpolate
 553                         (
 554                             rDeltaT*rA*rho.oldTime()
 555                         )*phi.oldTime()
 556                       - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
 557                       & mesh().Sf())
 558                     )
 559                 )
 560             );
 561         }
 562         else if
 563         (
 564             U.dimensions() == dimVelocity
 565          && phi.dimensions() == dimDensity*dimVelocity*dimArea
 566         )
 567         {
 568             return tmp<fluxFieldType>
 569             (
 570                 new fluxFieldType
 571                 (
 572                     ddtIOobject,
 573                     this->fvcDdtPhiCoeff
 574                     (
 575                         U.oldTime(),
 576                         phi.oldTime()/fvc::interpolate(rho.oldTime())
 577                     )*
 578                     (
 579                         fvc::interpolate(rDeltaT*rA*rho.oldTime())
 580                        *phi.oldTime()/fvc::interpolate(rho.oldTime())
 581                       - (
 582                             fvc::interpolate
 583                             (
 584                                 rDeltaT*rA*rho.oldTime()*U.oldTime()
 585                             ) & mesh().Sf()
 586                         )
 587                     )
 588                 )
 589             );
 590         }
 591         else if
 592         (
 593             U.dimensions() == dimDensity*dimVelocity
 594          && phi.dimensions() == dimDensity*dimVelocity*dimArea
 595         )
 596         {
 597             return tmp<fluxFieldType>
 598             (
 599                 new fluxFieldType
 600                 (
 601                     ddtIOobject,
 602                     this->fvcDdtPhiCoeff(rho.oldTime(), U.oldTime(), phi.oldTime())*
 603                    (
 604                         fvc::interpolate(rDeltaT*rA)*phi.oldTime()
 605                       - (
 606                             fvc::interpolate(rDeltaT*rA*U.oldTime())
 607                           & mesh().Sf()
 608                         )
 609                     )
 610                 )
 611             );
 612         }
 613         else
 614         {
 615             FatalErrorIn
 616             (
 617                 "EulerLocalDdtScheme<Type>::fvcDdtPhiCorr"
 618             )   << "dimensions of phi are not correct"
 619                 << abort(FatalError);
 620
 621             return fluxFieldType::null();
 622         }
 623     }
 624 }
 625
 626
 627 template<class Type>
 628 tmp<surfaceScalarField> EulerLocalDdtScheme<Type>::meshPhi
 629 (
 630     const GeometricField<Type, fvPatchField, volMesh>&
 631 )
 632 {
 633     return tmp<surfaceScalarField>
 634     (
 635         new surfaceScalarField
 636         (
 637             IOobject
 638             (
 639                 "meshPhi",
 640                 mesh().time().timeName(),
 641                 mesh()
 642             ),
 643             mesh(),
 644             dimensionedScalar("0", dimVolume/dimTime, 0.0)
 645         )
 646     );
 647 }
 648
 649
 650 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 651
 652 } // End namespace fv
 653
 654 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 655
 656 } // End namespace Foam
 657
 658 // ************************************************************************* //