sequentialIHTP.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
    sequentialIHTP
Description
    A general full order implementation of an inverse problem
SourceFiles
    sequentialIHTP.C
\*---------------------------------------------------------------------------*/

#ifndef sequentialIHTP_H
#define sequentialIHTP_H
#include <iostream>
#include "fvCFD.H"
#include "interpolation.H"
#include "fvOptions.H"
#include "simpleControl.H"
#include "IOmanip.H"
#include "Time.H"
#include "laplacianProblem.H"
// #include "reducedLaplacian.H"
#include "ITHACAutilities.H"
#include <Eigen/Dense>
#include <cmath>
#include "Foam2Eigen.H"
#include "cnpy.H"
#include "mixedFvPatchFields.H"
#include "cellDistFuncs.H"
#include "ITHACAregularization.H"
 
#define _USE_MATH_DEFINES
 
using namespace SPLINTER;

class sequentialIHTP: public laplacianProblem
{
 
    public:
        // Constructors
        sequentialIHTP();
        sequentialIHTP(int argc, char* argv[]);
        //sequentialIHTP(int argc, char* argv[], scalar _diffusivity);
        virtual ~sequentialIHTP() {};
 
        ITHACAparameters* para;
 
        // Dummy variables to transform inverseLaplacianFoam into a class
        autoPtr<volScalarField> _T;
        PtrList<volScalarField> Ttime;
        List<PtrList<volScalarField>> Tbasis;

        autoPtr<volScalarField> _Tad;
        PtrList<volScalarField> Tad_time;
        PtrList<volScalarField> T0_time;

        mutable autoPtr<fvMesh> _mesh;

        autoPtr<simpleControl> _simple;

        autoPtr<fv::options> _fvOptions;

        autoPtr<Time> _runTime;
        autoPtr<Time> _runTimeShort;

        dimensionedScalar DT;

        scalar diffusivity = 0.0;

        PtrList<volScalarField> T0field;

        volScalarModes T0modes;

        Eigen::MatrixXd T0explicitMatrix_red;

        Eigen::MatrixXd T0implicitMatrix_red;

        Eigen::MatrixXd Tbasis_projectionMat;

        Eigen::VectorXd Tad_projected;

        List<label> magicPoints;

        Eigen::MatrixXd pointsReconstructMatrix;

        Eigen::MatrixXd pointTbasis_reconstructionMat;

        Eigen::VectorXd pointTad_reconstructed;

        PtrList<volScalarField> projectionErrorTbasis;

        //double projectionErrorTad;
        PtrList<volScalarField> projectionErrorTad;

        int NmodesT0 = 0;

        scalar startTime;
        scalar deltaTime;
        scalar endTime;
        label Ntimes;
        List<scalar> timeSteps;

        List<scalar> samplingTime;

        scalar timeSamplesDeltaT;

        scalar timeSamplesT0;

        label timeSamplesNum;

        label NtimeStepsBetweenSamples;

        label timeSampleI;

        label basisDeltaSample = 2;

        List<label> samplingSteps;

        bool thermocouplesRead = 0;

        int thermocouplesNum;

        int gBasisSize;

        double J;

        Eigen::VectorXd Jlist;

        double thermalCond;

        double HTC;

        double density;

        double specificHeat;

        word folderOffline = "./ITHACAoutput/offlineParamBC/";

        label offlineTimestepsSize = 0;

        scalar offlineEndTime = 0.0;
 
        bool offlineFlag = 0;
        bool interpolationFlag = 0;

        label NtimestepsInSequence = 0;
 
        scalar homogeneousBC = 0.0;
        List<scalar> homogeneousBCcoldSide;
        List<scalar> Tf; //temperature at coldSide [K]
        List<scalar> refGrad;
        List<scalar> valueFraction;

        label hotSide_ind;

        label coldSide_ind;

        List<List<scalar>> g;

        List<List<scalar>> heatFluxSpaceBasis;

        List<List<List<scalar>>> gBaseFunctions;

        label Nbasis = 0;

        label NbasisInTime = 0;

        label NbasisInSpace = 0;

        label NsamplesWindow = 0;

        List<scalar> gWeights;

        List<scalar> gWeightsOld;

        List<List<scalar>> gTrue;

        Eigen::VectorXd residual;
        Eigen::VectorXd addSol;
        Eigen::VectorXd T0_vector;
        Eigen::VectorXd Tcomp;
        Eigen::MatrixXd Theta;
        word timeBasisType = "None";
 
        word linSys_solver;
        label TSVD_filter;
        scalar Tikhonov_filter;
        label CG_Nsteps;
 
        List<vector> thermocouplesPos;
        List<int> thermocouplesCellID;
        List<int> thermocouplesCellProc;
        List<Foam::vector> thermocouplesCellC;
        Eigen::VectorXd Tmeas;   //Temperature at the thermocouples locations [K]
        Eigen::VectorXd
        TmeasShort; //Temperature at the thermocouples locations in only two sampling times [K]
        Eigen::VectorXd Tdirect;
        Eigen::VectorXd Tdiff;
 
        // Functions
 
        //--------------------------------------------------------------------------

        void setDiffusivity(scalar _diff);
 
        //--------------------------------------------------------------------------

        void setSpaceBasis(word type, scalar shapeParameter, label Npod = 0);
 
 
        //--------------------------------------------------------------------------

        void set_gParametrized(word spaceBaseFuncType,
                               scalar shapeParameter_space);
 
        //--------------------------------------------------------------------------

        volScalarField list2Field(List<scalar> list, scalar innerField = 0.0);
 
        //--------------------------------------------------------------------------

        List<List<scalar>> interpolateWeights(List<scalar> Wold, List<scalar> Wnew);
 
        //--------------------------------------------------------------------------

        void update_gParametrized(List<scalar> weights);
 
        //--------------------------------------------------------------------------

        void parameterizedBCoffline(bool force = 0);
 
        //--------------------------------------------------------------------------

        void reconstrucT(word outputFolder);
 
        //--------------------------------------------------------------------------

        Eigen::VectorXd reconstrucT(Eigen::VectorXi cells);
 
        //--------------------------------------------------------------------------

        void parameterizedBC(word folder, volScalarField initialField);
 
        //--------------------------------------------------------------------------

        void set_valueFraction();
 
        //--------------------------------------------------------------------------

        virtual void assignDirectBC(label timeI);
 
        //--------------------------------------------------------------------------

        virtual void assignT0_IF(volScalarField& T0_init) = 0;
 
        //--------------------------------------------------------------------------

        virtual void solveT0(volScalarField initialField);
 
        //--------------------------------------------------------------------------

        void getT0modes();
 
        //--------------------------------------------------------------------------

        void projectT0();
 
        //--------------------------------------------------------------------------

        void projectDirectOntoT0();
 
        //--------------------------------------------------------------------------

        void pointProjectionOffline();
 
        //--------------------------------------------------------------------------

        void projectionErrorOffline();
 
        //--------------------------------------------------------------------------

        void T0offline(int NmagicPoints);
 
 
        //--------------------------------------------------------------------------

        void solveAdditional();
 
        //--------------------------------------------------------------------------

        void solveDirect();
 
        //--------------------------------------------------------------------------

        virtual void readThermocouples();
 
        //--------------------------------------------------------------------------

        Eigen::VectorXd fieldValueAtThermocouples(volScalarField& field);
 
        //--------------------------------------------------------------------------

        Eigen::VectorXd fieldValueAtThermocouples(PtrList<volScalarField> fieldList,
                label fieldI);
 
        //--------------------------------------------------------------------------

        Eigen::VectorXd fieldValueAtThermocouples(PtrList<volScalarField> fieldList);
 
        //--------------------------------------------------------------------------

        void restart();
 
        //--------------------------------------------------------------------------

        void restartOffline();
 
        //--------------------------------------------------------------------------

        void restartT0();
 
        //--------------------------------------------------------------------------

        void sampling2symulationTime();
 
        //--------------------------------------------------------------------------

        void parameterizedBC_postProcess(List<Eigen::MatrixXd> linSys,
                                         Eigen::VectorXd weigths, word outputFolder, label verbose = 0);
 
        //--------------------------------------------------------------------------

        void findMagicPoints(int NmagicPoints);
 
};
 
#endif