src/POD/PODOrthoNormalBase/PODOrthoNormalBase.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     PODOrthoNormalBase
  27
  28 Description
  29     Establish POD ortho-normal base and interpolation coefficients give a list
  30     of fields. Size of ortho-normal base is calculated from the desired
  31     accuracy, e.g. 0.99-0.99999 (in energy terms)
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "PODOrthoNormalBase.H"
  36 #include "POD.H"
  37 #include "PODEigenBase.H"
  38 #include "IOmanip.H"
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 template<class Type>
  43 void Foam::PODOrthoNormalBase<Type>::calcOrthoBase
  44 (
  45     const PtrList<GeometricField<Type, fvPatchField, volMesh> >& snapshots,
  46     const scalar accuracy
  47 )
  48 {
  49     // Calculate ortho-normal base for each component
  50     PtrList<PODEigenBase> eigenBase(pTraits<Type>::nComponents);
  51
  52     const label nSnapshots = snapshots.size();
  53
  54     typename
  55     powProduct<Vector<label>, pTraits<Type>::rank>::type validComponents
  56     (
  57         pow
  58         (
  59             snapshots[0].mesh().solutionD(),
  60             pTraits
  61             <
  62                 typename powProduct<Vector<label>,
  63                 pTraits<Type>::rank
  64             >::type>::zero
  65         )
  66     );
  67
  68     label nValidCmpts = 0;
  69
  70     for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
  71     {
  72         // Component not valid, skipping
  73         if (validComponents[cmpt] == -1) continue;
  74
  75         // Create a list of snapshots
  76         PtrList<volScalarField> sf(nSnapshots);
  77
  78         forAll (snapshots, i)
  79         {
  80             sf.set(i, new volScalarField(snapshots[i].component(cmpt)));
  81         }
  82
  83         // Create eigen base
  84         eigenBase.set(cmpt, new PODEigenBase(sf));
  85
  86         Info<< "Cumulative eigen-values for component " << cmpt
  87             << ": " << setprecision(14)
  88             << eigenBase[nValidCmpts].cumulativeEigenValues() << endl;
  89
  90         nValidCmpts++;
  91     }
  92
  93     eigenBase.setSize(nValidCmpts);
  94
  95     Info << "Number of valid eigen components: " << nValidCmpts << endl;
  96
  97     label baseSize = 0;
  98     for (label snapI = 0; snapI < nSnapshots; snapI++)
  99     {
 100         baseSize++;
 101
 102         // Get minimum cumulative eigen value for all valid components
 103         scalar minCumEigen = 1.0;
 104
 105         nValidCmpts = 0;
 106
 107         for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
 108         {
 109             // Skip invalid components
 110             if (validComponents[cmpt] != -1)
 111             {
 112                 minCumEigen =
 113                     Foam::min
 114                     (
 115                         minCumEigen,
 116                         eigenBase[nValidCmpts].cumulativeEigenValues()[snapI]
 117                     );
 118
 119                 nValidCmpts++;
 120             }
 121         }
 122
 123         if (minCumEigen > accuracy)
 124         {
 125             break;
 126         }
 127     }
 128
 129     Info << "Base size: " << baseSize << endl;
 130
 131     // Establish orthonormal base
 132     orthoFields_.setSize(baseSize);
 133
 134     for (label baseI = 0; baseI < baseSize; baseI++)
 135     {
 136         GeometricField<Type, fvPatchField, volMesh>* onBasePtr
 137         (
 138             new GeometricField<Type, fvPatchField, volMesh>
 139             (
 140                 IOobject
 141                 (
 142                     snapshots[0].name() + "POD" + name(baseI),
 143                     snapshots[0].time().timeName(),
 144                     snapshots[0].mesh(),
 145                     IOobject::NO_READ,
 146                     IOobject::NO_WRITE
 147                 ),
 148                 snapshots[0].mesh(),
 149                 dimensioned<Type>
 150                 (
 151                     "zero",
 152                     snapshots[0].dimensions(),
 153                     pTraits<Type>::zero
 154                 )
 155             )
 156         );
 157         GeometricField<Type, fvPatchField, volMesh>& onBase = *onBasePtr;
 158
 159         nValidCmpts = 0;
 160
 161         for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
 162         {
 163             if (validComponents[cmpt] != -1)
 164             {
 165                 // Valid component, grab eigen-factors
 166
 167                 const scalarField& eigenVector =
 168                     eigenBase[nValidCmpts].eigenVectors()[baseI];
 169                 nValidCmpts++;
 170
 171                 volScalarField onBaseCmpt = onBase.component(cmpt);
 172
 173                 forAll (eigenVector, eigenI)
 174                 {
 175                     onBaseCmpt +=
 176                         eigenVector[eigenI]*snapshots[eigenI].component(cmpt);
 177                 }
 178
 179                 // Re-normalise ortho-normal vector
 180                 onBaseCmpt /= Foam::sqrt(sumSqr(onBaseCmpt));
 181
 182                 onBase.replace(cmpt, onBaseCmpt);
 183             }
 184             else
 185             {
 186                 // Component invalid.  Grab first snapshot.
 187                 onBase.replace
 188                 (
 189                     cmpt,
 190                     snapshots[0].component(cmpt)
 191                 );
 192             }
 193         }
 194
 195         orthoFields_.set(baseI, onBasePtr);
 196     }
 197
 198     // Calculate interpolation coefficients
 199     interpolationCoeffsPtr_ =
 200         new RectangularMatrix<Type>(snapshots.size(), orthoFields_.size());
 201     RectangularMatrix<Type>& coeffs = *interpolationCoeffsPtr_;
 202
 203     forAll (snapshots, snapshotI)
 204     {
 205         forAll (orthoFields_, baseI)
 206         {
 207             coeffs[snapshotI][baseI] =
 208                 POD::projection
 209                 (
 210                     snapshots[snapshotI],
 211                     orthoFields_[baseI]
 212                 );
 213         }
 214     }
 215 }
 216
 217
 218 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 219
 220 // given list of snapshots and accuracy
 221 template<class Type>
 222 Foam::PODOrthoNormalBase<Type>::PODOrthoNormalBase
 223 (
 224     const PtrList<GeometricField<Type, fvPatchField, volMesh> >& snapshots,
 225     const scalar accuracy
 226 )
 227 :
 228     orthoFields_(),
 229     interpolationCoeffsPtr_(NULL)
 230 {
 231     calcOrthoBase(snapshots, accuracy);
 232 }
 233
 234
 235 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 236
 237 template<class Type>
 238 Foam::PODOrthoNormalBase<Type>::~PODOrthoNormalBase()
 239 {
 240     deleteDemandDrivenData(interpolationCoeffsPtr_);
 241 }
 242
 243
 244 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 245
 246
 247 // ************************************************************************* //