.. _program_listing_file_wallModels_ODEWallModelFvPatchScalarField.C:

Program Listing for File ODEWallModelFvPatchScalarField.C
=========================================================

|exhale_lsh| :ref:`Return to documentation for file <file_wallModels_ODEWallModelFvPatchScalarField.C>` (``wallModels/ODEWallModelFvPatchScalarField.C``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /*---------------------------------------------------------------------------* \
   License
       This file is part of libWallModelledLES.
   
       libWallModelledLES 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.
   
       libWallModelledLES 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 libWallModelledLES.
       If not, see <http://www.gnu.org/licenses/>.
   
   \*---------------------------------------------------------------------------*/
   
   #include "ODEWallModelFvPatchScalarField.H"
   #include "addToRunTimeSelectionTable.H"
   #include "codeRules.H"
   #include "scalarListIOList.H"
   #include "helpers.H"
   #include "AdaptiveIntegrator.hpp"
   #include <functional>
   
   typedef std::function<Foam::scalar(const Foam::scalar)> IntegrandFunc;
   
   // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
   
   #if !defined(DOXYGEN_SHOULD_SKIP_THIS)
   namespace Foam
   {
       defineTypeNameAndDebug(ODEWallModelFvPatchScalarField, 0);
   }
   #endif
   
   
   // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
   
   void Foam::ODEWallModelFvPatchScalarField::writeLocalEntries(Ostream& os) const
   {
       wallModelFvPatchScalarField::writeLocalEntries(os);
       eddyViscosity_->write(os);
       sampler_->write(os);
       os.writeKeyword("eps") << eps_ << token::END_STATEMENT << endl;
       os.writeKeyword("maxIter") << maxIter_ << token::END_STATEMENT << endl;
   }
   
   
   Foam::tmp<Foam::scalarField>
   Foam::ODEWallModelFvPatchScalarField::calcNut() const
   {
       if (debug)
       {
           Info<< "Updating nut for patch " << patch().name() << nl;
       }
   
   
   
       const label patchi = patch().index();
   
       tmp<scalarField> nuw = this->nu(patchi);
   
       const scalarListIOList & wallGradU =
           sampler_->db().lookupObject<scalarListIOList>("wallGradU");
   
       scalarField magGradU(Helpers::mag(wallGradU));
   
       return max
       (
           scalar(0),
           sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
       );
   }
   
   
   Foam::tmp<Foam::scalarField>
   Foam::ODEWallModelFvPatchScalarField::
   calcUTau(const scalarField & magGradU) const
   {
   
       const label patchi = patch().index();
       const label patchSize = patch().size();
   
       // Grab global uTau field
       volScalarField & uTauField =
           const_cast<volScalarField &>
           (
               db().lookupObject<volScalarField>("uTauPredicted")
           );
   
       const scalarField & uTauFieldBoundary = uTauField.boundaryFieldRef()[patchi];
   
       tmp<scalarField> tnuw = this->nu(patchi);
       const scalarField& nuw = tnuw();
   
       // vectorField for storing the source term
       vectorField sourceField(patchSize, vector(0, 0, 0));
   
       // Compute the source term
       source(sourceField);
   
       const scalarListIOList & U = sampler().db().lookupObject<scalarListIOList>("U");
       scalarField magU(patch().size());
   
       forAll(magU, i)
       {
           magU[i] = mag(vector(U[i][0], U[i][1], U[i][2]));
       }
   
       // Turbulent viscosity
       const scalarField & nutw = *this;
   
       // Computed friction velocity
       tmp<scalarField> tuTau(new scalarField(patchSize, 0.0));
   
       scalarField & uTau = tuTau.ref();
   
       AdaptiveIntegrator<scalar (scalar)> quad;
   
       // Compute uTau for each face
       forAll(uTau, faceI)
       {
           // Starting guess using resolved wall gradient or last timestep
           scalar tau = 0;
           if (uTauFieldBoundary[faceI] > 0)
           {
               tau = sqr(uTauFieldBoundary[faceI]);
           }
           else
           {
               tau = (nutw[faceI] + nuw[faceI]) * magGradU[faceI];
           }
   
           for (int iterI=0; iterI<maxIter_; iterI++)
           {
   
               IntegrandFunc eddyViscosity =
                   eddyViscosity_->value(sampler(), faceI, sqrt(tau), nuw[faceI]);
   
               scalar nuwI = nuw[faceI];
   
               IntegrandFunc integrand1 = [eddyViscosity, nuwI](const scalar y){
                   return 1.0/(nuwI + eddyViscosity(y));
               };
   
               IntegrandFunc integrand2 = [eddyViscosity, nuwI](const scalar y){
                   return y/(nuwI + eddyViscosity(y));
               };
   
               scalar integral1 =
                   quad.integrate(integrand1, 0.0, sampler().h()[faceI], 1e-6);
               scalar integral2 =
                   quad.integrate(integrand2, 0.0, sampler().h()[faceI], 1e-6);
   
               vector UFaceI(U[faceI][0], U[faceI][1], U[faceI][2]);
   
               scalar newTau =
                       sqr(magU[faceI]) + sqr(mag(sourceField[faceI])*integral2) -
                       2*(UFaceI & sourceField[faceI])*integral2;
   
               newTau  = sqrt(newTau)/(integral1 + VSMALL);
   
               scalar error = mag(tau - newTau)/(mag(tau) + ROOTVSMALL);
   
               tau = newTau;
               uTau[faceI] = sqrt(max(tau, scalar(0)));
   
               if (error < eps())
               {
                   break;
               }
   
               if (iterI == maxIter_-1)
               {
                   WarningIn
                   (
                       "Foam::ODEWallModelFvPatchScalarField::calcUTau()"
                   )
                       << "tau_w did not converge to desired tolerance "
                       << eps_ << ". Error value: " << error << nl;
               }
           }
   
       }
   
       // Assign computed uTau to the boundary field of the global field
       uTauField.boundaryFieldRef()[patch().index()] == uTau;
   
       return tuTau;
   }
   
   // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
   
   Foam::ODEWallModelFvPatchScalarField::
   ODEWallModelFvPatchScalarField
   (
       const fvPatch& p,
       const DimensionedField<scalar, volMesh>& iF
   )
   :
       wallModelFvPatchScalarField(p, iF),
       sampler_(nullptr),
       maxIter_(10),
       eps_(1e-3)
   {
   
       if (debug)
       {
           Info<< "Constructing ODEWallModelfvPatchScalarField (o1) "
               << "from copy and DimensionedField for patch " << patch().name()
               << nl;
       }
   }
   
   Foam::ODEWallModelFvPatchScalarField::
   ODEWallModelFvPatchScalarField
   (
       const ODEWallModelFvPatchScalarField& orig,
       const fvPatch& p,
       const DimensionedField<scalar, volMesh>& iF,
       const fvPatchFieldMapper& mapper
   )
   :
       wallModelFvPatchScalarField(orig, p, iF, mapper),
   #ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
       eddyViscosity_(orig.eddyViscosity_.clone()),
   #else
       eddyViscosity_(orig.eddyViscosity_, false),
   #endif
       sampler_(new SingleCellSampler(orig.sampler())),
       maxIter_(orig.maxIter_),
       eps_(orig.eps_)
   {
       if (debug)
       {
           Info<< "Constructing ODEWallModelfvPatchScalarField (o2) "
               << "from copy and DimensionedField for patch " << patch().name()
               << nl;
       }
   
       eddyViscosity_->addFieldsToSampler(sampler());
   }
   
   
   Foam::ODEWallModelFvPatchScalarField::
   ODEWallModelFvPatchScalarField
   (
       const fvPatch& p,
       const DimensionedField<scalar, volMesh>& iF,
       const dictionary& dict
   )
   :
       wallModelFvPatchScalarField(p, iF, dict),
       eddyViscosity_(EddyViscosity::New(dict.subDict("EddyViscosity"))),
       sampler_
       (
           new SingleCellSampler
           (
               p,
               averagingTime(),
               dict.lookupOrDefault<word>("interpolationType", "cell"),
               dict.lookupOrDefault<word>("sampler", "Tree"),
               dict.lookupOrDefault<word>
               (
                   "lengthScale",
                   dict.lookupOrDefault<word>("lengthScaleType", "CubeRootVol")
               ),
               dict.lookupOrDefault<bool>("hIsIndex", false)
           )
       ),
       maxIter_(dict.lookupOrDefault<label>("maxIter", 10)),
       eps_(dict.lookupOrDefault<scalar>("eps", 1e-3))
   
   {
       if (debug)
       {
           Info<< "Constructing ODEWallModelfvPatchScalarField (o3) "
               << "from copy and DimensionedField for patch " << patch().name()
               << nl;
       }
   
       eddyViscosity_->addFieldsToSampler(sampler());
   }
   
   
   #ifdef FOAM_FVPATCHFIELD_NO_COPY
   #else
   Foam::ODEWallModelFvPatchScalarField::
   ODEWallModelFvPatchScalarField
   (
       const ODEWallModelFvPatchScalarField& orig
   )
   :
       wallModelFvPatchScalarField(orig),
   #ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
       eddyViscosity_(orig.eddyViscosity_.clone()),
   #else
       eddyViscosity_(orig.eddyViscosity_, false),
   #endif
       sampler_(new SingleCellSampler(orig.sampler())),
       maxIter_(orig.maxIter_),
       eps_(orig.eps_)
   {
   
       if (debug)
       {
           Info<< "Constructing ODEWallModelfvPatchScalarField (o4) "
               << "from copy and DimensionedField for patch " << patch().name()
               << nl;
       }
       eddyViscosity_->addFieldsToSampler(sampler());
   }
   #endif
   
   
   Foam::ODEWallModelFvPatchScalarField::
   ODEWallModelFvPatchScalarField
   (
       const ODEWallModelFvPatchScalarField& orig,
       const DimensionedField<scalar, volMesh>& iF
   )
   :
       wallModelFvPatchScalarField(orig, iF),
   #ifdef FOAM_AUTOPTR_HAS_CLONE_METHOD
       eddyViscosity_(orig.eddyViscosity_.clone()),
   #else
       eddyViscosity_(orig.eddyViscosity_, false),
   #endif
       sampler_(new SingleCellSampler(orig.sampler_())),
       maxIter_(orig.maxIter_),
       eps_(orig.eps_)
   {
   
       if (debug)
       {
           Info<< "Constructing ODEWallModelfvPatchScalarField (o5) "
               << "from copy and DimensionedField for patch " << patch().name()
               << nl;
       }
       eddyViscosity_->addFieldsToSampler(sampler());
   }
   
   
   // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
   
   void Foam::ODEWallModelFvPatchScalarField::write(Ostream& os) const
   {
       wallModelFvPatchScalarField::write(os);
   }
   
   
   void Foam::ODEWallModelFvPatchScalarField::updateCoeffs()
   {
       if (updated())
       {
           return;
       }
   
   
       sampler().recomputeFields();
       sampler().sample();
   
       wallModelFvPatchScalarField::updateCoeffs();
   }
   
   // ************************************************************************* //