incrementalPOD.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * In real Time Highly Advanced Computational Applications for Finite Volumes
 * Copyright (C) 2017 by the ITHACA-FV authors
-------------------------------------------------------------------------------
License
    This file is part of ITHACA-FV
    ITHACA-FV is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    ITHACA-FV is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU Lesser General Public License for more details.
    You should have received a copy of the GNU Lesser General Public License
    along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
Class
    incrementalPOD
Description
    Implementation of a incremental POD decomposition of a general field
SourceFiles
    incrementalPOD.C
\*---------------------------------------------------------------------------*/

 
#ifndef INCREMENTALPOD_H
#define INCREMENTALPOD_H
 
#include "fvCFD.H"
#include "ITHACAutilities.H"
#include "ITHACAstream.H"
#include "ITHACAparameters.H"
#include "Foam2Eigen.H"
#include "EigenFunctions.H"
#include "Modes.H"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
#pragma GCC diagnostic ignored "-Wnon-virtual-dtor"
#include <Spectra/GenEigsSolver.h>
#include <Spectra/SymEigsSolver.h>
#include <Eigen/Eigen>
#include <unsupported/Eigen/SparseExtra>
#pragma GCC diagnostic pop
 
/*---------------------------------------------------------------------------*\
                        Class incrementalPOD Declaration
\*---------------------------------------------------------------------------*/

 
template<class Type, template<class> class PatchField, class GeoMesh>
class incrementalPOD : public Modes<Type, PatchField, GeoMesh>
{
    public:

        double tolleranceSVD = 0;

        int rank = 0;

        word PODnorm;

        Eigen::VectorXd massVector;

        Eigen::VectorXd singularValues;

        word outputFolder = "./ITHACAoutput/incrementalPOD/";


        incrementalPOD() {};
 
        //--------------------------------------------------------------------------
        incrementalPOD(GeometricField<Type, PatchField, GeoMesh>& snapshot,
                       double _tol, word _PODnorm);

        ~incrementalPOD() {};
 
        //--------------------------------------------------------------------------
        void initialize(GeometricField<Type, PatchField, GeoMesh>& snapshot);
 
        //--------------------------------------------------------------------------
        void addSnapshot(GeometricField<Type, PatchField, GeoMesh>& snapshot);
 
        //--------------------------------------------------------------------------
        void fillPtrList();
 
        //--------------------------------------------------------------------------
        Eigen::VectorXd project(
            GeometricField<Type, PatchField, GeoMesh>& inputField,
            label numberOfModes = 0);
 
        //--------------------------------------------------------------------------
        GeometricField<Type, PatchField, GeoMesh> reconstruct(
            GeometricField<Type, PatchField, GeoMesh>& inputField,
            Eigen::MatrixXd Coeff,
            word Name);
 
        //--------------------------------------------------------------------------
        GeometricField<Type, PatchField, GeoMesh> projectSnapshot(
            GeometricField<Type, PatchField, GeoMesh>& snapshot,
            label numberOfModes = 0);
 
        //--------------------------------------------------------------------------
        void projectSnapshots(
            PtrList<GeometricField<Type, PatchField, GeoMesh >> snapshots,
            PtrList<GeometricField<Type, PatchField, GeoMesh >>& projSnapshots,
            label numberOfModes = 0);
 
        //--------------------------------------------------------------------------
        void writeModes();
};
 
typedef incrementalPOD<scalar, fvPatchField, volMesh> scalarIncrementalPOD;
typedef incrementalPOD<vector, fvPatchField, volMesh> vectorIncrementalPOD;
 
namespace ITHACAPOD
{
//--------------------------------------------------------------------------
void weightedGramSchmidt(
    Eigen::MatrixXd& matrix,
    Eigen::VectorXd& weights)
{
    M_Assert(matrix.rows() == weights.rows(),
             "Matrix and weights must have the same number of rows");
    Eigen::MatrixXd Ortho = matrix;
    Ortho = matrix;
 
    for (label i = 0; i < matrix.cols(); i++)
    {
        for (label k = 0; k < i; k++)
        {
            double num = Ortho.col(k).transpose() * weights.asDiagonal()
                         * matrix.col(i);
            double den = (Ortho.col(k).transpose() * weights.asDiagonal()
                          * Ortho.col(k));
            double fact = num / den;
            Ortho.col(i) -= fact * Ortho.col(k) ;
        }
 
        Ortho.col(i).normalize();
    }
 
    matrix = Ortho;
};
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif