src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "pressureGradientExplicitSource.H"
  28 #include "volFields.H"
  29 #include "IFstream.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 void Foam::pressureGradientExplicitSource::writeGradP() const
  34 {
  35     // Only write on output time
  36     if (mesh_.time().outputTime())
  37     {
  38         IOdictionary propsDict
  39         (
  40             IOobject
  41             (
  42                 sourceName_ + "Properties",
  43                 mesh_.time().timeName(),
  44                 "uniform",
  45                 mesh_,
  46                 IOobject::NO_READ,
  47                 IOobject::NO_WRITE
  48             )
  49         );
  50         propsDict.add("gradient", gradP_);
  51         propsDict.regIOobject::write();
  52     }
  53 }
  54
  55
  56 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  57
  58 Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
  59 (
  60     const word& sourceName,
  61     volVectorField& U
  62 )
  63 :
  64     sourceName_(sourceName),
  65     mesh_(U.mesh()),
  66     U_(U),
  67     dict_
  68     (
  69         IOobject
  70         (
  71             sourceName + "Properties",
  72             mesh_.time().constant(),
  73             mesh_,
  74             IOobject::MUST_READ,
  75             IOobject::NO_WRITE
  76         )
  77     ),
  78     Ubar_(dict_.lookup("Ubar")),
  79     gradPini_(dict_.lookup("gradPini")),
  80     gradP_(gradPini_),
  81     flowDir_(Ubar_/mag(Ubar_)),
  82     cellSource_(dict_.lookup("cellSource")),
  83     cellSelector_
  84     (
  85         topoSetSource::New
  86         (
  87             cellSource_,
  88             mesh_,
  89             dict_.subDict(cellSource_ + "Coeffs")
  90         )
  91     ),
  92     selectedCellSet_
  93     (
  94         mesh_,
  95         sourceName_ + "CellSet",
  96         mesh_.nCells()/10 + 1  // Reasonable size estimate.
  97     )
  98 {
  99     // Create the cell set
 100     cellSelector_->applyToSet
 101     (
 102         topoSetSource::NEW,
 103         selectedCellSet_
 104     );
 105
 106     // Give some feedback
 107     Info<< "    Selected "
 108         << returnReduce(selectedCellSet_.size(), sumOp<label>())
 109         << " cells" << endl;
 110
 111     // Read the initial pressure gradient from file if it exists
 112     IFstream propsFile
 113     (
 114         mesh_.time().timeName()/"uniform"/(sourceName_ + "Properties")
 115     );
 116
 117     if (propsFile.good())
 118     {
 119         Info<< "    Reading pressure gradient from file" << endl;
 120         dictionary propsDict(dictionary::null, propsFile);
 121         propsDict.lookup("gradient") >> gradP_;
 122     }
 123
 124     Info<< "    Initial pressure gradient = " << gradP_ << nl << endl;
 125 }
 126
 127
 128 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 129
 130 Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
 131 Foam::pressureGradientExplicitSource::Su() const
 132 {
 133     tmp<DimensionedField<vector, volMesh> > tSource
 134     (
 135         new DimensionedField<vector, volMesh>
 136         (
 137             IOobject
 138             (
 139                 sourceName_,
 140                 mesh_.time().timeName(),
 141                 mesh_,
 142                 IOobject::NO_READ,
 143                 IOobject::NO_WRITE
 144             ),
 145             mesh_,
 146             dimensionedVector("zero", gradP_.dimensions(), vector::zero)
 147         )
 148     );
 149
 150     DimensionedField<vector, volMesh>& sourceField = tSource();
 151
 152     forAllConstIter(cellSet, selectedCellSet_, iter)
 153     {
 154         label cellI = iter.key();
 155
 156         sourceField[cellI] = flowDir_*gradP_.value();
 157     }
 158
 159     return tSource;
 160 }
 161
 162
 163 void Foam::pressureGradientExplicitSource::update()
 164 {
 165     const volScalarField& rUA =
 166         mesh_.lookupObject<volScalarField>("(1|A(" + U_.name() + "))");
 167
 168     // Integrate flow variables over cell set
 169     scalar volTot = 0.0;
 170     scalar magUbarAve = 0.0;
 171     scalar rUAave = 0.0;
 172     forAllConstIter(cellSet, selectedCellSet_, iter)
 173     {
 174         label cellI = iter.key();
 175
 176         scalar volCell = mesh_.V()[cellI];
 177         volTot += volCell;
 178
 179         magUbarAve += (flowDir_ & U_[cellI])*volCell;
 180         rUAave += rUA[cellI]*volCell;
 181     }
 182
 183     // Collect across all processors
 184     reduce(volTot, sumOp<scalar>());
 185     reduce(magUbarAve, sumOp<scalar>());
 186     reduce(rUAave, sumOp<scalar>());
 187
 188     // Volume averages
 189     magUbarAve /= volTot;
 190     rUAave /= volTot;
 191
 192     // Calculate the pressure gradient increment needed to adjust the average
 193     // flow-rate to the desired value
 194     scalar gradPplus = (mag(Ubar_) - magUbarAve)/rUAave;
 195
 196     // Apply correction to velocity field
 197     forAllConstIter(cellSet, selectedCellSet_, iter)
 198     {
 199         label cellI = iter.key();
 200         U_[cellI] += flowDir_*rUA[cellI]*gradPplus;
 201     }
 202
 203     // Update pressure gradient
 204     gradP_.value() += gradPplus;
 205
 206     Info<< "Uncorrected Ubar = " << magUbarAve << tab
 207         << "Pressure gradient = " << gradP_.value() << endl;
 208
 209     writeGradP();
 210 }
 211
 212
 213 // ************************************************************************* //