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