src/POD/PODEigenBase/PODEigenBase.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     PODEigenBase
  27
  28 Description
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "PODEigenBase.H"
  33 #include "volFields.H"
  34 #include "POD.H"
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 void Foam::PODEigenBase::calcEigenBase(const scalarSquareMatrix& orthMatrix)
  39 {
  40     // Calculate eigen-values
  41
  42     EigenSolver<scalar> eigenSolver(orthMatrix);
  43
  44     // Sort and assemble
  45
  46     SortableList<scalar> sortedList(orthMatrix.n());
  47
  48     forAll (sortedList, i)
  49     {
  50         sortedList[i] = eigenSolver.eigenValue(i);
  51     }
  52
  53     // Sort will sort the values in descending order and insert values
  54     sortedList.sort();
  55
  56     label n = 0;
  57     forAllReverse(sortedList, i)
  58     {
  59         eigenValues_[n] = sortedList[i];
  60         eigenVectors_.set
  61         (
  62             n,
  63             new scalarField(eigenSolver.eigenVector(sortedList.indices()[i]))
  64         );
  65
  66         n++;
  67     }
  68
  69     // Assemble cumulative relative eigen-values
  70     cumEigenValues_[0] = eigenValues_[0];
  71
  72     // Assemble accumulated normalised eigenvalues
  73     for (label i = 1; i < cumEigenValues_.size(); i++)
  74     {
  75         cumEigenValues_[i] = cumEigenValues_[i - 1] + eigenValues_[i];
  76     }
  77
  78     // Renormalise
  79     cumEigenValues_ /= sum(eigenValues_);
  80
  81 //     // Check products
  82 //     for (label i = 0; i < orthMatrix.m(); i++)
  83 //     {
  84 //         const scalarField& eVector = eigenVectors_[i];
  85
  86 //         Info<< "Scalar product: "
  87 //             << eigenValues_[i]*eVector
  88 //             << endl;
  89
  90 //         scalarField vp(orthMatrix.m(), 0);
  91
  92 //         forAll (vp, i)
  93 //         {
  94 //             forAll (vp, j)
  95 //             {
  96 //                 vp[i] += orthMatrix[i][j]*eVector[j];
  97 //             }
  98 //         }
  99
 100 //         Info << "Vector product: " << vp << nl << endl;
 101 //     }
 102 }
 103
 104
 105 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 106
 107 // Construct given a list of fields
 108 Foam::PODEigenBase::PODEigenBase(const PtrList<volScalarField>& snapshots)
 109 :
 110     eigenValues_(snapshots.size()),
 111     cumEigenValues_(snapshots.size()),
 112     eigenVectors_(snapshots.size())
 113 {
 114     // Calculate the snapshot of the field with all available fields
 115     label nSnapshots = snapshots.size();
 116
 117     scalarSquareMatrix orthMatrix(nSnapshots);
 118
 119     for (label snapI = 0; snapI < nSnapshots; snapI++)
 120     {
 121         for (label snapJ = 0; snapJ <= snapI; snapJ++)
 122         {
 123             // Calculate the inner product and insert it into the matrix
 124             orthMatrix[snapI][snapJ] =
 125                 POD::projection
 126                 (
 127                     snapshots[snapI],
 128                     snapshots[snapJ]
 129                 );
 130
 131             if (snapI != snapJ)
 132             {
 133                 orthMatrix[snapJ][snapI] = orthMatrix[snapI][snapJ];
 134
 135 //                 Info << "Product: " << orthMatrix[snapI][snapJ]
 136 //                     << " relative: "
 137 //                     <<
 138 //                     orthMatrix[snapI][snapJ]/
 139 //                     (
 140 //                         Foam::sqrt(sumSqr(snapshots[snapI]))*
 141 //                         Foam::sqrt(sumSqr(snapshots[snapJ]))
 142 //                       + SMALL
 143 //                     ) << endl;
 144             }
 145         }
 146     }
 147
 148     calcEigenBase(orthMatrix);
 149 }
 150
 151
 152 // ************************************************************************* //