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