text: /*---------------------------------------------------------------------------*\
  Copyright (C) 2011-2015 OpenFOAM Foundation
  Copyright (C) 2014 Applied CCM
  -------------------------------------------------------------------------------
  License
      This file is part of CAELUS.
  
      CAELUS is free software: you can redistribute it and/or modify it
      under the terms of the GNU General Public License as published by
      the Free Software Foundation, either version 3 of the License, or
      (at your option) any later version.
  
      CAELUS is distributed in the hope that it will be useful, but WITHOUT
      ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
      for more details.
  
      You should have received a copy of the GNU General Public License
      along with CAELUS.  If not, see <http://www.gnu.org/licenses/>.
  
  \*---------------------------------------------------------------------------*/
  
  #include "compressibleRASModel.hpp"
  #include "wallFvPatch.hpp"
  #include "addToRunTimeSelectionTable.hpp"
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  
  #ifdef windows
  // Symbol to force loading at runtime
  extern "C"
  void compressibleRASModelsLoad()
  {}
  #endif 
  
  namespace CML
  {
  namespace compressible
  {
  
  // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  
  defineTypeNameAndDebug(RASModel, 0);
  defineRunTimeSelectionTable(RASModel, dictionary);
  addToRunTimeSelectionTable(turbulenceModel, RASModel, turbulenceModel);
  
  // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
  
  void RASModel::printCoeffs()
  {
      if (printCoeffs_)
      {
          Info<< type() << "Coeffs" << coeffDict_ << endl;
      }
  }
  
  
  // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  
  RASModel::RASModel
  (
      const word& type,
      const volScalarField& rho,
      const volVectorField& U,
      const surfaceScalarField& phi,
      const fluidThermo& thermophysicalModel,
      const word& turbulenceModelName
  )
  :
      turbulenceModel(rho, U, phi, thermophysicalModel, turbulenceModelName),
  
      IOdictionary
      (
          IOobject
          (
              "RASProperties",
              U.time().constant(),
              U.db(),
              IOobject::MUST_READ_IF_MODIFIED,
              IOobject::NO_WRITE
          )
      ),
  
      turbulence_(lookup("turbulence")),
      printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
      coeffDict_(subOrEmptyDict(type + "Coeffs")),
  
      kMin_("kMin", sqr(dimVelocity), SMALL),
      epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
      omegaMin_("omegaMin", dimless/dimTime, SMALL),
  
      y_(mesh_)
  {
      kMin_.readIfPresent(*this);
      epsilonMin_.readIfPresent(*this);
      omegaMin_.readIfPresent(*this);
  
      // Force the construction of the mesh deltaCoeffs which may be needed
      // for the construction of the derived models and BCs
      mesh_.deltaCoeffs();
  }
  
  
  // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
  
  autoPtr<RASModel> RASModel::New
  (
      const volScalarField& rho,
      const volVectorField& U,
      const surfaceScalarField& phi,
      const fluidThermo& thermophysicalModel,
      const word& turbulenceModelName
  )
  {
      // get model name, but do not register the dictionary
      // otherwise it is registered in the database twice
      const word modelType
      (
          IOdictionary
          (
              IOobject
              (
                  "RASProperties",
                  U.time().constant(),
                  U.db(),
                  IOobject::MUST_READ_IF_MODIFIED,
                  IOobject::NO_WRITE,
                  false
              )
          ).lookup("RASModel")
      );
  
      Info<< "Selecting RAS turbulence model " << modelType << endl;
  
      dictionaryConstructorTable::iterator cstrIter =
          dictionaryConstructorTablePtr_->find(modelType);
  
      if (cstrIter == dictionaryConstructorTablePtr_->end())
      {
          FatalErrorInFunction
              << "Unknown RASModel type "
              << modelType << nl << nl
              << "Valid RASModel types:" << endl
              << dictionaryConstructorTablePtr_->sortedToc()
              << exit(FatalError);
      }
  
      return autoPtr<RASModel>
      (
          cstrIter()(rho, U, phi, thermophysicalModel, turbulenceModelName)
      );
  }
  
  
  // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  
  scalar RASModel::yPlusLam(const scalar kappa, const scalar E) const
  {
      scalar ypl = 11.0;
  
      for (int i=0; i<10; i++)
      {
          ypl = log(max(E*ypl, 1))/kappa;
      }
  
      return ypl;
  }
  
  
  tmp<scalarField> RASModel::yPlus(const label patchNo, const scalar Cmu) const
  {
      const fvPatch& curPatch = mesh_.boundary()[patchNo];
  
      tmp<scalarField> tYp(new scalarField(curPatch.size()));
      scalarField& Yp = tYp();
  
      if (isA<wallFvPatch>(curPatch))
      {
          Yp = pow025(Cmu)
              *y_[patchNo]
              *sqrt(k()().boundaryField()[patchNo].patchInternalField())
             /(
                  mu().boundaryField()[patchNo].patchInternalField()
                 /rho_.boundaryField()[patchNo]
              );
      }
      else
      {
          WarningInFunction
              << "Patch " << patchNo << " is not a wall. Returning null field"
              << nl << endl;
  
          Yp.setSize(0);
      }
  
      return tYp;
  }
  
  
  void RASModel::correct()
  {
      if (mesh_.changing())
      {
          y_.correct();
      }
  }
  
  
  bool RASModel::read()
  {
      //if (regIOobject::read())
  
      // Bit of trickery : we are both IOdictionary ('RASProperties') and
      // an regIOobject from the turbulenceModel level. Problem is to distinguish
      // between the two - we only want to reread the IOdictionary.
  
      bool ok = IOdictionary::readData
      (
          IOdictionary::readStream
          (
              IOdictionary::type()
          )
      );
      IOdictionary::close();
  
      if (ok)
      {
          lookup("turbulence") >> turbulence_;
  
          if (const dictionary* dictPtr = subDictPtr(type() + "Coeffs"))
          {
              coeffDict_ <<= *dictPtr;
          }
  
          kMin_.readIfPresent(*this);
          epsilonMin_.readIfPresent(*this);
          omegaMin_.readIfPresent(*this);
  
          return true;
      }
      else
      {
          return false;
      }
  }
  
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  
  } // End namespace compressible
  } // End namespace CML
  
  // ************************************************************************* //
  