main.cpp

/* This file is a part of Numsimulation project
Developed as part of Numerical Flow Simulation course
Author: Kyrylo Sovailo */

#include <potential_warped/potential_warped.h>
#include <numsimulation/common.h>
using namespace numsimulation;
using namespace potential_warped;
#include <stdexcept>
#include <iostream>
#include <math.h>

int _main(int argc, char **argv)
{
    //Parsing arguments
    if (argc != 3) throw std::runtime_error("Invalid argument number");

    //Creating parameters object
    ParameterReader reader(argv[1]);
    AbstractParameters *parameters;

    ///Linear unwarped space, for testing purposes
    struct ParallelParameters : AbstractParameters
    {
        
        double vinf[2];
        
        virtual Eigen::Vector2d get_xy(Eigen::Vector2d xieta) const { return xieta; }
        virtual double get_dxi_dx(Eigen::Vector2d xieta) const { return 1.0; }
        virtual double get_dxi_dy(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dy(Eigen::Vector2d xieta) const { return 1.0; }
        virtual double get_dxi_dx2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_dxi_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double potential(Eigen::Vector2d xy) const { return vinf[0] * xy.x() + vinf[1] * xy.y(); }
        ParallelParameters(ParameterReader &reader) : AbstractParameters(reader),
            vinf{reader.get_real("VINFX"),
            reader.get_real("VINFY")}
            { exact = true; }
    };

    ///X is linear, Y narrows down with X
    struct StagnationPointParameters : AbstractParameters
    {
        double a;
        virtual Eigen::Vector2d get_xy(Eigen::Vector2d xieta) const { return Eigen::Vector2d(xieta(0), xieta(1) * min[0] / xieta(0)); }
        virtual double get_dxi_dx(Eigen::Vector2d xieta) const { return 1.0; }
        virtual double get_dxi_dy(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx(Eigen::Vector2d xieta) const { return get_xy(xieta).y() / min[0]; }
        virtual double get_deta_dy(Eigen::Vector2d xieta) const { return get_xy(xieta).x() / min[0]; }
        virtual double get_dxi_dx2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_dxi_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double potential(Eigen::Vector2d xy) const { return a * (sqr(xy.x()) - sqr(xy.y())); }
        StagnationPointParameters(ParameterReader &reader) : AbstractParameters(reader),
            a(reader.get_real("A"))
            { exact = true; }
    };

    ///Pipe-like space. X is linear, Y's top and bottom are sinewave-shaped
    struct PipeParameters : AbstractParameters
    {
        double vinf[2];
        double ymin, ymax;
        virtual Eigen::Vector2d get_xy(Eigen::Vector2d xieta) const
        {
            double local_ymin = min[1] + (ymin - min[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            double local_ymax = max[1] + (ymax - max[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            return Eigen::Vector2d(
                xieta(0),
                local_ymin + (local_ymax - local_ymin) * (xieta(1) - min[1]) / (max[1] - min[1])
            );
        }
        virtual double get_dxi_dx(Eigen::Vector2d xieta) const { return 1.0; }
        virtual double get_dxi_dy(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx(Eigen::Vector2d xieta) const
        {
            double local_ymin = min[1] + (ymin - min[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            double local_ymax = max[1] + (ymax - max[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            double dymin_dx = (ymin - min[1]) * 0.5 * sin(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])) / (max[0] - min[0]);
            double dymax_dx = (ymax - max[1]) * 0.5 * sin(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])) / (max[0] - min[0]);
            return (max[1] - min[1]) * (-dymin_dx * (local_ymax - local_ymin) - (get_xy(xieta).y() - local_ymin) * (dymax_dx - dymin_dx)) / sqr(local_ymax - local_ymin);
        }
        virtual double get_deta_dy(Eigen::Vector2d xieta) const
        {
            double local_ymin = min[1] + (ymin - min[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            double local_ymax = max[1] + (ymax - max[1]) * (0.5 - 0.5 * cos(M_PI * (xieta(0) - min[0]) / (max[0] - min[0])));
            return (max[1] - min[1]) / (local_ymax - local_ymin);
        }
        virtual double get_dxi_dx2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_dxi_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double get_deta_dx2(Eigen::Vector2d xieta) const { return 0.0; } // Not exactly zero, but I can't calculate it by hand
        virtual double get_deta_dy2(Eigen::Vector2d xieta) const { return 0.0; }
        virtual double potential(Eigen::Vector2d xy) const { return vinf[0] * xy.x() + vinf[1] * xy.y(); }
        PipeParameters(ParameterReader &reader) : AbstractParameters(reader),
            vinf{reader.get_real("VINFX"), reader.get_real("VINFY")},
            ymin(reader.get_real("YMIN")),
            ymax(reader.get_real("YMAX"))
            { exact = false; }
    };
    if (reader.get_string("SCENARIO") == "PARALLEL") parameters = new ParallelParameters(reader);
    else if (reader.get_string("SCENARIO") == "STAGNATION_POINT") parameters = new StagnationPointParameters(reader);
    else if (reader.get_string("SCENARIO") == "PIPE") parameters = new PipeParameters(reader);
    else throw std::runtime_error("_main(): Unrecognized scenario");
    
    //Unused parameter warning
    reader.warn_unused();

    //Setting simulation up
    Solver solver(parameters, argv[2]);
    solver.solve();
    solver.output();
    return 0;
}

int main(int argc, char **argv)
{
    try
    {
        return _main(argc, argv);
    }
    catch (const std::exception &e)
    {
        std::cerr << e.what() << "\n";
        return 1;
    }
    
}