applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  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
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 int main(int argc, char *argv[])
  40 {
  41     #include "setRootCase.H"
  42     #include "createTime.H"
  43     #include "createMesh.H"
  44     #include "readGravitationalAcceleration.H"
  45     #include "createFields.H"
  46
  47     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48
  49     Info<< "\nStarting time loop\n" << endl;
  50
  51     while (runTime.loop())
  52     {
  53         Info<< "\n Time = " << runTime.timeName() << nl << endl;
  54
  55         #include "readPISOControls.H"
  56         #include "CourantNo.H"
  57
  58         for (int ucorr=0; ucorr<nOuterCorr; ucorr++)
  59         {
  60             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
  61
  62             fvVectorMatrix hUEqn
  63             (
  64                 fvm::ddt(hU)
  65               + fvm::div(phiv, hU)
  66             );
  67
  68             hUEqn.relax();
  69
  70             if (momentumPredictor)
  71             {
  72                 if (rotating)
  73                 {
  74                     solve(hUEqn + (F ^ hU) == -magg*h*fvc::grad(h + h0));
  75                 }
  76                 else
  77                 {
  78                     solve(hUEqn == -magg*h*fvc::grad(h + h0));
  79                 }
  80
  81                 // Constrain the momentum to be in the geometry if 3D geometry
  82                 if (mesh.nGeometricD() == 3)
  83                 {
  84                     hU -= (gHat & hU)*gHat;
  85                     hU.correctBoundaryConditions();
  86                 }
  87             }
  88
  89             // --- PISO loop
  90             for (int corr = 0; corr < nCorr; corr++)
  91             {
  92                 surfaceScalarField hf = fvc::interpolate(h);
  93                 volScalarField rUA = 1.0/hUEqn.A();
  94                 surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA);
  95
  96                 surfaceScalarField phih0 = ghrUAf*mesh.magSf()*fvc::snGrad(h0);
  97
  98                 if (rotating)
  99                 {
 100                     hU = rUA*(hUEqn.H() - (F ^ hU));
 101                 }
 102                 else
 103                 {
 104                     hU = rUA*hUEqn.H();
 105                 }
 106
 107                 phi = (fvc::interpolate(hU) & mesh.Sf())
 108                     + fvc::ddtPhiCorr(rUA, h, hU, phi)
 109                     - phih0;
 110
 111                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 112                 {
 113                     fvScalarMatrix hEqn
 114                     (
 115                         fvm::ddt(h)
 116                       + fvc::div(phi)
 117                       - fvm::laplacian(ghrUAf, h)
 118                     );
 119
 120                     if (ucorr < nOuterCorr - 1 || corr < nCorr - 1)
 121                     {
 122                         hEqn.solve();
 123                     }
 124                     else
 125                     {
 126                         hEqn.solve
 127                         (
 128                             mesh.solutionDict().solver(h.name() + "Final")
 129                         );
 130                     }
 131
 132                     if (nonOrth == nNonOrthCorr)
 133                     {
 134                         phi += hEqn.flux();
 135                     }
 136                 }
 137
 138                 hU -= rUA*h*magg*fvc::grad(h + h0);
 139
 140                 // Constrain the momentum to be in the geometry if 3D geometry
 141                 if (mesh.nGeometricD() == 3)
 142                 {
 143                     hU -= (gHat & hU)*gHat;
 144                 }
 145
 146                 hU.correctBoundaryConditions();
 147             }
 148         }
 149
 150         U == hU/h;
 151         hTotal == h + h0;
 152
 153         runTime.write();
 154
 155         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 156             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 157             << nl << endl;
 158     }
 159
 160     Info<< "End\n" << endl;
 161
 162     return 0;
 163 }
 164
 165
 166 // ************************************************************************* //