applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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     pUCoupledFoam
  26
  27 Description
  28     Steady-state solver for incompressible, turbulent flow, with implicit
  29     coupling between pressure and velocity achieved by BlockLduMatrix
  30     Turbulence is in this version solved using the existing turbulence
  31     structure.
  32
  33 Authors
  34     Klas Jareteg, Chalmers University of Technology,
  35     Vuko Vukcevic, FMENA Zagreb.
  36
  37 \*---------------------------------------------------------------------------*/
  38
  39 #include "fvCFD.H"
  40 #include "singlePhaseTransportModel.H"
  41 #include "RASModel.H"
  42
  43 #include "VectorNFieldTypes.H"
  44 #include "volVectorNFields.H"
  45 #include "blockLduSolvers.H"
  46 #include "blockVectorNMatrices.H"
  47
  48 #include "blockMatrixTools.H"
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 int main(int argc, char *argv[])
  53 {
  54
  55 #   include "setRootCase.H"
  56 #   include "createTime.H"
  57 #   include "createMesh.H"
  58 #   include "createFields.H"
  59 #   include "initContinuityErrs.H"
  60 #   include "readBlockSolverControls.H"
  61
  62     // Calculate coupling matrices only once since the mesh doesn't change and
  63     // implicit div and grad operators are only dependant on Sf. Actually
  64     // coupling terms (div(U) and grad(p)) in blockMatrix do not change, so they
  65     // could be inserted only once, resetting other parts of blockMatrix to zero
  66     // at the end of each time step. VV, 30/April/2014
  67 #   include "calculateCouplingMatrices.H"
  68
  69     Info<< "\nStarting time loop\n" << endl;
  70     while (runTime.loop())
  71     {
  72         Info<< "Time = " << runTime.timeName() << nl << endl;
  73
  74         p.storePrevIter();
  75
  76         // Initialize block matrix
  77 #       include "initializeBlockMatrix.H"
  78
  79         // Assemble and insert momentum equation
  80 #       include "UEqn.H"
  81
  82         // Assemble and insert pressure equation
  83 #       include "pEqn.H"
  84
  85         // Insert coupling, updating the boundary contributions
  86         // Last argument in insertBlockCoupling says if the first location
  87         // should be incremented. This is needed for arbitrary positioning
  88         // of U and p in the system. This could be better. VV, 30/April/2014
  89         blockMatrixTools::insertBlockCoupling(3, 0, UInp, U, A, b, false);
  90         blockMatrixTools::insertBlockCoupling(0, 3, pInU, p, A, b, true);
  91
  92         // Solve the block matrix
  93         BlockSolverPerformance<vector4> solverPerf =
  94             BlockLduSolver<vector4>::New
  95             (
  96                 word("Up"),
  97                 A,
  98                 mesh.solutionDict().solver("Up")
  99             )->solve(Up, b);
 100
 101         solverPerf.print();
 102
 103         // Retrieve solution
 104         blockMatrixTools::retrieveSolution(0, U.internalField(), Up);
 105         blockMatrixTools::retrieveSolution(3, p.internalField(), Up);
 106
 107         U.correctBoundaryConditions();
 108         p.correctBoundaryConditions();
 109
 110         phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
 111
 112 #       include "continuityErrs.H"
 113
 114         p.relax();
 115
 116         turbulence->correct();
 117         runTime.write();
 118
 119         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 120             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 121             << nl << endl;
 122     }
 123
 124     Info<< "End\n" << endl;
 125
 126     return 0;
 127 }