src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedCubic/limitedCubic.H

   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 Class
  25     Foam::limitedCubicLimiter
  26
  27 Description
  28     Class with limiter function which returns the limiter for the
  29     TVD limited centred-cubic differencing scheme based on r obtained from
  30     the LimiterFunc class.
  31
  32     Used in conjunction with the template class LimitedScheme.
  33
  34 SourceFiles
  35     limitedCubic.C
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #ifndef limitedCubic_H
  40 #define limitedCubic_H
  41
  42 #include "vector.H"
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 namespace Foam
  47 {
  48
  49 /*---------------------------------------------------------------------------*\
  50                       Class limitedCubicLimiter Declaration
  51 \*---------------------------------------------------------------------------*/
  52
  53 template<class LimiterFunc>
  54 class limitedCubicLimiter
  55 :
  56     public LimiterFunc
  57 {
  58     scalar k_;
  59     scalar twoByk_;
  60
  61 public:
  62
  63     limitedCubicLimiter(Istream& is)
  64     :
  65         k_(readScalar(is))
  66     {
  67         if (k_ < 0 || k_ > 1)
  68         {
  69             FatalIOErrorIn("limitedCubicLimiter(Istream& is)", is)
  70                 << "coefficient = " << k_
  71                 << " should be >= 0 and <= 1"
  72                 << exit(FatalIOError);
  73         }
  74
  75         // Avoid the /0 when k_ = 0
  76         twoByk_ = 2.0/max(k_, SMALL);
  77     }
  78
  79     scalar limiter
  80     (
  81         const scalar cdWeight,
  82         const scalar faceFlux,
  83         const typename LimiterFunc::phiType& phiP,
  84         const typename LimiterFunc::phiType& phiN,
  85         const typename LimiterFunc::gradPhiType& gradcP,
  86         const typename LimiterFunc::gradPhiType& gradcN,
  87         const vector& d
  88     ) const
  89     {
  90         scalar twor = twoByk_*LimiterFunc::r
  91         (
  92             faceFlux, phiP, phiN, gradcP, gradcN, d
  93         );
  94
  95         scalar phiU;
  96
  97         if (faceFlux > 0)
  98         {
  99             phiU = phiP;
 100         }
 101         else
 102         {
 103             phiU = phiN;
 104         }
 105
 106         // Calculate the face value using cubic interpolation
 107         scalar phif =
 108             cdWeight*(phiP - 0.25*(d & gradcN))
 109           + (1 - cdWeight)*(phiN + 0.25*(d & gradcP));
 110
 111         scalar phiCD = cdWeight*phiP + (1 - cdWeight)*phiN;
 112
 113         // Calculate the effective limiter for the cubic interpolation
 114         scalar cubicLimiter = (phif - phiU)/stabilise(phiCD - phiU, SMALL);
 115
 116         /*
 117         if (twor < 0.05)
 118         {
 119             cubicLimiter = twor;
 120         }
 121
 122         return max(min(cubicLimiter, 2), 0);
 123         */
 124
 125         // Limit the limiter to obey the TVD constraint
 126         return max(min(min(twor, cubicLimiter), 2), 0);
 127     }
 128 };
 129
 130
 131 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 132
 133 } // End namespace Foam
 134
 135 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 136
 137 #endif
 138
 139 // ************************************************************************* //