main.cpp

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

#include <simple_partial/simple_partial.h>
#include <numsimulation/common.h>
using namespace numsimulation;
using namespace simple_partial;
#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;
    struct PoiseuilleParameters : AbstractParameters
    {
        double min[2], max[2];
        double p[2];
        double u[2];
        double density, viscosity;
        double a, b, c;
        virtual Eigen::Vector2d velocity(Eigen::Vector2d coord, double time) const
        {
            return Eigen::Vector2d(a + b * coord.y() + c * sqr(coord.y()), 0);
        }
        virtual double pressure(Eigen::Vector2d coord, double time) const
        {
            return p[0] + (p[1] - p[0]) * (max[0] - coord.x()) / (max[0] - min[0]);
        }
        PoiseuilleParameters(ParameterReader &reader) : AbstractParameters(reader),
            min{reader.get_real("XMIN"), reader.get_real("YMIN")},
            max{reader.get_real("XMAX"), reader.get_real("YMAX")},
            p{reader.get_real("PMIN"), reader.get_real("PMAX")},
            u{reader.get_real("VMIN"), reader.get_real("VMAX")},
            density(reader.get_real("DENSITY")),
            viscosity(reader.get_real("VISCOSITY"))
        {
            c = (p[1] - p[0]) / (2 * viscosity * (max[1] - min[0]));
            b = ((u[1] - u[0]) - c * (sqr(max[1]) - sqr(min[1]))) / (max[1] - min[1]);
            a = u[0] - b * min[1] - c * sqr(min[1]);
            boundaries.resize(4);
            boundaries[0].figure = new Line(Eigen::Vector2d(max[0], max[1]), Eigen::Vector2d(max[0], min[1]), false);  //Right
            boundaries[1].figure = new Line(Eigen::Vector2d(max[0], min[1]), Eigen::Vector2d(min[0], min[1]), false);  //Bottom
            boundaries[2].figure = new Line(Eigen::Vector2d(min[0], min[1]), Eigen::Vector2d(min[0], max[1]), false);  //Left
            boundaries[3].figure = new Line(Eigen::Vector2d(min[0], max[1]), Eigen::Vector2d(max[0], max[1]), false);  //Top
            grid_origin = Eigen::Vector2d((min[0] + max[0]) / 2, (min[1] + max[1]) / 2);
            exact = true;
        }
    };
    struct RotationParameters : AbstractParameters
    {
        double w;
        virtual Eigen::Vector2d velocity(Eigen::Vector2d coord, double time) const
        {
            return Eigen::Vector2d(-coord.y() * w, coord.x() * w);
        }
        virtual double pressure(Eigen::Vector2d coord, double time) const
        {
            return coord.squaredNorm() * sqr(w) * density / 2;
        }
        RotationParameters(ParameterReader &reader) : AbstractParameters(reader),
        w(reader.get_real("W"))
        {
            boundaries.resize(1);
            boundaries[0].figure = new Circle(Eigen::Vector2d::Zero(), reader.get_real("R"), true);
            grid_origin = Eigen::Vector2d::Zero();
            exact = true;
        }
    };
    if (reader.get_string("SCENARIO") == "POISEUILLE") parameters = new PoiseuilleParameters(reader);
    else if (reader.get_string("SCENARIO") == "ROTATION") parameters = new RotationParameters(reader);
    else throw std::runtime_error("_main(): Unrecognized scenario");
    
    //Read the rest of parameters
    unsigned int steps = reader.get_integer("STEPS");
    unsigned int divider = reader.get_integer("DIVIDER");
    reader.warn_unused();

    //Simulation
    Solver solver(parameters, argv[2]);
    solver.output();
    for (unsigned int step = 0, divider_step = 0; step < steps; step++)
    {
        solver.step();
        divider_step++;
        if (divider_step == divider)
        {
            solver.output();
            divider_step = 0;
        }
    }
    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;
    }
}