applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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 \*---------------------------------------------------------------------------*/
  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     #include "readGravitationalAcceleration.H"
  46     #include "createFields.H"
  47
  48     pimpleControl pimple(mesh);
  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         for (pimple.start(); pimple.loop(); pimple++)
  62         {
  63             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
  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             // --- PISO loop
  93             for (int corr=0; corr<pimple.nCorr(); corr++)
  94             {
  95                 volScalarField rAU(1.0/hUEqn.A());
  96                 surfaceScalarField ghrAUf(magg*fvc::interpolate(h*rAU));
  97
  98                 surfaceScalarField phih0(ghrAUf*mesh.magSf()*fvc::snGrad(h0));
  99
 100                 if (rotating)
 101                 {
 102                     hU = rAU*(hUEqn.H() - (F ^ hU));
 103                 }
 104                 else
 105                 {
 106                     hU = rAU*hUEqn.H();
 107                 }
 108
 109                 phi = (fvc::interpolate(hU) & mesh.Sf())
 110                     + fvc::ddtPhiCorr(rAU, h, hU, phi)
 111                     - phih0;
 112
 113                 for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 114                 {
 115                     fvScalarMatrix hEqn
 116                     (
 117                         fvm::ddt(h)
 118                       + fvc::div(phi)
 119                       - fvm::laplacian(ghrAUf, h)
 120                     );
 121
 122                     hEqn.solve
 123                     (
 124                         mesh.solver
 125                         (
 126                             h.select(pimple.finalInnerIter(corr, nonOrth))
 127                         )
 128                     );
 129
 130                     if (nonOrth == pimple.nNonOrthCorr())
 131                     {
 132                         phi += hEqn.flux();
 133                     }
 134                 }
 135
 136                 hU -= 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 // ************************************************************************* //