main.cpp

Go to the documentation of this file.

/**
 * @file main.cpp
 * @brief Test driver for the 1D heat conduction solver.
 * Supports DuFort-Frankel, Richardson, Laasonen, and Crank-Nicolson schemes.
 */
 
#include "grid.hpp"
#include "io.hpp"
#include "method.hpp"
#include "methods/analytical.hpp"
#include "solver.hpp"
#include "types.hpp"
#include "user_input.h"
#include <filesystem>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace fs = std::filesystem;
 
/**
 * @brief Run a single simulation with the specified scheme.
 *
 * This helper function avoids code duplication by handling the solver setup,
 * directory creation, and output callback for any given scheme.
 *
 * @param phys       Physical parameters.
 * @param num        Numerical parameters.
 * @param scheme     The numerical scheme to use.
 * @param schemeName String name for output directory (e.g. "Crank-Nicolson").
 */
void run_simulation(const PhysParams& phys,
                    const NumParams& num,
                    SchemeKind scheme,
                    const std::string& schemeName)
{
    std::cout << "\n=== Running " << schemeName << " Solver ===\n";
 
    // --- Create solver ---
    HeatSolver solver(phys, num, scheme);
 
    // --- Prepare output directory ---
    // Results will be saved in: results/<schemeName>/
    fs::path outDir = fs::path("results") / schemeName;
    if (!fs::exists(outDir))
    {
        fs::create_directories(outDir);
    }
 
    // --- Output callback ---
    // Called by the solver at t=0 and every 'outEvery' hours.
    auto onDump = [&](double t, const Grid& g, const std::vector<double>& T)
    {
        std::ostringstream name;
        // Filename format: profile_t_0.10h.csv
        name << "profile_t_" << std::fixed << std::setprecision(2) << t << "h.csv";
 
        fs::path fullPath = outDir / name.str();
        save_profile_csv(fullPath.string(), g.x, T);
    };
 
    // --- Run ---
    solver.run(onDump);
    std::cout << "Results saved to " << outDir << "\n";
}
 
/**
 * @brief Main entry point for the 1D Heat Equation Solver.
 *
 * This function orchestrates the entire simulation workflow:
 * 1. Defines physical and numerical parameters.
 * 2. Validates the grid.
 * 3. Prompts the user for method selection.
 * 4. Executes the selected simulation(s).
 *
 * @return int 0 on success, 1 on error.
 */
int main()
{
    try
    {
        std::cout << "=== 1D Heat Equation Solver ===\n";
 
        // ---------------------------------------------------------------------
        // 1. Define Physical Parameters
        // ---------------------------------------------------------------------
        PhysParams phys;
        phys.L_cm = 31.0;
        phys.Tin = 38.0;
        phys.Tsur = 149.0;
        phys.D_cm2h = 93.0;
 
        // ---------------------------------------------------------------------
        // 2. Define Numerical Parameters
        // ---------------------------------------------------------------------
        NumParams num;
 
        // Ask for dx once
        num.dx = askForDouble(
                "Enter Δx [cm] (suggested 0.05): ",
                [](double v) { return v > 0.0 && v <= 1.0; },
                "Δx must be positive and <= 1.0 cm.");
 
        num.dt = 0.01;      // Time step [h]
        num.tEnd = 0.5;     // Total duration [h]
        num.outEvery = 0.1; // Output interval [h]
 
        // ---------------------------------------------------------------------
        // 3. Grid Validation (Early Check)
        // ---------------------------------------------------------------------
        Grid g(phys.L_cm, num.dx);
        validate_grid(g);
        std::cout << "Grid created ✅  Nx = " << g.Nx << ", dx = " << g.dx << " cm\n";
 
        // ---------------------------------------------------------------------
        // 4. Method Selection
        // ---------------------------------------------------------------------
        std::cout << "\nSelect numerical method:\n";
        std::cout << "1. DuFort-Frankel\n";
        std::cout << "2. Richardson\n";
        std::cout << "3. Laasonen (Simple Implicit)\n";
        std::cout << "4. Crank-Nicolson\n";
        std::cout << "5. Run All (Batch + Analytical)\n";
        std::cout << "6. Laasonen Stability Study (dt = 0.01, 0.025, 0.05, 0.1)\n";
 
        // We use askForDouble (or askForInteger) to get a valid choice between 1 and 6
        int choice = static_cast<int>(askForDouble(
                "Enter choice (1-6): ",
                [](double v) { return v >= 1.0 && v <= 6.0; },
                "Please enter a number between 1 and 6."));
 
        // ---------------------------------------------------------------------
        // 5. Execution Logic
        // ---------------------------------------------------------------------
 
        // Case 1: DuFort-Frankel
        if (choice == 1 || choice == 5)
        {
            run_simulation(phys, num, SchemeKind::DuFortFrankel, "DuFort-Frankel");
        }
 
        // Case 2: Richardson
        if (choice == 2 || choice == 5)
        {
            run_simulation(phys, num, SchemeKind::Richardson, "Richardson");
        }
 
        // Case 3: Laasonen
        if (choice == 3 || choice == 5)
        {
            run_simulation(phys, num, SchemeKind::Laasonen, "Laasonen");
        }
 
        // Case 4: Crank-Nicolson
        if (choice == 4 || choice == 5)
        {
            run_simulation(phys, num, SchemeKind::CrankNicolson, "Crank-Nicolson");
        }
 
        // Case 5: Analytical Solution (Only in Batch or if we added a specific option)
        if (choice == 5)
        {
            std::cout << "\n=== Computing Analytical Solution ===\n";
            std::string outDir = "results/Analytical";
            fs::create_directories(outDir);
 
            // Generate profiles at same intervals as numerical methods
            // t = 0, outEvery, 2*outEvery, ... tEnd
            int steps = static_cast<int>(num.tEnd / num.outEvery);
 
            // Explicitly export t=0.01 for comparison
            {
                double t = 0.01;
                std::vector<double> T_analytical = Analytical::get_profile(t, g.x, phys);
                std::ostringstream name;
                name << outDir << "/profile_t_" << std::fixed << std::setprecision(2) << t
                     << "h.csv";
                save_profile_csv(name.str(), g.x, T_analytical);
            }
 
            for (int i = 0; i <= steps; ++i)
            {
                double t = i * num.outEvery;
                // Avoid tiny floating point errors
                if (t > num.tEnd)
                    t = num.tEnd;
 
                std::vector<double> T_analytical = Analytical::get_profile(t, g.x, phys);
 
                std::ostringstream name;
                name << outDir << "/profile_t_" << std::fixed << std::setprecision(2) << t
                     << "h.csv";
                save_profile_csv(name.str(), g.x, T_analytical);
            }
            std::cout << "Results saved to " << outDir << "\n";
        }
 
        // Case 6: Laasonen Stability Study
        if (choice == 6)
        {
            std::cout << "\n=== Running Laasonen Stability Study ===\n";
            std::vector<double> dts = {0.01, 0.025, 0.05, 0.1};
 
            for (double dt : dts)
            {
                // Create a local copy of params to modify dt
                NumParams localNum = num;
                localNum.dt = dt;
 
                std::ostringstream name;
                name << "Laasonen_dt_" << std::fixed << std::setprecision(3) << dt;
                // Remove trailing zeros/dot if desired, or keep fixed precision for sorting
                // "Laasonen_dt_0.010", etc.
 
                run_simulation(phys, localNum, SchemeKind::Laasonen, name.str());
            }
        }
 
        std::cout << "\nAll requested simulations completed successfully ✅\n";
        return 0;
    }
    catch (const std::exception& e)
    {
        std::cerr << "Fatal error: " << e.what() << "\n";
        return 1;
    }
}