1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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
29 Estimates error for the incompressible laminar CFD application icoFoam.
31 \*---------------------------------------------------------------------------*/
35 #include "gaussConvectionScheme.H"
36 #include "gaussLaplacianScheme.H"
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 int main(int argc, char *argv[])
42 timeSelector::addOptions();
44 # include "setRootCase.H"
45 # include "createTime.H"
47 instantList timeDirs = timeSelector::select0(runTime, args);
49 # include "createMesh.H"
51 Info<< "\nEstimating error in the incompressible momentum equation\n"
52 << "Reading transportProperties\n" << endl;
54 IOdictionary transportProperties
58 "transportProperties",
68 transportProperties.lookup("nu")
71 forAll(timeDirs, timeI)
73 runTime.setTime(timeDirs[timeI], timeI);
75 Info<< "Time = " << runTime.timeName() << endl;
95 if (pHeader.headerOk() && Uheader.headerOk())
97 Info << "Reading p" << endl;
98 volScalarField p(pHeader, mesh);
100 Info << "Reading U" << endl;
101 volVectorField U(Uheader, mesh);
103 # include "createPhi.H"
105 volScalarField ek = 0.5*magSqr(U);
106 volTensorField gradU = fvc::grad(U);
108 // Divergence of the error in U squared
115 mesh.time().timeName(),
121 dimensionedScalar("one", dimLength, 1.0)
125 mesh.V()/fvc::surfaceSum(mesh.magSf())().internalField();
127 // Warning: 4th row of this equation specially modified
128 // for the momentum equation. The "real" formulation would
129 // have diffusivity*(gradV && gradV)
130 volScalarField momError
134 "momErrorL" + U.name(),
135 mesh.time().timeName(),
145 fv::gaussConvectionScheme<scalar>
149 tmp<surfaceInterpolationScheme<scalar> >
151 new linear<scalar>(mesh)
156 fv::gaussLaplacianScheme<scalar, scalar>(mesh)
162 // + nu*(gradU && gradU)
165 gradU && (gradU + gradU.T())
172 momError.boundaryField() = 0.0;
177 Info<< " No p or U" << endl;
183 Info<< "End\n" << endl;
189 // ************************************************************************* //