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 the error in the solution for a scalar transport equation in the
32 \*---------------------------------------------------------------------------*/
36 #include "gaussConvectionScheme.H"
37 #include "gaussLaplacianScheme.H"
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 int main(int argc, char *argv[])
43 timeSelector::addOptions();
45 # include "setRootCase.H"
46 # include "createTime.H"
48 instantList timeDirs = timeSelector::select0(runTime, args);
50 # include "createMesh.H"
52 Info<< "\nEstimating error in scalar transport equation\n"
53 << "Reading transportProperties\n" << endl;
55 IOdictionary transportProperties
59 "transportProperties",
68 Info<< "Reading diffusivity DT\n" << endl;
72 transportProperties.lookup("DT")
76 forAll(timeDirs, timeI)
78 runTime.setTime(timeDirs[timeI], timeI);
80 Info<< "Time = " << runTime.timeName() << endl;
100 if (THeader.headerOk() && Uheader.headerOk())
102 Info << "Reading T" << endl;
103 volScalarField T(THeader, mesh);
105 Info << "Reading U" << endl;
106 volVectorField U(Uheader, mesh);
108 # include "createPhi.H"
110 volVectorField gradT = fvc::grad(T);
112 volScalarField TE = 0.5*sqr(T);
119 mesh.time().timeName(),
125 dimensionedScalar("one", dimLength, 1.0)
129 mesh.V()/fvc::surfaceSum(mesh.magSf())().internalField();
131 // Divergence of the error in the T squared
132 volScalarField momError
136 "momErrorL" + T.name(),
137 mesh.time().timeName(),
147 fv::gaussConvectionScheme<scalar>
151 tmp<surfaceInterpolationScheme<scalar> >
153 new linear<scalar>(mesh)
158 fv::gaussLaplacianScheme<scalar, scalar>(mesh)
169 momError.boundaryField() = 0.0;
174 Info<< " No T or U" << endl;
180 Info<< "End\n" << endl;
186 // ************************************************************************* //