shallowWaterFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 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 "argList.H"
36 #include "timeSelector.H"
37 #include "pimpleControl.H"
38 
39 #include "fvcDdt.H"
40 #include "fvcGrad.H"
41 #include "fvcSnGrad.H"
42 #include "fvcFlux.H"
43 
44 #include "fvmDdt.H"
45 #include "fvmDiv.H"
46 #include "fvmLaplacian.H"
47 
48 using namespace Foam;
49 
50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
51 
52 int main(int argc, char *argv[])
53 {
54  #include "postProcess.H"
55 
56  #include "setRootCase.H"
57  #include "createTime.H"
58  #include "createMesh.H"
59  #include "createControl.H"
60  #include "createFields.H"
61 
62  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
63 
64  Info<< "\nStarting time loop\n" << endl;
65 
66  while (pimple.loop(runTime))
67  {
68  Info<< "\n Time = " << runTime.name() << nl << endl;
69 
70  #include "CourantNo.H"
71 
72  // --- Pressure-velocity PIMPLE corrector loop
73  while (pimple.loop())
74  {
75  surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
76 
77  fvVectorMatrix hUEqn
78  (
79  fvm::ddt(hU)
80  + fvm::div(phiv, hU)
81  );
82 
83  hUEqn.relax();
84 
86  {
87  if (rotating)
88  {
89  solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));
90  }
91  else
92  {
93  solve(hUEqn == -magg*h*fvc::grad(h + h0));
94  }
95 
96  // Constrain the momentum to be in the geometry if 3D geometry
97  if (mesh.nGeometricD() == 3)
98  {
99  hU -= (gHat & hU)*gHat;
100  hU.correctBoundaryConditions();
101  }
102  }
103 
104  // --- Pressure corrector loop
105  while (pimple.correct())
106  {
107  volScalarField rAU(1.0/hUEqn.A());
109 
110  surfaceScalarField phih0(ghrAUf*mesh.magSf()*fvc::snGrad(h0));
111 
112  volVectorField HbyA("HbyA", hU);
113  if (rotating)
114  {
115  HbyA = rAU*(hUEqn.H() - (F ^ hU));
116  }
117  else
118  {
119  HbyA = rAU*hUEqn.H();
120  }
121 
123  (
124  "phiHbyA",
125  fvc::flux(HbyA)
126  + fvc::interpolate(rAU)*fvc::ddtCorr(h, hU, phi)
127  - phih0
128  );
129 
130  while (pimple.correctNonOrthogonal())
131  {
132  fvScalarMatrix hEqn
133  (
134  fvm::ddt(h)
135  + fvc::div(phiHbyA)
136  - fvm::laplacian(ghrAUf, h)
137  );
138 
139  hEqn.solve();
140 
142  {
143  phi = phiHbyA + hEqn.flux();
144  }
145  }
146 
147  hU = HbyA - rAU*h*magg*fvc::grad(h + h0);
148 
149  // Constrain the momentum to be in the geometry if 3D geometry
150  if (mesh.nGeometricD() == 3)
151  {
152  hU -= (gHat & hU)*gHat;
153  }
154 
155  hU.correctBoundaryConditions();
156  }
157  }
158 
159  U == hU/h;
160  hTotal == h + h0;
161 
162  runTime.write();
163 
164  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
165  << " ClockTime = " << runTime.elapsedClockTime() << " s"
166  << nl << endl;
167  }
168 
169  Info<< "End\n" << endl;
170 
171  return 0;
172 }
173 
174 
175 // ************************************************************************* //
Calculates and outputs the maximum Courant Number.
Generic GeometricField class.
bool momentumPredictor() const
Flag to indicate to solve for momentum.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:603
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
tmp< volScalarField > A() const
Return the central coefficient.
Definition: fvMatrix.C:808
tmp< VolField< Type > > H() const
Return the H operation source.
Definition: fvMatrix.C:868
tmp< SurfaceField< Type > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:964
bool finalNonOrthogonalIter() const
Flag to indicate the last non-orthogonal iteration.
bool correct()
Piso loop within outer loop.
bool correctNonOrthogonal()
Non-orthogonal corrector loop.
pimpleControl pimple(mesh)
Calculate the first temporal derivative.
Calculate the face-flux of the given field.
Calculate the gradient of the given field.
Calculate the snGrad of the given volField.
Calculate the matrix for the first temporal derivative.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
volScalarField rAU(1.0/UEqn.A())
U
Definition: pEqn.H:72
volVectorField & HbyA
Definition: pEqn.H:13
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho) *fvc::flux(HbyA))
const dimensionedScalar F
Faraday constant: default SI units: [C/kmol].
const dimensionedScalar h
Planck constant.
tmp< SurfaceField< typename innerProduct< vector, Type >::type > > flux(const VolField< Type > &vf)
Return the face-flux field obtained from the given volVectorField.
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
tmp< SurfaceField< Type > > snGrad(const VolField< Type > &vf, const word &name)
Definition: fvcSnGrad.C:45
tmp< SurfaceField< typename Foam::flux< Type >::type > > ddtCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
Definition: fvcDdt.C:196
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
Namespace for OpenFOAM.
scalar h0(const fluidMulticomponentThermo &thermo, const scalarList &Y, const scalar p, const scalar T)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
messageStream Info
static const char nl
Definition: Ostream.H:266
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
Execute application functionObjects to post-process existing results.
int main(int argc, char *argv[])