hyperReduction.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
    HyperReduction
Description
    Implementation of the hyper-reduction methods.
SourceFiles
    hyperReduction.C
\*---------------------------------------------------------------------------*/
#ifndef hyperReduction_H
#define hyperReduction_H
 
#include "fvCFD.H"
#include "ITHACAPOD.H"
#include "Foam2Eigen.H"
#include "EigenFunctions.H"
#include "ITHACAutilities.H"
#include "fvMeshSubset.H"
#include <set>
#include "redsvd"
 
 
 
// TODO uniform style for separating offline and onlie stages
template <typename... SnapshotsLists>
class HyperReduction
{
public:
    word methodName;
 
    using SnapshotsListTuple = std::tuple<std::decay_t<SnapshotsLists>...>;
    using FieldsTuple = std::tuple<typename std::decay_t<SnapshotsLists>::value_type ...>;
    static constexpr auto n_fields = sizeof...(SnapshotsLists);
 
    // TODO get vectorial dim from class template parameters
    unsigned int vectorial_dim;
    unsigned int sumFieldsDim;
 
    template <std::size_t N>
    using NthFieldListType = typename std::tuple_element<N, SnapshotsListTuple>::type;
 
    template <std::size_t N>
    using NthFieldType = typename NthFieldListType<N>::value_type;
 
    //----------------------------------------------------------------------
    HyperReduction(label n_modes,
                   label n_nodes,
                   Eigen::VectorXi initialSeeds,
                   word problemName,
                   SnapshotsLists &&...snapshotsLists);
    
    //----------------------------------------------------------------------
    HyperReduction(label n_modes,
                   label n_nodes,
                   unsigned int vectorial_dim, 
                   label n_cells,
                   Eigen::VectorXi initialSeeds,
                   word problemName);
    
    // ITHACA parameters dict
    ITHACAparameters* para;
 
    //TODO
    label n_boundary_patches;
 
    //TODO
    List<label> n_boundary_cells_list;
 
    //TODO
    label n_boundary_cells;
 
    label n_modes;
 
    label n_nodes;
 
    label n_snapshots;
 
    // Initial nodes
    Eigen::VectorXi initialSeeds;
 
    // nodes
    Eigen::VectorXi nodes;
 
    word problemName;
 
    SnapshotsListTuple snapshotsListTuple;
 
    List<word> fieldNames;
 
    List<unsigned int> fieldDims;
 
    label n_cells;
 
    label n_cellsSubfields;
 
    autoPtr<IOList<label>> nodePoints;
 
    autoPtr<IOList<labelList>> totalNodePoints;
 
    autoPtr<IOList<label>> uniqueNodePoints;
 
    word folderProblem;
 
    word folderMethod;
 
    List<label> localNodePoints;
    
    // Outputs
    Eigen::MatrixXd basisMatrix;
 
    // Normalizing weights of SVD decomposition
    Eigen::VectorXd normalizingWeights;
 
    Eigen::MatrixXd renormalizedBasisMatrix;
 
    // TODO
    Eigen::MatrixXd pinvPU;
 
    // TODO
    Eigen::MatrixXd wPU;
 
    // TODO
    Eigen::VectorXd eigenValueseig;
 
    Eigen::SparseMatrix<double> P;
 
    // Mask field to submesh
    Eigen::SparseMatrix<double> field2submesh;
 
    // Mask submesh to nodes
    Eigen::SparseMatrix<double> submesh2nodes;
 
    // Indices submesh to nodes
    Eigen::VectorXi submesh2nodesMask;
 
    Eigen::VectorXd quadratureWeights;
 
    autoPtr<fvMeshSubset> submesh;
    
    autoPtr<volVectorField> submesh_field;
 
    Eigen::MatrixXd MatrixOnline;
    
    //----------------------------------------------------------------------
    void offlineGappyDEIM(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& normalizingWeights){
        methodName = "GappyDEIM";
        word name = "ITHACAoutput/" + problemName + "/" + methodName;
        offlineGappyDEIM(snapshotsModes, normalizingWeights, name);
    }
 
    void offlineGappyDEIM(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& normalizingWeights, word folderMethodName);
 
    //----------------------------------------------------------------------
    void offlineECP(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& normalizingWeights){
        methodName = "ECP";
        word name = "ITHACAoutput/" + problemName + "/" + methodName;
        offlineECP(snapshotsModes, normalizingWeights, name);
    }
 
    void offlineECP(Eigen::MatrixXd& snapshotsModes, Eigen::VectorXd& normalizingWeights, word folderMethodName);
    
    //----------------------------------------------------------------------
    void initSeeds(Eigen::VectorXd mp_not_mask, std::set<label> nodePointsSet);
 
    //----------------------------------------------------------------------
    void computeLS(Eigen::MatrixXd& J, Eigen::MatrixXd& JWhole, Eigen::VectorXd& b, Eigen::VectorXd& q);
 
    //----------------------------------------------------------------------
    void getModesSVD(SnapshotsListTuple &SnapshotsListTuple, Eigen::MatrixXd &modesSVD, Eigen::VectorXd &fieldWeights, bool saveModesFlag = false);
 
    //----------------------------------------------------------------------
    void getModesSVD(SnapshotsListTuple& snapshotsListTuple, Eigen::MatrixXd& modesSVD, Eigen::VectorXd& fieldWeights, Eigen::MatrixXd& modesSVDBoundary, Eigen::VectorXd& fieldWeightsBoundary, bool saveModesFlag);
    
    //----------------------------------------------------------------------
    void updateNodes(Eigen::SparseMatrix<double>& P, label& ind_max, Eigen::VectorXd& mp_not_mask);
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    void stackSnapshots(SnapshotsList sList, Eigen::MatrixXd &snapshotsMatrix, Eigen::VectorXd &fieldWeights);
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    void stackSnapshotsBoundary(SnapshotsList sList, List<Eigen::MatrixXd>& snapshotsMatrixBoundary, List<Eigen::VectorXd>& fieldWeightsBoundary);
    
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    void saveModes(SnapshotsList sList, Eigen::MatrixXd& snapshotsMatrix, unsigned int &rowIndex, unsigned int &modeIndex, word folder);
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    void saveModes(SnapshotsList sList, Eigen::MatrixXd& snapshotsMatrix, Eigen::MatrixXd& snapshotsMatrixBoundary, unsigned int &rowIndex,  unsigned int& rowIndexBoundary, unsigned int &modeIndex, word folder);
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    void stackNames(SnapshotsList sList){
        fieldNames.append(sList[0].name());
    }
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    inline void stackDimensions(SnapshotsList sList){
        fieldDims.append(get_field_dim<typename SnapshotsList::value_type>());}
 
    //----------------------------------------------------------------------
    template <typename SnapshotsList>
    inline void sumDimensions(double sum, SnapshotsList sList){
        sum+=get_field_dim<typename SnapshotsList::value_type>();}
 
    //----------------------------------------------------------------------
    void evaluatePinv(Eigen::SparseMatrix<double> &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd& fieldWeights);
 
     //----------------------------------------------------------------------
    void evaluateWPU(Eigen::SparseMatrix<double> &Projector, Eigen::MatrixXd &Modes, Eigen::VectorXd& fieldWeights, Eigen::VectorXd& quadratureWeights);
 
    //----------------------------------------------------------------------
    void generateSubmesh(label layers, const fvMesh &mesh);
 
    //----------------------------------------------------------------------
    List<label> global2local(List<label> &points, fvMeshSubset &submesh);
 
    template <typename FieldType>
    inline autoPtr<FieldType> interpolateField(const FieldType &field)
    {
        return autoPtr<FieldType>(new FieldType(submesh->interpolate(field)));
    }
 
    // TODO expand with sphericalTensor, symmTensor, tensor
    template <typename Field>
    inline constexpr unsigned int get_field_dim()
    {
        return std::is_same<Field, volScalarField>::value ? 1 : 3;
    }
 
    template <typename LastList>
    inline constexpr unsigned int compute_vectorial_dim(LastList x) { return get_field_dim<typename std::decay_t<LastList>::value_type>(); }
 
    template <typename List, typename... RemainingLists>
    inline constexpr unsigned int compute_vectorial_dim(List &&head, RemainingLists &&...tail)
    {
        return get_field_dim<typename std::decay_t<List>::value_type>() + compute_vectorial_dim<RemainingLists...>(tail...);
    }
 
    static inline double s_optimality(Eigen::MatrixXd& A)
    {   
        return std::pow(
            std::sqrt((A.transpose()*A).determinant())/\
            A.colwise().lpNorm<2>().prod(),
            1.0/A.cols());
    }
 
    template <typename... FieldsArgs>
    void eigen2fields(Eigen::VectorXd& eFields, FieldsArgs&&... oFields)
    {   
        unsigned int cumulativeSize{0}, ith_field{0};
        
        ([&](auto& oFieldsArg)
            {
                double ith_fieldSize = fieldDims[ith_field]*oFieldsArg.size();
                Eigen::VectorXd vec = eFields.segment(cumulativeSize, ith_fieldSize);
                oFieldsArg = Foam2Eigen::Eigen2field(oFieldsArg, vec);
                cumulativeSize+=ith_fieldSize;
                ith_field++;
            }(oFields), ...);
 
    }
 
    void initReshapeMat(Eigen::SparseMatrix<double> &reshapeMat);
 
    void createMasks(bool offlineStage = true);
 
    void getSnapMatrix(Eigen::MatrixXd& snapMatrix, Eigen::VectorXd& fieldWeights);
 
    void getSnapMatrix(Eigen::MatrixXd& snapMatrix, Eigen::VectorXd& fieldWeights, List<Eigen::MatrixXd>& snapMatrixBoundary, List<Eigen::VectorXd>& fieldWeightsBoundary);
};
 
#endif