src/finiteVolume/finiteVolume/gradSchemes/fourthGrad/fourthGrad.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 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 "fourthGrad.H"
  28 #include "leastSquaresGrad.H"
  29 #include "gaussGrad.H"
  30 #include "fvMesh.H"
  31 #include "volMesh.H"
  32 #include "surfaceMesh.H"
  33 #include "GeometricField.H"
  34 #include "zeroGradientFvPatchField.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace fv
  44 {
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 template<class Type>
  49 tmp
  50 <
  51     GeometricField
  52     <
  53         typename outerProduct<vector, Type>::type, fvPatchField, volMesh
  54     >
  55 >
  56 fourthGrad<Type>::grad
  57 (
  58     const GeometricField<Type, fvPatchField, volMesh>& vsf
  59 ) const
  60 {
  61     // The fourth-order gradient is calculated in two passes.  First,
  62     // the standard least-square gradient is assembled.  Then, the
  63     // fourth-order correction is added to the second-order accurate
  64     // gradient to complete the accuracy.
  65
  66     typedef typename outerProduct<vector, Type>::type GradType;
  67
  68     const fvMesh& mesh = vsf.mesh();
  69
  70     // Assemble the second-order least-square gradient
  71     // Calculate the second-order least-square gradient
  72     tmp<GeometricField<GradType, fvPatchField, volMesh> > tsecondfGrad
  73         = leastSquaresGrad<Type>(mesh).grad(vsf);
  74     const GeometricField<GradType, fvPatchField, volMesh>& secondfGrad =
  75         tsecondfGrad();
  76
  77     tmp<GeometricField<GradType, fvPatchField, volMesh> > tfGrad
  78     (
  79         new GeometricField<GradType, fvPatchField, volMesh>
  80         (
  81             IOobject
  82             (
  83                 "grad("+vsf.name()+')',
  84                 vsf.instance(),
  85                 mesh,
  86                 IOobject::NO_READ,
  87                 IOobject::NO_WRITE
  88             ),
  89             secondfGrad
  90         )
  91     );
  92     GeometricField<GradType, fvPatchField, volMesh>& fGrad = tfGrad();
  93
  94     const vectorField& C = mesh.C();
  95
  96     const surfaceScalarField& lambda = mesh.weights();
  97
  98     // Get reference to least square vectors
  99     const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
 100     const surfaceVectorField& ownLs = lsv.pVectors();
 101     const surfaceVectorField& neiLs = lsv.nVectors();
 102
 103     // owner/neighbour addressing
 104     const unallocLabelList& own = mesh.owner();
 105     const unallocLabelList& nei = mesh.neighbour();
 106
 107     // Assemble the fourth-order gradient
 108
 109     // Internal faces
 110     forAll(own, facei)
 111     {
 112         Type dDotGradDelta = 0.5*
 113         (
 114            (C[nei[facei]] - C[own[facei]])
 115          & (secondfGrad[nei[facei]] - secondfGrad[own[facei]])
 116         );
 117
 118         fGrad[own[facei]] -= lambda[facei]*ownLs[facei]*dDotGradDelta;
 119         fGrad[nei[facei]] -= (1.0 - lambda[facei])*neiLs[facei]*dDotGradDelta;
 120     }
 121
 122     // Boundary faces
 123     forAll(vsf.boundaryField(), patchi)
 124     {
 125         if (secondfGrad.boundaryField()[patchi].coupled())
 126         {
 127             const fvsPatchVectorField& patchOwnLs =
 128                 ownLs.boundaryField()[patchi];
 129
 130             const scalarField& lambdap = lambda.boundaryField()[patchi];
 131
 132             // Build the d-vectors
 133             vectorField pd =
 134                 mesh.Sf().boundaryField()[patchi]
 135                /(
 136                    mesh.magSf().boundaryField()[patchi]
 137                   *mesh.deltaCoeffs().boundaryField()[patchi]
 138                 );
 139
 140             if (!mesh.orthogonal())
 141             {
 142                 pd -= mesh.correctionVectors().boundaryField()[patchi]
 143                      /mesh.deltaCoeffs().boundaryField()[patchi];
 144             }
 145
 146             const unallocLabelList& faceCells =
 147                 fGrad.boundaryField()[patchi].patch().faceCells();
 148
 149             Field<GradType> neighbourSecondfGrad =
 150                 secondfGrad.boundaryField()[patchi].patchNeighbourField();
 151
 152             forAll(faceCells, patchFaceI)
 153             {
 154                 fGrad[faceCells[patchFaceI]] -=
 155                     0.5*lambdap[patchFaceI]*patchOwnLs[patchFaceI]
 156                    *(
 157                         pd[patchFaceI]
 158                       & (
 159                             neighbourSecondfGrad[patchFaceI]
 160                           - secondfGrad[faceCells[patchFaceI]]
 161                         )
 162                     );
 163             }
 164         }
 165     }
 166
 167     fGrad.correctBoundaryConditions();
 168     gaussGrad<Type>::correctBoundaryConditions(vsf, fGrad);
 169
 170     return tfGrad;
 171 }
 172
 173
 174 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 175
 176 } // End namespace fv
 177
 178 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 179
 180 } // End namespace Foam
 181
 182 // ************************************************************************* //