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) 2011-2015 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
34 
35 #include "fvCFD.H"
36 #include "pimpleControl.H"
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 int main(int argc, char *argv[])
41 {
42  #include "setRootCase.H"
43  #include "createTime.H"
44  #include "createMesh.H"
45 
46  pimpleControl pimple(mesh);
47 
48  #include "readGravitationalAcceleration.H"
49  #include "createFields.H"
50 
51  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
52 
53  Info<< "\nStarting time loop\n" << endl;
54 
55  while (runTime.loop())
56  {
57  Info<< "\n Time = " << runTime.timeName() << nl << endl;
58 
59  #include "CourantNo.H"
60 
61  // --- Pressure-velocity PIMPLE corrector loop
62  while (pimple.loop())
63  {
65 
66  fvVectorMatrix hUEqn
67  (
68  fvm::ddt(hU)
69  + fvm::div(phiv, hU)
70  );
71 
72  hUEqn.relax();
73 
74  if (pimple.momentumPredictor())
75  {
76  if (rotating)
77  {
78  solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));
79  }
80  else
81  {
82  solve(hUEqn == -magg*h*fvc::grad(h + h0));
83  }
84 
85  // Constrain the momentum to be in the geometry if 3D geometry
86  if (mesh.nGeometricD() == 3)
87  {
88  hU -= (gHat & hU)*gHat;
89  hU.correctBoundaryConditions();
90  }
91  }
92 
93  // --- Pressure corrector loop
94  while (pimple.correct())
95  {
96  volScalarField rAU(1.0/hUEqn.A());
98 
99  surfaceScalarField phih0(ghrAUf*mesh.magSf()*fvc::snGrad(h0));
100 
101  volVectorField HbyA("HbyA", hU);
102  if (rotating)
103  {
104  HbyA = rAU*(hUEqn.H() - (F ^ hU));
105  }
106  else
107  {
108  HbyA = rAU*hUEqn.H();
109  }
110 
112  (
113  "phiHbyA",
114  (fvc::interpolate(HbyA) & mesh.Sf())
116  - phih0
117  );
118 
119  while (pimple.correctNonOrthogonal())
120  {
121  fvScalarMatrix hEqn
122  (
123  fvm::ddt(h)
124  + fvc::div(phiHbyA)
125  - fvm::laplacian(ghrAUf, h)
126  );
127 
128  hEqn.solve(mesh.solver(h.select(pimple.finalInnerIter())));
129 
130  if (pimple.finalNonOrthogonalIter())
131  {
132  phi = phiHbyA + hEqn.flux();
133  }
134  }
135 
136  hU = HbyA - rAU*h*magg*fvc::grad(h + h0);
137 
138  // Constrain the momentum to be in the geometry if 3D geometry
139  if (mesh.nGeometricD() == 3)
140  {
141  hU -= (gHat & hU)*gHat;
142  }
143 
144  hU.correctBoundaryConditions();
145  }
146  }
147 
148  U == hU/h;
149  hTotal == h + h0;
150 
151  runTime.write();
152 
153  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
154  << " ClockTime = " << runTime.elapsedClockTime() << " s"
155  << nl << endl;
156  }
157 
158  Info<< "End\n" << endl;
159 
160  return 0;
161 }
162 
163 
164 // ************************************************************************* //
rhoEqn solve()
phiHbyA
Definition: pcEqn.H:74
tmp< GeometricField< typename flux< Type >::type, fvsPatchField, surfaceMesh > > ddtCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: fvcDdt.C:155
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
dynamicFvMesh & mesh
scalar h0
const dictionary & pimple
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
surfaceScalarField & phi
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
volScalarField rAU(1.0/UEqn.A())
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
int main(int argc, char *argv[])
Definition: postCalc.C:54
const dimensionedScalar F
Faraday constant: default SI units: [C/mol].
U
Definition: pEqn.H:82
const dimensionedScalar h
Planck constant.
Definition: createFields.H:6
HbyA
Definition: pEqn.H:7