hyperReduction.H

/*---------------------------------------------------------------------------*\
     ██╗████████╗██╗  ██╗ █████╗  ██████╗ █████╗       ███████╗██╗   ██╗
     ██║╚══██╔══╝██║  ██║██╔══██╗██╔════╝██╔══██╗      ██╔════╝██║   ██║
     ██║   ██║   ███████║███████║██║     ███████║█████╗█████╗  ██║   ██║
     ██║   ██║   ██╔══██║██╔══██║██║     ██╔══██║╚════╝██╔══╝  ╚██╗ ██╔╝
     ██║   ██║   ██║  ██║██║  ██║╚██████╗██║  ██║      ██║      ╚████╔╝
     ╚═╝   ╚═╝   ╚═╝  ╚═╝╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝      ╚═╝       ╚═══╝
 
 * 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"
#include <unsupported/Eigen/NNLS>
 
 
 
// 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,
                       Eigen::VectorXd volumes = Eigen::VectorXd::Constant(1,0.0));
 
        // 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;

        Eigen::VectorXd volumes;

        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;
        List<Eigen::VectorXd> quadratureWeightsSubspaces;

        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)
        {
            List<Eigen::MatrixXd> listSnapshotsModes(1);
            listSnapshotsModes[0] = snapshotsModes.leftCols(n_modes);
            offlineECP(listSnapshotsModes, normalizingWeights, folderMethodName);
        }
 
        //----------------------------------------------------------------------
        void offlineECP(List<Eigen::MatrixXd>& listSnapshotsModes,
                        Eigen::VectorXd& weights, 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, bool NNLS);
 
        //----------------------------------------------------------------------
        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);
 
        //----------------------------------------------------------------------
        void updateNodesRemove(Eigen::SparseMatrix<double>& P, Eigen::VectorXd weights);
 
        //----------------------------------------------------------------------
        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);

        Eigen::SparseMatrix<double> maskToOtherDim(int newDim);
};
 
#endif