diff --git a/.gitignore b/.gitignore index 154bef8..3bdb058 100644 --- a/.gitignore +++ b/.gitignore @@ -25,4 +25,7 @@ .temp/ # Logs -*.log \ No newline at end of file +*.log + +# Resources +.temp \ No newline at end of file diff --git a/include/methods/crank_nicolson.hpp b/include/methods/crank_nicolson.hpp new file mode 100644 index 0000000..5129e90 --- /dev/null +++ b/include/methods/crank_nicolson.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include "method.hpp" +#include + +/** + * @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& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) override; +}; \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 6318292..b07e2ed 100644 --- a/src/main.cpp +++ b/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(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; } diff --git a/src/method.cpp b/src/method.cpp index 867dbac..3bb195b 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -5,17 +5,115 @@ #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 #include +// ============================================================================= +// 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& /*Tprev*/, // Unused + const std::vector& Tcurr, + std::vector& 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& Tprev, + const std::vector& Tcurr, + std::vector& 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& Tprev, + const std::vector& Tcurr, + std::vector& 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& /*Tprev*/, // Unused + const std::vector& Tcurr, + std::vector& 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 Owning pointer to the created scheme. + * @return std::unique_ptr Owning pointer to the created scheme wrapper. * @throws std::invalid_argument if the scheme is not supported. */ std::unique_ptr make_method(SchemeKind scheme) @@ -23,13 +121,16 @@ std::unique_ptr make_method(SchemeKind scheme) switch (scheme) { case SchemeKind::Laasonen: - return std::make_unique(); + return std::make_unique(); case SchemeKind::DuFortFrankel: - return std::make_unique(); + return std::make_unique(); case SchemeKind::Richardson: - return std::make_unique(); + return std::make_unique(); + + case SchemeKind::CrankNicolson: + return std::make_unique(); default: throw std::invalid_argument("make_method: unsupported or unknown scheme"); diff --git a/src/methods/crank_nicolson.cpp b/src/methods/crank_nicolson.cpp new file mode 100644 index 0000000..7c0c72a --- /dev/null +++ b/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 +#include + +void CrankNicolson::step(const Grid& g, + double D, + double dt, + const std::vector& /*Tprev*/, + const std::vector& Tcurr, + std::vector& 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 a(M, -alpha); // Lower diagonal + std::vector b(M, 1.0 + r); // Main diagonal + std::vector c(M, -alpha); // Upper diagonal + std::vector 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 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]; + } +} \ No newline at end of file