UnsteadyNSTurb.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
 *     UnsteadyNSTurb
 *
 * Description
 *     Reduction class for a NON-Stationary turbulent NS problem
 *
 * SourceFiles
 *     UnsteadyNSTurb.C
\*---------------------------------------------------------------------------*/
 
#ifndef UnsteadyNSTurb_H
#define UnsteadyNSTurb_H
 
#include "IOMRFZoneList.H"
#include "IOporosityModelList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
#include "fvCFD.H"
#include "fvOptions.H"
#include "pimpleControl.H"
#include "singlePhaseTransportModel.H"
#include "steadyNS.H"
#include "turbulentTransportModel.H"
#include "unsteadyNS.H"
 
#include <bspline.h>
#include <bsplinebuilder.h>
#include <datatable.h>
#include <rbfspline.h>
#include <spline.h>
 
#include <iostream>
#include <vector>
 
#include "Smagorinsky.H"
#include "kOmegaSSTBase.H"
#include "DESModelBase.H"
#include "DESModel.H"
#include "kOmegaSSTDES.H"
#include "StoredParameters.H"
#include <optional>
 
/*---------------------------------------------------------------------------*\
                        Class UnsteadyNSTurb Declaration
\*---------------------------------------------------------------------------*/

class UnsteadyNSTurb : public unsteadyNS
{
    public:
        // Constructors
        UnsteadyNSTurb();
        UnsteadyNSTurb(int argc, char* argv[]);
 
        // Small construct with Parameters to access member functions linked to Smagorinsky diffusion 
        UnsteadyNSTurb(const Parameters* myParameters);
 
        // Specific variable for the unstationary case
        // Member Functions

        PtrList<volScalarField> nutFields;

        volScalarModes nutModes;

        PtrList<volScalarField> nutAve;
 
        // === Added for dual RBF interpolation (average and fluctuation νt) ===
        int nNutAvgModes;
        int nNutFluctModes;
        int nphiNutAvg;
        int nphiNutFluct;
 
        PtrList<volScalarField> nutAvgModes;
        PtrList<volScalarField> nutFluctModes;
        PtrList<volScalarField> nutAvgFields;
        PtrList<volScalarField> nutFluctFields;
 
        List<SPLINTER::RBFSpline*> rbfSplinesNutAvg;
        List<SPLINTER::RBFSpline*> rbfSplinesNutFluct;
 
        List<SPLINTER::DataTable*> samplesNutAvg;
        List<SPLINTER::DataTable*> samplesNutFluct;
 
        // (optional but recommended for compatibility)
        std::vector<SPLINTER::RBFSpline*> rbfSplines;

        std::vector<SPLINTER::DataTable*> samples;

        Eigen::MatrixXd btMatrix;
 
        // New: PPE RHS boundary vector L (size NPmodes x 1)
        Eigen::MatrixXd L_vector;
 
        // Builder prototype
        Eigen::MatrixXd pressurePPE_L(label NPmodes);

        Eigen::MatrixXd z;

        Eigen::Tensor<double, 3> ct1Tensor;

        Eigen::Tensor<double, 3> ct1AveTensor;

        Eigen::Tensor<double, 3> ct1PPETensor;

        Eigen::Tensor<double, 3> ct1PPEAveTensor;

        Eigen::Tensor<double, 3> ct2Tensor;

        Eigen::Tensor<double, 3> ct2AveTensor;

        Eigen::Tensor<double, 3> ct2PPETensor;

        Eigen::Tensor<double, 3> ct2PPEAveTensor;

        Eigen::Tensor<double, 3> ct1FluctTensor;

        Eigen::Tensor<double, 3> ct1PPEFluctTensor;

        Eigen::Tensor<double, 3> ct2FluctTensor;

        Eigen::Tensor<double, 3> ct2PPEFluctTensor;

        Eigen::Tensor<double, 3> cTotalTensor;

        Eigen::Tensor<double, 3> cTotalAveTensor;
 
        // ================= NEW FOR SUP: continuity matrix =================
        // P_{ij} = ( χ_i , ∇·φ_j ), size NP × (lift + NU + NSUP)
        Eigen::MatrixXd P_matrix;
        Eigen::MatrixXd continuity_matrix(label NUmodes, label NSUPmodes,
                                          label NPmodes);

        Eigen::Tensor<double, 3> cTotalPPETensor;

        Eigen::Tensor<double, 3> cTotalPPEAveTensor;
 
        Eigen::Tensor<double, 3> cTotalFluctTensor;
        Eigen::Tensor<double, 3> cTotalPPEFluctTensor;

        Eigen::MatrixXd bTotalMatrix;

        Eigen::MatrixXd coeffL2;

        Eigen::MatrixXd velRBF;

        Eigen::VectorXd radii;

        double e = 1.0;

        label nNutModes;

        label interChoice = 1;

        autoPtr<volScalarField> _nut;

        autoPtr<surfaceVectorField> _Uf;

        label _pRefCell;

        scalar _pRefValue;
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd btTurbulence(label NUmodes, label NSUPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceTensor1(label NUmodes, label NSUPmodes,
                label nNutModes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceTensor2(label NUmodes, label NSUPmodes,
                label nNutModes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPETensor1(
            label NUmodes, label NSUPmodes, label NPmodes, label nNutModes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPETensor2(
            label NUmodes, label NSUPmodes, label NPmodes, label nNutModes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceAveTensor1(label NUmodes, label NSUPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceAveTensor2(label NUmodes, label NSUPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPEAveTensor1(label NUmodes, label NSUPmodes,
                label NPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPEAveTensor2(label NUmodes, label NSUPmodes,
                label NPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceFluctTensor1(label NUmodes, label NSUPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulenceFluctTensor2(label NUmodes, label NSUPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPEFluctTensor1(label NUmodes,
                label NSUPmodes, label NPmodes);
 
        //--------------------------------------------------------------------------
        Eigen::Tensor<double, 3> turbulencePPEFluctTensor2(label NUmodes,
                label NSUPmodes, label NPmodes);
 
        //--------------------------------------------------------------------------
        void truthSolve(List<scalar> mu_now, std::string& offlinepath);
 
        //----------------------------------------------------------------------
        void projectSUP(
            fileName folder,
            label NUmodes,
            label NPmodes,
            label NSUPmodes,
            label nNutModes,
            bool rbfInterp = true);
 
        //----------------------------------------------------------------------
        void projectPPE(
            fileName folder,
            label NUmodes,
            label NPmodes,
            label NSUPmodes,
            label nNutModes,
            bool rbfInterp = true);
 
        //--------------------------------------------------------------------------
        List<Eigen::MatrixXd> velDerivativeCoeff(
            Eigen::MatrixXd A,
            Eigen::MatrixXd G,
            Eigen::VectorXd initSnapInd,
            Eigen::VectorXd timeSnap);
 
        //--------------------------------------------------------------------------
        List<Eigen::MatrixXd> velParCoeff(Eigen::MatrixXd A, Eigen::MatrixXd G);
 
        //--------------------------------------------------------------------------
        List<Eigen::MatrixXd> velParDerivativeCoeff(
            Eigen::MatrixXd A,
            Eigen::MatrixXd G,
            Eigen::VectorXd initSnapInd,
            Eigen::VectorXd timeSnap);
 
        //--------------------------------------------------------------------------
        Eigen::MatrixXd velParDerivativeCoeff(Eigen::MatrixXd A, Eigen::VectorXd par,
                                              double timeSnap);
 
private :

        const StoredParameters* m_parameters;
 
        // /// Turbulence model
        // autoPtr<incompressible::turbulenceModel> turbModel;
    
 public :
 
    // =========================================================================
    // ============================= Smagorinsky ===============================
    // =========================================================================
 
        //--------------------------------------------------------------------------
        void computeNonLinearSnapshot_at_time(const word& snap_time, volVectorField& Smag, 
                        std::optional<PtrList<volVectorField>> modesU = std::nullopt);
 
        //--------------------------------------------------------------------------
        volVectorField initSmagFunction();
 
        //--------------------------------------------------------------------------
        volScalarField initSmagPhiFunction(const volScalarField template_field_phi);
 
        //--------------------------------------------------------------------------
        volVectorField computeSmagTerm_at_time(const word& snap_time, std::optional<PtrList<volVectorField>> modesU = std::nullopt);
 
        //--------------------------------------------------------------------------
        volScalarField computeSmagTermPhi_at_time(const word& snap_time, const volScalarField template_field_phi);
 
        //--------------------------------------------------------------------------
        volVectorField computeSmagTerm_fromU(const volVectorField& snapshotj);
 
        //--------------------------------------------------------------------------
        volScalarField computeSmagTermPhi_fromUPhi(const volVectorField& snapshotj, const volScalarField& phij);
 
        //--------------------------------------------------------------------------
        volVectorField computeSmagTermPhi_fromUPhi(const volVectorField& snapshotj, const volVectorField& phij);
 
        //--------------------------------------------------------------------------
        template<typename T>
        static volScalarField diffusion(const T& coefDiff, const volScalarField& phi);
 
        //--------------------------------------------------------------------------
        template<typename T>
        static volVectorField diffusion(const T& coefDiff, const volVectorField& u);
 
 
    // =========================================================================
    // ======================= projFullStressFunction ==========================
    // =========================================================================
 
        //--------------------------------------------------------------------------
        void computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& phi, volVectorField& modeU);
 
        void computeNonLinearSnapshot_at_time(const word& snap_time, volVectorField& phi, volVectorField& modeU)
        {
            Info << "Output field with this call to computeNonLinearSnapshot_at_time is scalar" << endl;
            abort();
        };
 
        //--------------------------------------------------------------------------
        volScalarField initProjSmagFunction();
 
        //--------------------------------------------------------------------------
        volScalarField computeProjSmagTerm_at_time_fromMode(const word& snap_time, const volVectorField& mode);
 
        //--------------------------------------------------------------------------
        void computeProjSmagTerm_fromSmag_fromMode(volScalarField& phi, const volVectorField& Smag, const volVectorField& mode);
 
        void computeProjSmagTerm_fromSmag_fromMode(volScalarField phi,const volScalarField& Smag, const volVectorField& mode)
        {
            Info << "Input fields for this call to computeProjSmagTerm_fromSmag_fromMode are vectors" << endl;
            abort();
        };
 
        void computeProjSmagTerm_fromSmag_fromMode(volVectorField phi, const volVectorField& Smag, const volVectorField& mode)
        {
            Info << "Output field for this call to computeProjSmagTerm_fromSmag_fromMode is scalar" << endl;
            abort();
        };
 
 
    // =========================================================================
    // =========================== projSmagFromNut =============================
    // =========================================================================
 
        //--------------------------------------------------------------------------
        void computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& phi, 
                                volVectorField& modeU_proj , volVectorField& modeU_grad);
 
        void computeNonLinearSnapshot_at_time(const word& snap_time, volVectorField& phi, 
                                volVectorField& modeU_proj , volVectorField& modeU_grad)
        {
            Info << "Output field with this call to computeNonLinearSnapshot_at_time is scalar" << endl;
            abort();
        };
 
        //--------------------------------------------------------------------------
        volScalarField initProjSmagFromNutFunction();
 
        //--------------------------------------------------------------------------
        volScalarField computeProjSmagFromNut_at_time_fromModes(const word& snap_time, const volVectorField& modeU_proj, 
                                                           const volVectorField& modeU_grad);
 
        //--------------------------------------------------------------------------
        void computeProjSmagFromNut_fromNut_fromModes(volScalarField& phi, const volScalarField& Nut, 
                                const volVectorField& modeU_proj, const volVectorField& modeU_grad);
 
        void computeProjSmagFromNut_fromNut_fromModes(volVectorField& phi, const volVectorField& Nut, 
                                const volVectorField& modeU_proj, const volVectorField& modeU_grad)
        {
            Info << "Output field for this call to computeProjSmagFromNut_fromNut_fromModes is a scalar field" << endl;
            abort();
        };
 
        void computeProjSmagFromNut_fromNut_fromModes(volScalarField& phi, const volVectorField& Nut, 
                                const volVectorField& modeU_proj, const volVectorField& modeU_grad)
        {
            Info << "Input field for this call to computeProjSmagFromNut_fromNut_fromModes is a scalar field" << endl;
            abort();
        };
 
        //--------------------------------------------------------------------------
        static volScalarField projDiffusionIBP(const volScalarField& coefDiff, 
                             const volVectorField& u, const volVectorField& v);
 
 
 
    // =========================================================================
    // ============================= NUT =======================================
    // =========================================================================
 
        //--------------------------------------------------------------------------
        void computeNonLinearSnapshot_at_time(const word& snap_time, volScalarField& nut, 
              std::optional<PtrList<volVectorField>> modesU = std::nullopt);
 
        //--------------------------------------------------------------------------
        volScalarField initNutFunction();
 
        //--------------------------------------------------------------------------
        volScalarField computeNut_at_time(const word& snap_time, std::optional<PtrList<volVectorField>> modesU = std::nullopt);
 
        //--------------------------------------------------------------------------
        volScalarField computeNut_fromU(const volVectorField& snapshotj);
 
        //--------------------------------------------------------------------------
        volScalarField computeNut_fromS(const volTensorField& S);
 
        // void initTurbModel();
 
        //--------------------------------------------------------------------------
        static volTensorField computeS_fromU(const volVectorField& snapshotj);
 
        //--------------------------------------------------------------------------
        static volVectorField computeS_fromU(const volScalarField& phij);
 
};
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif