src/POD/scalarTransportPOD/scalarTransportPOD.C

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