Feat/richardson

.gitignore

@@ -15,7 +15,7 @@
 .cache/
 
 # Outputs
-/results/*
+/result/*
 !.gitkeep
 
 #MacOs Files

include/method.hpp

@@ -14,8 +14,9 @@
  */
 enum class SchemeKind
 {
-    DuFortFrankel
-    // Richardson, Laasonen, CrankNicolson (to be added later)
+    DuFortFrankel,
+    Richardson
+    // Laasonen, CrankNicolson (to be added later)
 };
 
 /**

include/methods/richardson.hpp

@@ -0,0 +1,60 @@
+#pragma once
+/**
+ * @file richardson.hpp
+ * @brief Richardson explicit (central time - central space) scheme for the 1D heat equation.
+ */
+
+#include <vector>
+#include <string>
+#include "../grid.hpp"
+#include "../method.hpp"
+
+/**
+ * @class Richardson
+ * @brief Implementation of the explicit Richardson scheme (CTCS).
+ *
+ * The scheme reads:
+ * \f[
+ *   T_i^{n+1} = T_i^{n-1} + 2 r \left( T_{i+1}^n - 2 T_i^n + T_{i-1}^n \right),
+ *   \quad r = \frac{D \, \Delta t}{\Delta x^2}
+ * \f]
+ *
+ * It is a three-level scheme, therefore it requires both \f$T^{n-1}\f$ and \f$T^n\f$
+ * to produce \f$T^{n+1}\f$. Boundary conditions are not enforced here, but left
+ * to the solver.
+ */
+class Richardson : public Method {
+public:
+    /**
+     * @brief Indicate that this is a three-level scheme.
+     * @return true always, since Richardson needs \f$T^{n-1}\f$.
+     */
+    bool uses_previous_step() const noexcept override { return true; }
+
+    /**
+     * @brief Get human-readable name of the scheme.
+     * @return "Richardson"
+     */
+    std::string name() const override { return "Richardson"; }
+
+    /**
+     * @brief Advance one time step using the Richardson explicit scheme.
+     *
+     * This computes \f$T^{n+1}\f$ from \f$T^{n}\f$ and \f$T^{n-1}\f$.
+     * Interior nodes \f$i = 1, \dots, N-2\f$ are updated, while boundary
+     * nodes are left untouched (the solver will impose Dirichlet BCs).
+     *
+     * @param g Grid metadata (spacing and number of nodes)
+     * @param D Diffusivity [cm^2/h]
+     * @param dt Time step [h]
+     * @param Tprev Temperature field at previous time layer \f$T^{n-1}\f$
+     * @param Tcurr Temperature field at current time layer \f$T^{n}\f$
+     * @param Tnext Output temperature field to be filled with \f$T^{n+1}\f$
+     */
+    void step(const Grid& g,
+              double D,
+              double dt,
+              const std::vector<double>& Tprev,
+              const std::vector<double>& Tcurr,
+              std::vector<double>& Tnext) override;
+};

src/io.cpp

@@ -1,8 +1,9 @@
 #include "io.hpp"
-#include <fstream>
+
 #include <filesystem>
-#include <stdexcept>
+#include <fstream>
 #include <iomanip>
+#include <stdexcept>
 
 namespace fs = std::filesystem;
 
@@ -23,9 +24,12 @@ void save_profile_csv(const std::string& path,
     if (x.size() != T.size())
         throw std::runtime_error("save_profile_csv: vector size mismatch");
 
-    // --- Ensure the result directory exists ---
-    fs::path fullPath = fs::path("result") / path;
-    fs::create_directories(fullPath.parent_path());
+    // --- Ensure the directory exists ---
+    fs::path fullPath(path);
+    if (fullPath.has_parent_path())
+    {
+        fs::create_directories(fullPath.parent_path());
+    }
 
     // --- Open file for writing ---
     std::ofstream file(fullPath);

src/main.cpp

@@ -1,67 +1,120 @@
 /**
  * @file main.cpp
- * @brief Run DuFort-Frankel simulation with user-chosen dx, dt; dump CSV every 0.1 h.
+ * @brief Test driver for the 1D heat conduction solver.
+ *        Supports DuFort-Frankel and Richardson schemes.
  */
+
 #include "grid.hpp"
 #include "io.hpp"
+#include "method.hpp"
 #include "solver.hpp"
 #include "types.hpp"
 #include "user_input.h"
 #include <filesystem>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
+#include <string>
+
+namespace fs = std::filesystem;
+
+/**
+ * @brief Run a single simulation with the specified scheme.
+ *
+ * @param phys Physical parameters
+ * @param num Numerical parameters
+ * @param scheme The numerical scheme to use
+ * @param schemeName String name for output directory
+ */
+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 ---
+    fs::path outDir = fs::path("results") / schemeName;
+    if (!fs::exists(outDir))
+    {
+        fs::create_directories(outDir);
+    }
+
+    // --- Output callback ---
+    auto onDump = [&](double t, const Grid& g, const std::vector<double>& T)
+    {
+        std::ostringstream name;
+        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);
+        // Optional: reduce verbosity for batch runs
+        // std::cout << "Saved " << fullPath.string() << "\n";
+    };
+
+    // --- Run ---
+    solver.run(onDump);
+    std::cout << "Results saved to " << outDir << "\n";
+}
 
 int main()
 {
     try
     {
-        PhysParams phys; // L=31 cm, D=93 cm^2/h, Tin=38°C, Tsur=149°C
-        NumParams num;   // defaults: dx=0.05, dt=0.01, tEnd=0.5, outEvery=0.1
+        std::cout << "=== 1D Heat Equation Solver ===\n";
 
-        // Ask user for dx, dt (dt choices per assignment)
-        std::cout << "=== DuFort-Frankel simulation ===\n";
+        // --- Physical parameters ---
+        PhysParams phys;
+        phys.L_cm = 31.0;
+        phys.Tin = 38.0;
+        phys.Tsur = 149.0;
+        phys.D_cm2h = 93.0;
+
+        // --- 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 in (0, 1].");
-        // --- Run DuFort-Frankel Solver ---
-        std::cout << "=== Running DuFort-Frankel Solver ===\n";
-        // Use default dt and tEnd from NumParams or ask user?
-        // For now use defaults: dt=0.01, tEnd=0.5
-
-        HeatSolver solver(phys, num, SchemeKind::DuFortFrankel);
-
-        // --- Prepare output directory ---
-        namespace fs = std::filesystem;
-        std::string schemeName =
-                "DuFort-Frankel"; // Could get this from solver.method_name() if exposed, but
-                                  // hardcoding for now is fine or better yet, expose it.
-        // Actually, let's just hardcode "DuFort-Frankel" as requested or use the scheme enum.
-        // The user said: results/"method name"/...
-
-        fs::path outDir = fs::path("results") / schemeName;
-        if (!fs::exists(outDir))
+                "Δx must be positive and <= 1.0 cm.");
+        num.dt = 0.01;
+        num.tEnd = 0.5;
+        num.outEvery = 0.1;
+
+        // --- Grid validation ---
+        Grid g(phys.L_cm, num.dx);
+        validate_grid(g);
+        std::cout << "Grid created ✅  Nx = " << g.Nx << ", dx = " << g.dx << " cm\n";
+
+        // --- Method Selection ---
+        std::cout << "\nSelect numerical method:\n";
+        std::cout << "1. DuFort-Frankel\n";
+        std::cout << "2. Richardson\n";
+        std::cout << "3. Run All (Batch)\n";
+
+        int choice = static_cast<int>(askForDouble(
+                "Enter choice (1-3): ",
+                [](double v) { return v >= 1.0 && v <= 3.0; },
+                "Please enter 1, 2, or 3."));
+
+        if (choice == 1 || choice == 3)
         {
-            fs::create_directories(outDir);
+            run_simulation(phys, num, SchemeKind::DuFortFrankel, "DuFort-Frankel");
         }
 
-        auto onDump = [&](double t, const Grid& g, const std::vector<double>& T)
+        if (choice == 2 || choice == 3)
         {
-            std::ostringstream name;
-            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);
-            std::cout << "Saved " << fullPath.string() << "\n";
-        };
-
-        solver.run(onDump);
-        std::cout << "Simulation complete.\n";
+            run_simulation(phys, num, SchemeKind::Richardson, "Richardson");
+        }
+
+        std::cout << "\nAll requested simulations completed successfully ✅\n";
         return 0;
     }
     catch (const std::exception& e)
     {
-        std::cerr << "Error: " << e.what() << "\n";
+        std::cerr << "Fatal error: " << e.what() << "\n";
         return 1;
     }
 }

src/method.cpp

@@ -1,16 +1,25 @@
 #include "method.hpp"
 
 #include "methods/dufort_frankel.hpp"
+#include "methods/richardson.hpp"
 #include <memory>
 #include <stdexcept>
 
+/**
+ * @brief Build a concrete numerical method from the requested scheme kind.
+ *
+ * @param scheme Identifier of the requested scheme.
+ * @return std::unique_ptr<Method> Owning pointer to the created scheme.
+ * @throws std::invalid_argument if the scheme is not available in this build.
+ */
 std::unique_ptr<Method> make_method(SchemeKind scheme)
 {
     switch (scheme)
     {
-        // @brief Factory implementation.
         case SchemeKind::DuFortFrankel:
             return std::make_unique<DuFortFrankel>();
+        case SchemeKind::Richardson:
+            return std::make_unique<Richardson>();
         default:
             throw std::invalid_argument("make_method: unsupported scheme");
     }

src/methods/richardson.cpp

@@ -0,0 +1,25 @@
+/**
+ * @file richardson.cpp
+ * @brief Richardson explicit scheme implementation.
+ */
+
+#include "methods/richardson.hpp"
+#include <cstddef>  // for std::size_t
+
+void Richardson::step(const Grid& g,
+                      double D,
+                      double dt,
+                      const std::vector<double>& Tprev,
+                      const std::vector<double>& Tcurr,
+                      std::vector<double>& Tnext)
+{
+    const double r        = D * dt / (g.dx * g.dx);
+    const std::size_t N   = g.Nx;
+
+    // Update interior nodes only; boundaries are imposed by the solver.
+    for (std::size_t i = 1; i + 1 < N; ++i) {
+        Tnext[i] = Tprev[i]
+                 + 2.0 * r * ( Tcurr[i + 1] - 2.0 * Tcurr[i] + Tcurr[i - 1] );
+    }
+    // Tnext[0] and Tnext[N-1] left as-is
+}

src/solver.cpp

@@ -1,11 +1,16 @@
+/**
+ * @file solver.cpp
+ * @brief Time-marching driver for the 1D heat equation.
+ */
+
 #include "solver.hpp"
 
-#include <algorithm>
+#include <algorithm> // std::fill
 
 HeatSolver::HeatSolver(const PhysParams& phys, const NumParams& num, SchemeKind scheme)
     : phys_(phys)
     , num_(num)
-    , grid_(phys.L_cm, num.dx) // grid construction
+    , grid_(phys.L_cm, num.dx)
 {
     // Create the calculation method (factory based on scheme)
     method_ = make_method(scheme);
@@ -35,10 +40,7 @@ void HeatSolver::apply_dirichlet(std::vector<double>& T) const
 
 void HeatSolver::init_field()
 {
-    // Uniform initial field at Tin
     std::fill(T_.begin(), T_.end(), phys_.Tin);
-
-    // Boundary conditions
     apply_dirichlet(T_);
 }
 
@@ -59,7 +61,6 @@ void HeatSolver::bootstrap_if_needed()
     const double r = phys_.D_cm2h * num_.dt / (dx * dx);
     const std::size_t N = grid_.Nx;
 
-    // Start from T^0
     next_ = T_;
 
     // FTCS scheme (conditional stability: r <= 0.5)
@@ -76,11 +77,8 @@ void HeatSolver::bootstrap_if_needed()
         next_ = T_;
     }
 
-    // Boundary conditions
     apply_dirichlet(next_);
-
-    // 3) Update T_ to contain T^1
-    T_ = next_;
+    T_ = next_; // T^1
 }
 
 void HeatSolver::run(
@@ -90,18 +88,17 @@ void HeatSolver::run(
     const double dt = num_.dt;
     const double tEnd = num_.tEnd;
     const double outEvery = num_.outEvery;
-
     double nextDump = 0.0;
 
-    // Initial dump
+    // Initial output
     onDump(t, grid_, T_);
 
     // Main time loop
     while (t < tEnd - 1e-12)
     {
         if (uses_prev_step_)
         {
-            // Three-level method (e.g. DuFort-Frankel)
+            // Three-level method (e.g. DuFort-Frankel, Richardson)
             method_->step(grid_, phys_.D_cm2h, dt, Tprev_, T_, next_);
         }
         else