FreeFOAM The Cross-Platform CFD Toolkit
shallowWaterFoam.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  shallowWaterFoam
26 
27 Description
28  Transient solver for inviscid shallow-water equations with rotation.
29 
30  If the geometry is 3D then it is assumed to be one layers of cells and the
31  component of the velocity normal to gravity is removed.
32 
33 Usage
34  - shallowWaterFoam [OPTION]
35 
36  @param -case <dir> \n
37  Specify the case directory
38 
39  @param -parallel \n
40  Run the case in parallel
41 
42  @param -help \n
43  Display short usage message
44 
45  @param -doc \n
46  Display Doxygen documentation page
47 
48  @param -srcDoc \n
49  Display source code
50 
51 \*---------------------------------------------------------------------------*/
52 
53 #include <finiteVolume/fvCFD.H>
54 
55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
56 
57 int main(int argc, char *argv[])
58 {
59  #include <OpenFOAM/setRootCase.H>
60  #include <OpenFOAM/createTime.H>
61  #include <OpenFOAM/createMesh.H>
63  #include "createFields.H"
64 
65  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
66 
67  Info<< "\nStarting time loop\n" << endl;
68 
69  while (runTime.loop())
70  {
71  Info<< "\n Time = " << runTime.timeName() << nl << endl;
72 
74  #include "CourantNo.H"
75 
76  for (int ucorr=0; ucorr<nOuterCorr; ucorr++)
77  {
79 
80  fvVectorMatrix hUEqn
81  (
82  fvm::ddt(hU)
83  + fvm::div(phiv, hU)
84  );
85 
86  hUEqn.relax();
87 
89  {
90  if (rotating)
91  {
92  solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));
93  }
94  else
95  {
96  solve(hUEqn == -magg*h*fvc::grad(h + h0));
97  }
98 
99  // Constrain the momentum to be in the geometry if 3D geometry
100  if (mesh.nGeometricD() == 3)
101  {
102  hU -= (gHat & hU)*gHat;
103  hU.correctBoundaryConditions();
104  }
105  }
106 
107  // --- PISO loop
108  for (int corr=0; corr<nCorr; corr++)
109  {
111  volScalarField rUA = 1.0/hUEqn.A();
112  surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA);
113 
114  surfaceScalarField phih0 = ghrUAf*mesh.magSf()*fvc::snGrad(h0);
115 
116  if (rotating)
117  {
118  hU = rUA*(hUEqn.H() - (F ^ hU));
119  }
120  else
121  {
122  hU = rUA*hUEqn.H();
123  }
124 
125  phi = (fvc::interpolate(hU) & mesh.Sf())
126  + fvc::ddtPhiCorr(rUA, h, hU, phi)
127  - phih0;
128 
129  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
130  {
131  fvScalarMatrix hEqn
132  (
133  fvm::ddt(h)
134  + fvc::div(phi)
135  - fvm::laplacian(ghrUAf, h)
136  );
137 
138  if (ucorr < nOuterCorr-1 || corr < nCorr-1)
139  {
140  hEqn.solve();
141  }
142  else
143  {
144  hEqn.solve(mesh.solver(h.name() + "Final"));
145  }
146 
147  if (nonOrth == nNonOrthCorr)
148  {
149  phi += hEqn.flux();
150  }
151  }
152 
153  hU -= rUA*h*magg*fvc::grad(h + h0);
154 
155  // Constrain the momentum to be in the geometry if 3D geometry
156  if (mesh.nGeometricD() == 3)
157  {
158  hU -= (gHat & hU)*gHat;
159  }
160 
161  hU.correctBoundaryConditions();
162  }
163  }
164 
165  U == hU/h;
166  hTotal == h + h0;
167 
168  runTime.write();
169 
170  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
171  << " ClockTime = " << runTime.elapsedClockTime() << " s"
172  << nl << endl;
173  }
174 
175  Info<< "End\n" << endl;
176 
177  return 0;
178 }
179 
180 
181 // ************************ vim: set sw=4 sts=4 et: ************************ //