main.cpp

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <cctype>

#include "io/lookUpInFile.h"
#include "io/PBMWriter.h"

#include "numerics/TimeIntegrator.h"

#include "pbm/Pbm0d.h"
#include "pbm/PBMState.h"


static std::string to_key(std::string s)
{
    std::transform
    (
        s.begin(), s.end(), s.begin(),
        [](unsigned char c) { return char(std::tolower(c)); }
    );

    s.erase
    (
        std::remove_if
        (
            s.begin(), s.end(),
            [](unsigned char c) { return !std::isalnum(c); }
        ),
        s.end()
    );

    return s;
}


int main()
{
    const std::string inputFile  = "SetPbm";
    const std::string outputFile = "SetOutput";


    // -------------------------------------------------------------------------
    // Read PBM configuration
    // -------------------------------------------------------------------------

    PBMConfig cfg = setPBMConfig(inputFile);

    PBM0D pbm(cfg);


    // -------------------------------------------------------------------------
    // Read global simulation parameters
    // -------------------------------------------------------------------------

          double dt      = lookUpInFile<double>(inputFile, "dt");
    const double epsilon = lookUpInFile<double>(inputFile, "epsilon");
    const double rho_d   = lookUpInFile<double>(inputFile, "rho_d");
    const double tFinal  = lookUpInFile<double>(inputFile, "t_final");
    const double Ymin    = lookUpInFile<double>(inputFile, "Ymin");

    const double Y0 =
        lookUpInFile<double>(inputFile, "MassContent");


    // -------------------------------------------------------------------------
    // Construct initial PBM state
    // -------------------------------------------------------------------------

    std::vector<double> y(cfg.M + 1, 0.0);

    {
        const std::string raw =
            lookUpInFile<std::string>(inputFile, "Initial");

        const std::string key = to_key(raw);

        PBMState st0
        (
            std::vector<double>(cfg.M, 0.0),
            0.0,
            cfg.P,
            rho_d
        );

        if (key == "uniform")
        {
            const double Dmin =
                lookUpInFile<double>(inputFile, "Dmin");

            const double Dmax =
                lookUpInFile<double>(inputFile, "Dmax");

            st0 =
                PBMState::Uniform
                (
                    Dmin,
                    Dmax,
                    cfg.P,
                    Y0,
                    rho_d
                );
        }
        else if (key == "rosinrammler")
        {
            const double shape =
                lookUpInFile<double>(inputFile, "shape");

            const double scale =
                lookUpInFile<double>(inputFile, "scale");

            st0 =
                PBMState::RosinRammler
                (
                    shape,
                    scale,
                    cfg.P,
                    Y0,
                    rho_d
                );
        }
        else
        {
            std::cerr
                << "Unknown Initial = '" << raw << "'.\n"
                << "Supported types: Uniform, RosinRammler.\n";

            return 2;
        }

        for (int i = 0; i < cfg.M; ++i)
        {
            y[size_t(i)] = st0.N()[size_t(i)];
        }

        y[size_t(cfg.M)] = st0.Y();
    }


    // -------------------------------------------------------------------------
    // Create output writer
    // -------------------------------------------------------------------------

    PBMWriter writer(outputFile, cfg.M);


    // -------------------------------------------------------------------------
    // Time integration loop
    // -------------------------------------------------------------------------

    double time = 0.0;

    int nSteps = 0;

    while (time < tFinal)
    {
        // ---------------------------------------------------------------------
        // Write diagnostics from the current state
        // ---------------------------------------------------------------------

        {
            PBMState state(y, cfg.P, rho_d);

            writer.write(time, state);
        }

        ++nSteps;


        // ---------------------------------------------------------------------
        // Advance one time step
        // ---------------------------------------------------------------------

        const double dtUsed = dt;

        const bool accepted =
            Integrator::step(dt, epsilon, cfg.M + 1, pbm, rho_d, y);

        if (!accepted)
        {
            continue;
        }

        time += dtUsed;

        // ---------------------------------------------------------------------
        // Stop if the dispersed phase vanished
        // ---------------------------------------------------------------------

        if (y[size_t(cfg.M)] <= Ymin)
        {
            time += dtUsed;

            std::cout
                << "Dispersed phase vanished at step "
                << nSteps << ".\n";

            break;
        }
    }


    // -------------------------------------------------------------------------
    // Write final state
    // -------------------------------------------------------------------------

    {
        PBMState finalState(y, cfg.P, rho_d);

        writer.write(time, finalState);
    }


    // -------------------------------------------------------------------------
    // Final report
    // -------------------------------------------------------------------------

    std::cout
        << "Simulation completed in "
        << nSteps
        << " time steps.\n";


    return 0;
}