src/POD/scalarTransportPOD/scalarTransportPOD.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Class
  26     scalarTransportPOD
  27
  28 Description
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "scalarTransportPOD.H"
  33 #include "fvc.H"
  34 #include "addToRunTimeSelectionTable.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40     defineTypeNameAndDebug(scalarTransportPOD, 0);
  41
  42     addToRunTimeSelectionTable
  43     (
  44         PODODE,
  45         scalarTransportPOD,
  46         dictionary
  47     );
  48 }
  49
  50
  51 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  52
  53 void Foam::scalarTransportPOD::calcOrthoBase() const
  54 {
  55     if (orthoBasePtr_)
  56     {
  57         FatalErrorIn
  58         (
  59             "scalarTransportPOD::calcOrthoBase()"
  60         )   << "Orthogonal base already calculated"
  61             << abort(FatalError);
  62     }
  63
  64     // Create ortho-normal base
  65     scalar accuracy = readScalar(dict().lookup("accuracy"));
  66
  67     // Get times list
  68     Time& runTime = const_cast<Time&>(mesh().time());
  69
  70     // Remember time index to restore it after the scan
  71     label origTimeIndex = runTime.timeIndex();
  72
  73     instantList Times = runTime.times();
  74
  75     // Create a list of snapshots
  76     PtrList<volScalarField> fields(Times.size());
  77
  78     label nSnapshots = 0;
  79
  80     forAll (Times, i)
  81     {
  82         if (Times[i].value() < SMALL || Times[i] == runTime.constant())
  83         {
  84             Info << "Skipping time " << Times[i] << endl;
  85
  86             continue;
  87         }
  88
  89         runTime.setTime(Times[i], i);
  90
  91         Info<< "Time = " << runTime.timeName() << endl;
  92
  93         IOobject phiHeader
  94         (
  95             phiName_,
  96             runTime.timeName(),
  97             mesh(),
  98             IOobject::MUST_READ
  99         );
 100
 101         if (phiHeader.headerOk())
 102         {
 103             Info<< "    Reading " << phiName_ << endl;
 104
 105             fields.set(nSnapshots, new volScalarField(phiHeader, mesh()));
 106
 107             // Rename the field
 108             fields[nSnapshots].rename(phiName_ + name(i));
 109             nSnapshots++;
 110         }
 111         else
 112         {
 113             Info<< "    No " << phiName_ << endl;
 114         }
 115     }
 116
 117     // Reset time index to initial state
 118     runTime.setTime(Times[origTimeIndex], origTimeIndex);
 119
 120     // Resize snapshots
 121     if (nSnapshots < 2)
 122     {
 123         FatalErrorIn
 124         (
 125             "scalarTransportPOD::calcOrthoBase()"
 126         )   << "Insufficient number of snapshots: " << nSnapshots
 127             << abort(FatalError);
 128     }
 129
 130     Info << "Number of snapshots: " << nSnapshots << endl;
 131
 132     fields.setSize(nSnapshots);
 133
 134     // Create ortho-normal base for transported variable
 135     orthoBasePtr_ = new scalarPODOrthoNormalBase(fields, accuracy);
 136 }
 137
 138
 139 void Foam::scalarTransportPOD::calcDerivativeCoeffs() const
 140 {
 141     if (derivativeMatrixPtr_)
 142     {
 143         FatalErrorIn
 144         (
 145             "void scalarTransportPOD::calcDerivativeCoeffs() const"
 146         )   << "Derivative matrix already calculated"
 147             << abort(FatalError);
 148     }
 149
 150     // Calculate coefficients for differential equation
 151     // Get times list
 152     Time& runTime = const_cast<Time&>(this->mesh().time());
 153
 154     // Remember time index to restore it
 155     label origTimeIndex = runTime.timeIndex();
 156
 157     // Read diffusivity
 158
 159     IOdictionary transportProperties
 160     (
 161         IOobject
 162         (
 163             "transportProperties",
 164             runTime.constant(),
 165             this->mesh(),
 166             IOobject::MUST_READ,
 167             IOobject::NO_WRITE
 168         )
 169     );
 170
 171     Info<< "Reading diffusivity D\n" << endl;
 172
 173     dimensionedScalar DT
 174     (
 175         transportProperties.lookup("DT")
 176     );
 177
 178     // Find velocity field
 179
 180     word Uname(this->dict().lookup("velocity"));
 181
 182     instantList Times = runTime.times();
 183
 184     volVectorField* Uptr = NULL;
 185
 186     forAll (Times, i)
 187     {
 188         runTime.setTime(Times[i], i);
 189
 190         Info<< "Time = " << runTime.timeName() << endl;
 191
 192         IOobject Uheader
 193         (
 194             Uname,
 195             runTime.timeName(),
 196             this->mesh(),
 197             IOobject::MUST_READ
 198         );
 199
 200         if (Uheader.headerOk())
 201         {
 202             Info<< "    Reading " << Uname << endl;
 203
 204             Uptr = new volVectorField(Uheader, this->mesh());
 205             break;
 206         }
 207         else
 208         {
 209             Info<< "    No " << Uname << endl;
 210         }
 211     }
 212
 213     // Reset time index to initial state
 214     runTime.setTime(Times[origTimeIndex], origTimeIndex);
 215
 216     if (!Uptr)
 217     {
 218         FatalErrorIn
 219         (
 220             "void scalarTransportPOD::calcDerivativeCoeffs() const"
 221         )   << "Cannot find velocity: " << Uname
 222             << abort(FatalError);
 223     }
 224
 225     volVectorField& U = *Uptr;
 226
 227     // Create derivative matrix
 228
 229     const scalarPODOrthoNormalBase& b = orthoBase();
 230
 231     derivativeMatrixPtr_ = new scalarSquareMatrix(b.baseSize(), 0.0);
 232     scalarSquareMatrix& derivative = *derivativeMatrixPtr_;
 233
 234     for (label i = 0; i < b.baseSize(); i++)
 235     {
 236         const volScalarField& snapI = b.orthoField(i);
 237
 238         volVectorField gradSnapI = fvc::grad(snapI);
 239
 240         for (label j = 0; j < b.baseSize(); j++)
 241         {
 242             const volScalarField& snapJ = b.orthoField(j);
 243
 244             volVectorField gradSnapJ = fvc::grad(snapJ);
 245
 246             derivative[i][j] =
 247                 DT.value()*POD::projection(fvc::div(gradSnapJ), snapI)
 248               - POD::projection((U & gradSnapJ), snapI);
 249         }
 250     }
 251 }
 252
 253 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 254
 255 // Construct from components
 256 Foam::scalarTransportPOD::scalarTransportPOD
 257 (
 258     const fvMesh& mesh,
 259     const dictionary& dict
 260 )
 261 :
 262     PODODE(mesh, dict),
 263     phiName_(dict.lookup("field")),
 264     coeffs_(),
 265     derivativeMatrixPtr_(NULL),
 266     orthoBasePtr_(NULL),
 267     fieldPtr_(NULL)
 268 {
 269     // Grab coefficients from the first snapshot of the ortho-normal base
 270     coeffs_.setSize(orthoBase().baseSize());
 271
 272     const scalarRectangularMatrix& orthoBaseCoeffs =
 273         orthoBase().interpolationCoeffs();
 274
 275     forAll (coeffs_, i)
 276     {
 277         coeffs_[i] = orthoBaseCoeffs[0][i];
 278     }
 279 }
 280
 281
 282 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 283
 284 Foam::scalarTransportPOD::~scalarTransportPOD()
 285 {
 286     deleteDemandDrivenData(derivativeMatrixPtr_);
 287
 288     clearBase();
 289     clearFields();
 290 }
 291
 292
 293 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 294
 295 Foam::label Foam::scalarTransportPOD::nEqns() const
 296 {
 297     return coeffs().size();
 298 }
 299
 300
 301 Foam::scalarField& Foam::scalarTransportPOD::coeffs()
 302 {
 303     return coeffs_;
 304 }
 305
 306
 307 const Foam::scalarField& Foam::scalarTransportPOD::coeffs() const
 308 {
 309     return coeffs_;
 310 }
 311
 312
 313 void Foam::scalarTransportPOD::derivatives
 314 (
 315     const scalar x,
 316     const scalarField& y,
 317     scalarField& dydx
 318 ) const
 319 {
 320     if (!derivativeMatrixPtr_)
 321     {
 322         calcDerivativeCoeffs();
 323     }
 324
 325    const scalarSquareMatrix& derivative = *derivativeMatrixPtr_;
 326
 327     forAll (dydx, i)
 328     {
 329         dydx[i] = 0;
 330
 331         forAll (y, j)
 332         {
 333             dydx[i] += derivative[i][j]*y[j];
 334         }
 335     }
 336 }
 337
 338
 339 void Foam::scalarTransportPOD::jacobian
 340 (
 341     const scalar x,
 342     const scalarField& y,
 343     scalarField& dfdx,
 344     scalarSquareMatrix& dfdy
 345 ) const
 346 {
 347     derivatives(x, y, dfdx);
 348     dfdy = 0;
 349 }
 350
 351
 352 const Foam::scalarPODOrthoNormalBase&
 353 Foam::scalarTransportPOD::orthoBase() const
 354 {
 355     if (!orthoBasePtr_)
 356     {
 357         calcOrthoBase();
 358     }
 359
 360     return *orthoBasePtr_;
 361 }
 362
 363
 364 void Foam::scalarTransportPOD::clearBase() const
 365 {
 366     deleteDemandDrivenData(orthoBasePtr_);
 367 }
 368
 369
 370 const Foam::volScalarField& Foam::scalarTransportPOD::field() const
 371 {
 372     if (!fieldPtr_)
 373     {
 374         updateFields();
 375     }
 376
 377     return *fieldPtr_;
 378 }
 379
 380
 381 void Foam::scalarTransportPOD::updateFields() const
 382 {
 383     if (!fieldPtr_)
 384     {
 385         // Allocate field
 386         fieldPtr_ =
 387             new volScalarField
 388             (
 389                 IOobject
 390                 (
 391                     phiName_ + "POD",
 392                     mesh().time().timeName(),
 393                     mesh(),
 394                     IOobject::NO_READ,
 395                     IOobject::NO_WRITE
 396                 ),
 397                 mesh(),
 398                 dimensionedScalar
 399                 (
 400                     "zero",
 401                     orthoBase().orthoField(0).dimensions(),
 402                     0
 403                 )
 404             );
 405     }
 406
 407     volScalarField& phi = *fieldPtr_;
 408
 409     const scalarPODOrthoNormalBase& b = orthoBase();
 410
 411     phi = dimensionedScalar("zero", b.orthoField(0).dimensions(), 0);
 412
 413     forAll (coeffs_, i)
 414     {
 415         phi += coeffs_[i]*b.orthoField(i);
 416     }
 417 }
 418
 419
 420 void Foam::scalarTransportPOD::clearFields() const
 421 {
 422     deleteDemandDrivenData(fieldPtr_);
 423 }
 424
 425
 426 void Foam::scalarTransportPOD::write() const
 427 {
 428     // Recalculate field and force a write
 429     updateFields();
 430     field().write();
 431 }
 432
 433
 434 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 435
 436
 437 // ************************************************************************* //