FreeFOAM The Cross-Platform CFD Toolkit
momentScalarError.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  momentScalarError
26 
27 Description
28  Estimates the error in the solution for a scalar transport equation in the
29  standard form
30 
31 Usage
32 
33  - momentScalarError [OPTIONS]
34 
35  @param -noZero \n
36  Ignore timestep 0.
37 
38  @param -constant \n
39  Include the constant directory.
40 
41  @param -time <time>\n
42  Apply only to specific time.
43 
44  @param -latestTime \n
45  Only apply to latest time step.
46 
47  @param -case <dir>\n
48  Case directory.
49 
50  @param -parallel \n
51  Run in parallel.
52 
53  @param -help \n
54  Display help message.
55 
56  @param -doc \n
57  Display Doxygen API documentation page for this application.
58 
59  @param -srcDoc \n
60  Display Doxygen source documentation page for this application.
61 
62 \*---------------------------------------------------------------------------*/
63 
64 #include <finiteVolume/fvCFD.H>
65 #include <finiteVolume/linear.H>
68 
69 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70 
71 int main(int argc, char *argv[])
72 {
73  timeSelector::addOptions();
74 
75 # include <OpenFOAM/setRootCase.H>
76 # include <OpenFOAM/createTime.H>
77 
78  instantList timeDirs = timeSelector::select0(runTime, args);
79 
80 # include <OpenFOAM/createMesh.H>
81 
82  Info<< "\nEstimating error in scalar transport equation\n"
83  << "Reading transportProperties\n" << endl;
84 
86  (
87  IOobject
88  (
89  "transportProperties",
90  runTime.constant(),
91  mesh,
92  IOobject::MUST_READ,
93  IOobject::NO_WRITE
94  )
95  );
96 
97 
98  Info<< "Reading diffusivity DT\n" << endl;
99 
101  (
103  );
104 
105 
106  forAll(timeDirs, timeI)
107  {
108  runTime.setTime(timeDirs[timeI], timeI);
109 
110  Info<< "Time = " << runTime.timeName() << endl;
111 
112  mesh.readUpdate();
113 
114  IOobject THeader
115  (
116  "T",
117  runTime.timeName(),
118  mesh,
119  IOobject::MUST_READ
120  );
121 
122  IOobject Uheader
123  (
124  "U",
125  runTime.timeName(),
126  mesh,
127  IOobject::MUST_READ
128  );
129 
130  if (THeader.headerOk() && Uheader.headerOk())
131  {
132  Info << "Reading T" << endl;
133  volScalarField T(THeader, mesh);
134 
135  Info << "Reading U" << endl;
136  volVectorField U(Uheader, mesh);
137 
138 # include <finiteVolume/createPhi.H>
139 
140  volVectorField gradT = fvc::grad(T);
141 
142  volScalarField TE = 0.5*sqr(T);
143 
145  (
146  IOobject
147  (
148  "L",
149  mesh.time().timeName(),
150  mesh,
151  IOobject::NO_READ,
152  IOobject::NO_WRITE
153  ),
154  mesh,
155  dimensionedScalar("one", dimLength, 1.0)
156  );
157 
158  L.internalField() =
159  mesh.V()/fvc::surfaceSum(mesh.magSf())().internalField();
160 
161  // Divergence of the error in the T squared
162  volScalarField momError
163  (
164  IOobject
165  (
166  "momErrorL" + T.name(),
167  mesh.time().timeName(),
168  mesh,
169  IOobject::NO_READ,
170  IOobject::NO_WRITE
171  ),
172  sqrt
173  (
174  2.0*mag
175  (
176  (
178  (
179  mesh,
180  phi,
182  (
183  new linear<scalar>(mesh)
184  )
185  ).fvcDiv(phi, TE)
186 
187  - DT*
189  .fvcLaplacian
190  (
191  TE
192  )
193  + DT*(gradT & gradT)
194  )*L/(mag(U) + DT/L)
195  )
196  )
197  );
198 
199  momError.boundaryField() = 0.0;
200  momError.write();
201  }
202  else
203  {
204  Info<< " No T or U" << endl;
205  }
206 
207  Info<< endl;
208  }
209 
210  Info<< "End\n" << endl;
211 
212  return 0;
213 }
214 
215 
216 // ************************ vim: set sw=4 sts=4 et: ************************ //