FreeFOAM The Cross-Platform CFD Toolkit
icoMomentError.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  icoMomentError
26 
27 Description
28  Estimates error for the incompressible laminar CFD application icoFoam.
29 
30 Usage
31 
32  - icoMomentError [OPTIONS]
33 
34  @param -noZero \n
35  Ignore timestep 0.
36 
37  @param -constant \n
38  Include the constant directory.
39 
40  @param -time <time>\n
41  Apply only to specific time.
42 
43  @param -latestTime \n
44  Only apply to latest time step.
45 
46  @param -case <dir>\n
47  Case directory.
48 
49  @param -parallel \n
50  Run in parallel.
51 
52  @param -help \n
53  Display help message.
54 
55  @param -doc \n
56  Display Doxygen API documentation page for this application.
57 
58  @param -srcDoc \n
59  Display Doxygen source documentation page for this application.
60 
61 \*---------------------------------------------------------------------------*/
62 
63 #include <finiteVolume/fvCFD.H>
64 #include <finiteVolume/linear.H>
67 
68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
69 
70 int main(int argc, char *argv[])
71 {
72  timeSelector::addOptions();
73 
74 # include <OpenFOAM/setRootCase.H>
75 # include <OpenFOAM/createTime.H>
76 
77  instantList timeDirs = timeSelector::select0(runTime, args);
78 
79 # include <OpenFOAM/createMesh.H>
80 
81  Info<< "\nEstimating error in the incompressible momentum equation\n"
82  << "Reading transportProperties\n" << endl;
83 
85  (
86  IOobject
87  (
88  "transportProperties",
89  runTime.constant(),
90  mesh,
91  IOobject::MUST_READ,
92  IOobject::NO_WRITE
93  )
94  );
95 
97  (
99  );
100 
101  forAll(timeDirs, timeI)
102  {
103  runTime.setTime(timeDirs[timeI], timeI);
104 
105  Info<< "Time = " << runTime.timeName() << endl;
106 
107  mesh.readUpdate();
108 
109  IOobject pHeader
110  (
111  "p",
112  runTime.timeName(),
113  mesh,
114  IOobject::MUST_READ
115  );
116 
117  IOobject Uheader
118  (
119  "U",
120  runTime.timeName(),
121  mesh,
122  IOobject::MUST_READ
123  );
124 
125  if (pHeader.headerOk() && Uheader.headerOk())
126  {
127  Info << "Reading p" << endl;
128  volScalarField p(pHeader, mesh);
129 
130  Info << "Reading U" << endl;
131  volVectorField U(Uheader, mesh);
132 
133 # include <finiteVolume/createPhi.H>
134 
135  volScalarField ek = 0.5*magSqr(U);
136  volTensorField gradU = fvc::grad(U);
137 
138  // Divergence of the error in U squared
139 
141  (
142  IOobject
143  (
144  "L",
145  mesh.time().timeName(),
146  mesh,
147  IOobject::NO_READ,
148  IOobject::NO_WRITE
149  ),
150  mesh,
151  dimensionedScalar("one", dimLength, 1.0)
152  );
153 
154  L.internalField() =
155  mesh.V()/fvc::surfaceSum(mesh.magSf())().internalField();
156 
157  // Warning: 4th row of this equation specially modified
158  // for the momentum equation. The "real" formulation would
159  // have diffusivity*(gradV && gradV)
160  volScalarField momError
161  (
162  IOobject
163  (
164  "momErrorL" + U.name(),
165  mesh.time().timeName(),
166  mesh,
167  IOobject::NO_READ,
168  IOobject::NO_WRITE
169  ),
170  sqrt
171  (
172  2.0*mag
173  (
174  (
176  (
177  mesh,
178  phi,
180  (
181  new linear<scalar>(mesh)
182  )
183  ).fvcDiv(phi, ek)
184 
185  - nu*
187  .fvcLaplacian
188  (
189  ek
190  )
191  - (U & fvc::grad(p))
192 // + nu*(gradU && gradU)
193  + 0.5*nu*
194  (
195  gradU && (gradU + gradU.T())
196  )
197  )*L/(mag(U) + nu/L)
198  )
199  )
200  );
201 
202  momError.boundaryField() = 0.0;
203  momError.write();
204  }
205  else
206  {
207  Info<< " No p or U" << endl;
208  }
209 
210  Info<< endl;
211  }
212 
213  Info<< "End\n" << endl;
214 
215  return 0;
216 }
217 
218 
219 // ************************ vim: set sw=4 sts=4 et: ************************ //