applications/utilities/errorEstimation/icoMomentError/icoMomentError.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     icoMomentError
  26
  27 Description
  28     Estimates error for the incompressible laminar CFD application icoFoam.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "fvCFD.H"
  33 #include "linear.H"
  34 #include "gaussConvectionScheme.H"
  35 #include "gaussLaplacianScheme.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 int main(int argc, char *argv[])
  40 {
  41     timeSelector::addOptions();
  42
  43 #   include "setRootCase.H"
  44 #   include "createTime.H"
  45
  46     instantList timeDirs = timeSelector::select0(runTime, args);
  47
  48 #   include "createMesh.H"
  49
  50     Info<< "\nEstimating error in the incompressible momentum equation\n"
  51         << "Reading transportProperties\n" << endl;
  52
  53     IOdictionary transportProperties
  54     (
  55         IOobject
  56         (
  57             "transportProperties",
  58             runTime.constant(),
  59             mesh,
  60             IOobject::MUST_READ,
  61             IOobject::NO_WRITE
  62         )
  63     );
  64
  65     dimensionedScalar nu
  66     (
  67         transportProperties.lookup("nu")
  68     );
  69
  70     forAll(timeDirs, timeI)
  71     {
  72         runTime.setTime(timeDirs[timeI], timeI);
  73
  74         Info<< "Time = " << runTime.timeName() << endl;
  75
  76         mesh.readUpdate();
  77
  78         IOobject pHeader
  79         (
  80             "p",
  81             runTime.timeName(),
  82             mesh,
  83             IOobject::MUST_READ
  84         );
  85
  86         IOobject Uheader
  87         (
  88             "U",
  89             runTime.timeName(),
  90             mesh,
  91             IOobject::MUST_READ
  92         );
  93
  94         if (pHeader.headerOk() && Uheader.headerOk())
  95         {
  96             Info << "Reading p" << endl;
  97             volScalarField p(pHeader, mesh);
  98
  99             Info << "Reading U" << endl;
 100             volVectorField U(Uheader, mesh);
 101
 102 #           include "createPhi.H"
 103
 104             volScalarField ek = 0.5*magSqr(U);
 105             volTensorField gradU = fvc::grad(U);
 106
 107             // Divergence of the error in U squared
 108
 109             volScalarField L
 110             (
 111                 IOobject
 112                 (
 113                     "L",
 114                     mesh.time().timeName(),
 115                     mesh,
 116                     IOobject::NO_READ,
 117                     IOobject::NO_WRITE
 118                 ),
 119                 mesh,
 120                 dimensionedScalar("one", dimLength, 1.0)
 121             );
 122
 123             L.internalField() =
 124                 mesh.V()/fvc::surfaceSum(mesh.magSf())().internalField();
 125
 126             // Warning: 4th row of this equation specially modified
 127             // for the momentum equation. The "real" formulation would
 128             // have diffusivity*(gradV && gradV)
 129             volScalarField momError
 130             (
 131                 IOobject
 132                 (
 133                     "momErrorL" + U.name(),
 134                     mesh.time().timeName(),
 135                     mesh,
 136                     IOobject::NO_READ,
 137                     IOobject::NO_WRITE
 138                 ),
 139                 sqrt
 140                 (
 141                     2.0*mag
 142                     (
 143                         (
 144                             fv::gaussConvectionScheme<scalar>
 145                             (
 146                                 mesh,
 147                                 phi,
 148                                 tmp<surfaceInterpolationScheme<scalar> >
 149                                 (
 150                                     new linear<scalar>(mesh)
 151                                 )
 152                             ).fvcDiv(phi, ek)
 153
 154                           - nu*
 155                             fv::gaussLaplacianScheme<scalar, scalar>(mesh)
 156                            .fvcLaplacian
 157                             (
 158                                 ek
 159                             )
 160                           - (U & fvc::grad(p))
 161 //                        + nu*(gradU && gradU)
 162                           + 0.5*nu*
 163                             (
 164                                 gradU && (gradU + gradU.T())
 165                             )
 166                         )*L/(mag(U) + nu/L)
 167                     )
 168                 )
 169             );
 170
 171             momError.boundaryField() = 0.0;
 172             momError.write();
 173         }
 174         else
 175         {
 176             Info<< "    No p or U" << endl;
 177         }
 178
 179         Info<< endl;
 180     }
 181
 182     Info<< "End\n" << endl;
 183
 184     return 0;
 185 }
 186
 187
 188 // ************************************************************************* //