src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

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