Feat/crank–nicolson

.gitignore

@@ -25,4 +25,7 @@
 .temp/
 
 # Logs
-*.log
+*.log
+
+# Resources
+.temp

include/methods/crank_nicolson.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "method.hpp"
+#include <string>
+
+/**
+ * @brief Crank-Nicolson scheme implementation.
+ */
+class CrankNicolson : public Method
+{
+  public:
+    std::string name() const override
+    {
+        return "Crank-Nicolson";
+    }
+
+    bool uses_previous_step() const noexcept override
+    {
+        return false;
+    }
+
+    /**
+     * @brief Advance one time step using Crank-Nicolson scheme.
+     */
+    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/main.cpp

@@ -1,7 +1,7 @@
 /**
  * @file main.cpp
  * @brief Test driver for the 1D heat conduction solver.
- * Supports DuFort-Frankel, Richardson, and Laasonen schemes.
+ * Supports DuFort-Frankel, Richardson, Laasonen, and Crank-Nicolson schemes.
  */
 
 #include "grid.hpp"
@@ -29,7 +29,7 @@ namespace fs = std::filesystem;
  * @param phys       Physical parameters.
  * @param num        Numerical parameters.
  * @param scheme     The numerical scheme to use.
- * @param schemeName String name for output directory (e.g. "Laasonen").
+ * @param schemeName String name for output directory (e.g. "Crank-Nicolson").
  */
 void run_simulation(const PhysParams& phys,
                     const NumParams& num,
@@ -101,35 +101,42 @@ int main()
         std::cout << "\nSelect numerical method:\n";
         std::cout << "1. DuFort-Frankel\n";
         std::cout << "2. Richardson\n";
-        std::cout << "3. Laasonen\n";
-        std::cout << "4. Run All (Batch)\n";
+        std::cout << "3. Laasonen (Simple Implicit)\n";
+        std::cout << "4. Crank-Nicolson\n";
+        std::cout << "5. Run All (Batch)\n";
 
-        // We use askForDouble (or askForInteger) to get a valid choice between 1 and 4
+        // We use askForDouble (or askForInteger) to get a valid choice between 1 and 5
         int choice = static_cast<int>(askForDouble(
-                "Enter choice (1-4): ",
-                [](double v) { return v >= 1.0 && v <= 4.0; },
-                "Please enter a number between 1 and 4."));
+                "Enter choice (1-5): ",
+                [](double v) { return v >= 1.0 && v <= 5.0; },
+                "Please enter a number between 1 and 5."));
 
         // --- Execution Logic ---
 
         // 1. DuFort-Frankel
-        if (choice == 1 || choice == 4)
+        if (choice == 1 || choice == 5)
         {
             run_simulation(phys, num, SchemeKind::DuFortFrankel, "DuFort-Frankel");
         }
 
         // 2. Richardson
-        if (choice == 2 || choice == 4)
+        if (choice == 2 || choice == 5)
         {
             run_simulation(phys, num, SchemeKind::Richardson, "Richardson");
         }
 
         // 3. Laasonen
-        if (choice == 3 || choice == 4)
+        if (choice == 3 || choice == 5)
         {
             run_simulation(phys, num, SchemeKind::Laasonen, "Laasonen");
         }
 
+        // 4. Crank-Nicolson (New Addition!)
+        if (choice == 4 || choice == 5)
+        {
+            run_simulation(phys, num, SchemeKind::CrankNicolson, "Crank-Nicolson");
+        }
+
         std::cout << "\nAll requested simulations completed successfully ✅\n";
         return 0;
     }

src/method.cpp

@@ -5,31 +5,132 @@
 
 #include "method.hpp"
 
+// Include headers for ALL specific schemes
 #include "methods/dufort_frankel.hpp"
 #include "methods/laasonen.hpp"
 #include "methods/richardson.hpp"
+#include "methods/crank_nicolson.hpp"
+
 #include <memory>
 #include <stdexcept>
 
+// =============================================================================
+// ADAPTER CLASSES (Wrappers)
+// These classes adapt the specific scheme interfaces to the generic Method interface.
+// This allows the Solver to treat 2-level and 3-level schemes uniformly.
+// =============================================================================
+
+/**
+ * @brief Adapter for the Laasonen scheme (2-level, Implicit).
+ */
+class LaasonenMethod : public Method {
+public:
+    bool uses_previous_step() const override {
+        return false;
+    }
+
+    void step(const Grid& g, double D, double dt,
+              const std::vector<double>& /*Tprev*/, // Unused
+              const std::vector<double>& Tcurr,
+              std::vector<double>& Tnext) const override
+    {
+        // Delegate to concrete implementation
+        impl_.step(g, D, dt, Tnext, Tcurr);
+    }
+
+private:
+    Laasonen impl_;
+};
+
+/**
+ * @brief Adapter for the Richardson scheme (3-level, Explicit, Unstable).
+ */
+class RichardsonMethod : public Method {
+public:
+    bool uses_previous_step() const override {
+        return true;
+    }
+
+    void step(const Grid& g, double D, double dt,
+              const std::vector<double>& Tprev,
+              const std::vector<double>& Tcurr,
+              std::vector<double>& Tnext) const override
+    {
+        impl_.step(g, D, dt, Tnext, Tcurr, Tprev);
+    }
+
+private:
+    Richardson impl_;
+};
+
+/**
+ * @brief Adapter for the DuFort-Frankel scheme (3-level, Explicit, Stable).
+ */
+class DuFortFrankelMethod : public Method {
+public:
+    bool uses_previous_step() const override {
+        return true;
+    }
+
+    void step(const Grid& g, double D, double dt,
+              const std::vector<double>& Tprev,
+              const std::vector<double>& Tcurr,
+              std::vector<double>& Tnext) const override
+    {
+        impl_.step(g, D, dt, Tnext, Tcurr, Tprev);
+    }
+
+private:
+    DuFortFrankel impl_;
+};
+
+/**
+ * @brief Adapter for the Crank-Nicolson scheme (2-level, Implicit, Stable).
+ */
+class CrankNicolsonMethod : public Method {
+public:
+    bool uses_previous_step() const override {
+        return false;
+    }
+
+    void step(const Grid& g, double D, double dt,
+              const std::vector<double>& /*Tprev*/, // Unused
+              const std::vector<double>& Tcurr,
+              std::vector<double>& Tnext) const override
+    {
+        impl_.step(g, D, dt, Tnext, Tcurr);
+    }
+
+private:
+    CrankNicolson impl_;
+};
+
+// =============================================================================
+// FACTORY FUNCTION
+// =============================================================================
+
 /**
  * @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.
+ * @return std::unique_ptr<Method> Owning pointer to the created scheme wrapper.
  * @throws std::invalid_argument if the scheme is not supported.
  */
 std::unique_ptr<Method> make_method(SchemeKind scheme)
 {
     switch (scheme)
     {
         case SchemeKind::Laasonen:
-            return std::make_unique<Laasonen>();
+            return std::make_unique<LaasonenMethod>();
 
         case SchemeKind::DuFortFrankel:
-            return std::make_unique<DuFortFrankel>();
+            return std::make_unique<DuFortFrankelMethod>();
 
         case SchemeKind::Richardson:
-            return std::make_unique<Richardson>();
+            return std::make_unique<RichardsonMethod>();
+
+        case SchemeKind::CrankNicolson:
+            return std::make_unique<CrankNicolsonMethod>();
 
         default:
             throw std::invalid_argument("make_method: unsupported or unknown scheme");

src/methods/crank_nicolson.cpp

@@ -0,0 +1,85 @@
+/**
+ * @file crank_nicolson.cpp
+ * @brief Implementation of the Crank-Nicolson scheme using the Thomas Algorithm.
+ */
+
+#include "methods/crank_nicolson.hpp"
+
+#include <cstddef>
+#include <vector>
+
+void CrankNicolson::step(const Grid& g,
+                         double D,
+                         double dt,
+                         const std::vector<double>& /*Tprev*/,
+                         const std::vector<double>& Tcurr,
+                         std::vector<double>& Tnext)
+{
+    const std::size_t N = g.Nx;
+    const double dx = g.dx;
+
+    // Calculate the Fourier number r
+    const double r = D * dt / (dx * dx);
+    const double alpha = r / 2.0;
+
+    // Safety check for grid size
+    if (N < 3)
+        return;
+
+    // Resize Tnext if needed and copy current state (to preserve boundaries)
+    if (Tnext.size() != N)
+        Tnext.resize(N);
+    Tnext = Tcurr;
+
+    // Number of internal nodes (excluding boundaries 0 and N-1)
+    const std::size_t M = N - 2;
+
+    // --- 1. Set up the Tridiagonal Matrix Coefficients (LHS) ---
+    // The system is A * T^{n+1} = d
+    // A is constant: diag = (1+r), off-diag = -r/2
+    std::vector<double> a(M, -alpha);  // Lower diagonal
+    std::vector<double> b(M, 1.0 + r); // Main diagonal
+    std::vector<double> c(M, -alpha);  // Upper diagonal
+    std::vector<double> d(M, 0.0);     // Right-hand side vector
+
+    // --- 2. Construct the Right-Hand Side (RHS) ---
+    // d_i = (r/2)*T_{i-1}^n + (1-r)*T_i^n + (r/2)*T_{i+1}^n
+    // Note: j indexes the internal nodes (0 to M-1), corresponding to grid indices 1 to N-2.
+    for (std::size_t j = 0; j < M; ++j)
+    {
+        std::size_t i = j + 1; // Actual grid index
+        d[j] = alpha * Tcurr[i - 1] + (1.0 - r) * Tcurr[i] + alpha * Tcurr[i + 1];
+    }
+
+    // --- 3. Apply Boundary Conditions to RHS ---
+    // The terms -alpha * T^{n+1}_boundary are moved to the RHS.
+    // T[0] and T[N-1] already contain the fixed boundary values (Tsur).
+    d[0] += alpha * Tcurr[0];
+    d[M - 1] += alpha * Tcurr[N - 1];
+
+    // --- 4. Solve using Thomas Algorithm (TDMA) ---
+
+    // Forward elimination phase
+    for (std::size_t i = 1; i < M; ++i)
+    {
+        double m = a[i] / b[i - 1];
+        b[i] -= m * c[i - 1];
+        d[i] -= m * d[i - 1];
+    }
+
+    // Backward substitution phase
+    std::vector<double> x(M);
+    x[M - 1] = d[M - 1] / b[M - 1];
+
+    for (std::size_t i = M - 1; i-- > 0;)
+    {
+        x[i] = (d[i] - c[i] * x[i + 1]) / b[i];
+    }
+
+    // --- 5. Update Solution Vector ---
+    // Map the internal solution back to the full grid vector
+    for (std::size_t j = 0; j < M; ++j)
+    {
+        Tnext[j + 1] = x[j];
+    }
+}