src/vofDynamicMesh/CICSAM/CICSAM.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "CICSAM.H"
  28 #include "fvc.H"
  29 #include "volFields.H"
  30 #include "surfaceFields.H"
  31 #include "upwind.H"
  32
  33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  34
  35 Foam::scalar Foam::CICSAM::limiter
  36 (
  37     const scalar cdWeight,
  38     const scalar faceFlux,
  39     const scalar& phiP,
  40     const scalar& phiN,
  41     const vector& gradcP,
  42     const vector& gradcN,
  43     const scalar Cof,
  44     const vector d
  45 ) const
  46 {
  47     // Calculate CICSAM weight
  48     scalar w =
  49         weight
  50         (
  51             cdWeight,
  52             faceFlux,
  53             phiP,
  54             phiN,
  55             gradcP,
  56             gradcN,
  57             Cof,
  58             d
  59         );
  60
  61     // Calculate CICSAM face value
  62     scalar phif = w*phiP + (1 - w)*phiN;
  63
  64     // Calculate UD and CD face value
  65     scalar phiU = faceFlux >= 0 ? phiP : phiN;
  66     scalar phiCD = cdWeight*phiP + (1 - cdWeight)*phiN;
  67
  68     // Calculate the effective limiter for the QUICK interpolation
  69     scalar CLimiter = (phif - phiU)/stabilise(phiCD - phiU, SMALL);
  70
  71     // Limit the limiter between upwind and downwind
  72     return max(min(CLimiter, 2), 0);
  73 }
  74
  75
  76 Foam::scalar Foam::CICSAM::weight
  77 (
  78     const scalar cdWeight,
  79     const scalar faceFlux,
  80     const scalar& phiP,
  81     const scalar& phiN,
  82     const vector& gradcP,
  83     const vector& gradcN,
  84     const scalar Cof,
  85     const vector d
  86 ) const
  87 {
  88     // Calculate upwind value, faceFlux C tilde and do a stabilisation
  89
  90     scalar phict = 0;
  91     scalar phiupw = 0;
  92     scalar costheta = 0;
  93
  94     if (faceFlux > 0)
  95     {
  96         costheta = mag((gradcP & d)/(mag(gradcP)*mag(d) + SMALL));
  97
  98         phiupw = phiN - 2*(gradcP & d);
  99
 100         phiupw = max(min(phiupw, 1.0), 0.0);
 101
 102         if ((phiN - phiupw) > 0)
 103         {
 104             phict = (phiP - phiupw)/(phiN - phiupw + SMALL);
 105         }
 106         else
 107         {
 108             phict = (phiP - phiupw)/(phiN - phiupw - SMALL);
 109         }
 110     }
 111     else
 112     {
 113         costheta = mag((gradcN & d)/(mag(gradcN)*mag(d) + SMALL));
 114
 115         phiupw = phiP + 2*(gradcN & d);
 116
 117         phiupw = max(min(phiupw, 1.0), 0.0);
 118
 119         if ((phiP - phiupw) > 0)
 120         {
 121             phict = (phiN - phiupw)/(phiP - phiupw + SMALL);
 122         }
 123         else
 124         {
 125             phict = (phiN - phiupw)/(phiP - phiupw - SMALL);
 126         }
 127     }
 128
 129
 130     // Calculate the weighting factors for CICSAM
 131
 132     scalar cicsamFactor = (k_ + SMALL)/(1 - k_ + SMALL);
 133
 134     costheta = min(1.0, cicsamFactor*(costheta));
 135     costheta = (cos(2*(acos(costheta))) + 1)/2;
 136
 137     scalar k1 = (3*Cof*Cof - 3*Cof)/(2*Cof*Cof + 6*Cof - 8);
 138     scalar k2 = Cof;
 139     scalar k3 = (3*Cof + 5)/(2*Cof + 6);
 140     scalar weight;
 141
 142     if (phict > 0 && phict <= k1)             // use blended scheme 1
 143     {
 144         scalar phifCM = phict/(Cof + SMALL);
 145         weight = (phifCM - phict)/(1 - phict);
 146     }
 147     else if (phict > k1 && phict <= k2)     // use blended scheme 2
 148     {
 149         scalar phifHC = phict/(Cof + SMALL);
 150         scalar phifUQ = (8*Cof*phict + (1 - Cof)*(6*phict + 3))/8;
 151         scalar phifCM = costheta*phifHC + (1 - costheta)*phifUQ;
 152         weight = (phifCM - phict)/(1 - phict);
 153     }
 154     else if (phict > k2 && phict < k3)     // use blended scheme 3
 155     {
 156         scalar phifUQ = (8*Cof*phict + (1 - Cof)*(6*phict + 3))/8;
 157         scalar phifCM = costheta + (1 - costheta)*phifUQ;
 158         weight = (phifCM - phict)/(1 - phict);
 159     }
 160     else if (phict >= k3 && phict <= 1)     // use downwind
 161     {
 162         weight = 1;
 163     }
 164     else                                    // use upwind
 165     {
 166         weight = 0;
 167     }
 168
 169     if (faceFlux > 0)
 170     {
 171         return 1 - weight;
 172     }
 173     else
 174     {
 175         return weight;
 176     }
 177 }
 178
 179
 180 Foam::tmp<Foam::surfaceScalarField> Foam::CICSAM::limiter
 181 (
 182     const volScalarField& phi
 183 ) const
 184 {
 185     const fvMesh& mesh = this->mesh();
 186
 187     tmp<surfaceScalarField> tLimiter
 188     (
 189         new surfaceScalarField
 190         (
 191             IOobject
 192             (
 193                 type() + "Limiter(" + phi.name() + ')',
 194                 mesh.time().timeName(),
 195                 mesh
 196             ),
 197             mesh,
 198             dimless
 199         )
 200     );
 201     surfaceScalarField& lim = tLimiter();
 202
 203     volVectorField gradc = fvc::grad(phi);
 204
 205     surfaceScalarField Cof =
 206         mesh.time().deltaT()
 207        *upwind<scalar>(mesh, faceFlux_).interpolate
 208         (
 209             fvc::surfaceIntegrate(faceFlux_)
 210         );
 211
 212     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
 213
 214     const unallocLabelList& owner = mesh.owner();
 215     const unallocLabelList& neighbour = mesh.neighbour();
 216
 217     const vectorField& C = mesh.C();
 218
 219     scalarField& pLim = lim.internalField();
 220
 221     forAll(pLim, faceI)
 222     {
 223         label own = owner[faceI];
 224         label nei = neighbour[faceI];
 225
 226         pLim[faceI] = limiter
 227         (
 228             CDweights[faceI],
 229             this->faceFlux_[faceI],
 230             phi[own],
 231             phi[nei],
 232             gradc[own],
 233             gradc[nei],
 234             Cof[faceI],
 235             C[nei] - C[own]
 236         );
 237     }
 238
 239     surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
 240     surfaceScalarField::GeometricBoundaryField& bCof = Cof.boundaryField();
 241     forAll(bLim, patchi)
 242     {
 243         scalarField& pLim = bLim[patchi];
 244
 245         if (bLim[patchi].coupled())
 246         {
 247             const scalarField& pCDweights = CDweights.boundaryField()[patchi];
 248
 249             const scalarField& pFaceFlux =
 250                 this->faceFlux_.boundaryField()[patchi];
 251
 252             scalarField pphiP =
 253                 phi.boundaryField()[patchi].patchInternalField();
 254
 255             scalarField pphiN =
 256                 phi.boundaryField()[patchi].patchNeighbourField();
 257
 258             vectorField pGradcP =
 259                 gradc.boundaryField()[patchi].patchInternalField();
 260
 261             vectorField pGradcN =
 262                 gradc.boundaryField()[patchi].patchNeighbourField();
 263
 264             const scalarField& pCof = Cof.boundaryField()[patchi];
 265
 266             // Build the d-vectors
 267             // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
 268             vectorField pd = bLim[patchi].patch().delta();
 269
 270             forAll(pLim, faceI)
 271             {
 272                 pLim[faceI] = limiter
 273                 (
 274                     pCDweights[faceI],
 275                     pFaceFlux[faceI],
 276                     pphiP[faceI],
 277                     pphiN[faceI],
 278                     pGradcP[faceI],
 279                     pGradcN[faceI],
 280                     pCof[faceI],
 281                     pd[faceI]
 282                 );
 283             }
 284         }
 285         else
 286         {
 287             pLim = 1.0;
 288         }
 289     }
 290
 291     return tLimiter;
 292 }
 293
 294
 295
 296 Foam::tmp<Foam::surfaceScalarField> Foam::CICSAM::weights
 297 (
 298     const volScalarField& phi
 299 ) const
 300 {
 301     const fvMesh& mesh = this->mesh();
 302
 303     tmp<surfaceScalarField> tWeightingFactors
 304     (
 305         new surfaceScalarField(mesh.surfaceInterpolation::weights())
 306     );
 307     surfaceScalarField& weightingFactors = tWeightingFactors();
 308
 309     scalarField& weights = weightingFactors.internalField();
 310
 311
 312     volVectorField gradc = fvc::grad(phi);
 313
 314     surfaceScalarField Cof =
 315         mesh.time().deltaT()
 316        *upwind<scalar>(mesh, faceFlux_).interpolate
 317         (
 318             fvc::surfaceIntegrate(faceFlux_)
 319         );
 320
 321     const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
 322
 323     const unallocLabelList& owner = mesh.owner();
 324     const unallocLabelList& neighbour = mesh.neighbour();
 325
 326     const vectorField& C = mesh.C();
 327
 328     scalarField& w = weightingFactors.internalField();
 329
 330     forAll(w, faceI)
 331     {
 332         label own = owner[faceI];
 333         label nei = neighbour[faceI];
 334
 335         w[faceI] = weight
 336         (
 337             CDweights[faceI],
 338             this->faceFlux_[faceI],
 339             phi[own],
 340             phi[nei],
 341             gradc[own],
 342             gradc[nei],
 343             Cof[faceI],
 344             C[nei] - C[own]
 345         );
 346     }
 347
 348     surfaceScalarField::GeometricBoundaryField& bWeights =
 349         weightingFactors.boundaryField();
 350
 351     surfaceScalarField::GeometricBoundaryField& bCof = Cof.boundaryField();
 352
 353     forAll(bWeights, patchi)
 354     {
 355         scalarField& pWeights = bWeights[patchi];
 356
 357         if (bWeights[patchi].coupled())
 358         {
 359             const scalarField& pCDweights = CDweights.boundaryField()[patchi];
 360
 361             const scalarField& pFaceFlux =
 362                 this->faceFlux_.boundaryField()[patchi];
 363
 364             scalarField pphiP =
 365                 phi.boundaryField()[patchi].patchInternalField();
 366
 367             scalarField pphiN =
 368                 phi.boundaryField()[patchi].patchNeighbourField();
 369
 370             vectorField pGradcP =
 371                 gradc.boundaryField()[patchi].patchInternalField();
 372
 373             vectorField pGradcN =
 374                 gradc.boundaryField()[patchi].patchNeighbourField();
 375
 376             const scalarField& pCof = Cof.boundaryField()[patchi];
 377
 378             // Build the d-vectors
 379             // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
 380             vectorField pd = bWeights[patchi].patch().delta();
 381
 382             forAll(pWeights, faceI)
 383             {
 384                 pWeights[faceI] = weight
 385                 (
 386                     pCDweights[faceI],
 387                     pFaceFlux[faceI],
 388                     pphiP[faceI],
 389                     pphiN[faceI],
 390                     pGradcP[faceI],
 391                     pGradcN[faceI],
 392                     pCof[faceI],
 393                     pd[faceI]
 394                 );
 395             }
 396         }
 397         else
 398         {
 399             pWeights = 1.0;
 400         }
 401     }
 402
 403     return tWeightingFactors;
 404 }
 405
 406
 407 namespace Foam
 408 {
 409 defineNamedTemplateTypeNameAndDebug(CICSAM, 0);
 410
 411 surfaceInterpolationScheme<scalar>::addMeshConstructorToTable<CICSAM>
 412     addCICSAMMeshConstructorToTable_;
 413
 414 surfaceInterpolationScheme<scalar>::addMeshFluxConstructorToTable<CICSAM>
 415     addCICSAMMeshFluxConstructorToTable_;
 416
 417 limitedSurfaceInterpolationScheme<scalar>::addMeshConstructorToTable<CICSAM>
 418     addCICSAMMeshConstructorToLimitedTable_;
 419
 420 limitedSurfaceInterpolationScheme<scalar>::
 421 addMeshFluxConstructorToTable<CICSAM>
 422     addCICSAMMeshFluxConstructorToLimitedTable_;
 423 }
 424
 425 // ************************************************************************* //