03steadyNS.C

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/>.
Description
    Example of steady NS Reduction Problem
SourceFiles
    03steadyNS.C
\*---------------------------------------------------------------------------*/
 
#include "IOmanip.H"
#include "ITHACAPOD.H"
#include "ITHACAstream.H"
#include "ReducedSteadyNS.H"
#include "forces.H"
#include "steadyNS.H"
 
class tutorial03 : public steadyNS
{
    public:
        explicit tutorial03(int argc, char* argv[])
            : steadyNS(argc, argv), U(_U()), p(_p()), args(_args()) {}

        volVectorField& U;
        volScalarField& p;
        argList& args;

        void offlineSolve()
        {
            Vector<double> inl(0, 0, 0);
            List<scalar> mu_now(1);
 
            // if the offline solution is already performed read the fields
            if (offline)
            {
                ITHACAstream::read_fields(Ufield, U, "./ITHACAoutput/Offline/");
                ITHACAstream::read_fields(Pfield, p, "./ITHACAoutput/Offline/");
                mu_samples =
                    ITHACAstream::readMatrix("./ITHACAoutput/Offline/mu_samples_mat.txt");
            }
            else
            {
                Vector<double> Uinl(0, 0, 0);
 
                for (label i = 0; i < mu.cols(); i++)
                {
                    mu_now[0] = mu(0, i);
                    change_viscosity(mu(0, i));
                    truthSolve(mu_now);
                    restart();
                }
            }
        }
};
 
void offline_stage(tutorial03& example);
void online_stage(tutorial03& example);
 
int main(int argc, char* argv[])
{
#if OPENFOAM > 1812
    // load stage to perform
    argList::addOption("stage", "offline", "Perform offline stage");
    argList::addOption("stage", "online", "Perform online stage");
    // add options for parallel run
    HashTable<string> validParOptions;
    validParOptions.set
    (
        "stage",
        "offline"
    );
    validParOptions.set
    (
        "stage",
        "online"
    );
    Pstream::addValidParOptions(validParOptions);
    // Construct the tutorial object
    tutorial03 example(argc, argv);
 
    if (example._args().get("stage").match("offline"))
    {
        // perform the offline stage, extracting the modes from the snapshots'
        // dataset corresponding to parOffline
        offline_stage(example);
    }
    else if (example._args().get("stage").match("online"))
    {
        // load precomputed modes and reduced matrices
        offline_stage(example);
        // perform online solve with respect to the parameters in parOnline
        online_stage(example);
    }
    else
    {
        Info << "Pass '-stage offline', '-stage online'" << endl;
    }
 
#else
 
    if (argc == 1)
    {
        Info << "Pass 'offline' or 'online' as first arguments."
             << endl;
        exit(0);
    }
 
    // process arguments removing "offline" or "online" keywords
    int argc_proc = argc - 1;
    char* argv_proc[argc_proc];
    argv_proc[0] = argv[0];
 
    if (argc > 2)
    {
        std::copy(argv + 2, argv + argc, argv_proc + 1);
    }
 
    argc--;
    // Construct the tutorial object
    tutorial03 example(argc, argv);
 
    if (std::strcmp(argv[1], "offline") == 0)
    {
        // perform the offline stage, extracting the modes from the snapshots' dataset corresponding to parOffline
        offline_stage(example);
    }
    else if (std::strcmp(argv[1], "online") == 0)
    {
        // load precomputed modes and reduced matrices
        offline_stage(example);
        // perform online solve with respect to the parameters in parOnline
        online_stage(example);
    }
    else
    {
        Info << "Pass offline, online" << endl;
    }
 
#endif
    exit(0);
}
 
void offline_stage(tutorial03& example)
{
    // Read some parameters from file
    ITHACAparameters* para =
        ITHACAparameters::getInstance(example._mesh(), example._runTime());
    int NmodesUout = para->ITHACAdict->lookupOrDefault<int>("NmodesUout", 15);
    int NmodesPout = para->ITHACAdict->lookupOrDefault<int>("NmodesPout", 15);
    int NmodesSUPout = para->ITHACAdict->lookupOrDefault<int>("NmodesSUPout", 15);
    int NmodesUproj = para->ITHACAdict->lookupOrDefault<int>("NmodesUproj", 10);
    int NmodesPproj = para->ITHACAdict->lookupOrDefault<int>("NmodesPproj", 10);
    int NmodesSUPproj =
        para->ITHACAdict->lookupOrDefault<int>("NmodesSUPproj", 10);
    // Read the par file where the training parameters are stored
    word filename("./parOffline");
    example.mu = ITHACAstream::readMatrix(filename);
    // Set the inlet boundaries patch 0 directions x and y
    example.inletIndex.resize(1, 2);
    example.inletIndex(0, 0) = 0;
    example.inletIndex(0, 1) = 0;
    // Perform the offline solve
    example.offlineSolve();
    // Solve the supremizer problem
    example.solvesupremizer();
 
    if (example.bcMethod == "lift")
    {
        ITHACAstream::read_fields(example.liftfield, example._U(), "./lift/");
        ITHACAutilities::normalizeFields(example.liftfield);
        // Homogenize the snapshots
        example.computeLift(example.Ufield, example.liftfield, example.Uomfield);
        // Perform POD on the velocity snapshots
        ITHACAPOD::getModes(example.Uomfield, example.Umodes, example._U().name(),
                            example.podex, 0, 0, NmodesUout);
    }
    else
    {
        ITHACAPOD::getModes(example.Ufield, example.Umodes, example._U().name(),
                            example.podex, 0, 0, NmodesUout);
    }
 
    // Perform POD on velocity pressure and supremizers and store the first 10
    // modes
    ITHACAPOD::getModes(example.Pfield, example.Pmodes, example._p().name(),
                        example.podex, 0, 0, NmodesPout);
    ITHACAPOD::getModes(example.supfield, example.supmodes, example._U().name(),
                        example.podex, example.supex, 1, NmodesSUPout);
    // Perform the Galerkin Projection
    example.projectSUP("./Matrices", NmodesUproj, NmodesPproj, NmodesSUPproj);
}
 
void online_stage(tutorial03& example)
{
    // Create the reduced object
    reducedSteadyNS reduced(example);
    // Read the par file where the test parameters are stored
    word filename("./parOnline");
    example.mu = ITHACAstream::readMatrix(filename);
    // Set the inlet velocity
    Eigen::MatrixXd vel_now(1, 1);
    vel_now(0, 0) = 1;
    // used only for penalty approach
    reduced.tauU = Eigen::MatrixXd::Zero(1, 1);
    reduced.tauU(0, 0) = 1e-1;
 
    // Perform an online solve for the new values of inlet velocities
    for (label k = 0; k < example.mu.size(); k++)
    {
        Info << "Evaluation of the reduced order model on the test set" << endl;
        Info << "Inlet Ux = " << vel_now(0, 0) << " nu = " << example.mu(0, k)
             << endl;
        // Set the reduced viscosity
        reduced.nu = example.mu(0, k);
        reduced.solveOnline_sup(vel_now);
        Eigen::MatrixXd tmp_sol(reduced.y.rows() + 1, 1);
        tmp_sol(0) = k + 1;
        tmp_sol.col(0).tail(reduced.y.rows()) = reduced.y;
        reduced.online_solution.append(tmp_sol);
    }
 
    // Save the online solution
    ITHACAstream::exportMatrix(reduced.online_solution, "red_coeff", "python",
                               "./ITHACAoutput/red_coeff");
    ITHACAstream::exportMatrix(reduced.online_solution, "red_coeff", "matlab",
                               "./ITHACAoutput/red_coeff");
    ITHACAstream::exportMatrix(reduced.online_solution, "red_coeff", "eigen",
                               "./ITHACAoutput/red_coeff");
    // Reconstruct and export the solution
    reduced.reconstruct(true, "./ITHACAoutput/Reconstruction/");
 
    if (Pstream::parRun())
    {
        bool endedOnline = true;
        reduce(endedOnline, sumOp<label>());
    }
}
 
//--------