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