From 19bbf407cbccc19c2c77d3b308d09b75d58d39bb Mon Sep 17 00:00:00 2001 From: Broky64 Date: Mon, 1 Dec 2025 08:58:38 +0100 Subject: [PATCH 1/2] feat: Add Crank-Nicolson numerical scheme implementation and integrate it into the heat solver. --- .gitignore | 6 ++- include/method.hpp | 24 ++++++--- include/methods/crank_nicolson.hpp | 25 +++++++++ src/main.cpp | 45 +++++++++------- src/method.cpp | 87 ++++++++++++++++++++++++++---- src/methods/crank_nicolson.cpp | 73 +++++++++++++++++++++++++ 6 files changed, 223 insertions(+), 37 deletions(-) create mode 100644 include/methods/crank_nicolson.hpp create mode 100644 src/methods/crank_nicolson.cpp diff --git a/.gitignore b/.gitignore index 2787c44..0dcb83b 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ # IDE .vscode/ .idea/ +.cache/ # Outputs /results/* @@ -21,4 +22,7 @@ .DS_Store # Logs -*.log \ No newline at end of file +*.log + +# Resources +.temp \ No newline at end of file diff --git a/include/method.hpp b/include/method.hpp index 5488170..1ed16c3 100644 --- a/include/method.hpp +++ b/include/method.hpp @@ -5,23 +5,26 @@ */ #include "grid.hpp" -#include #include #include +#include /** * @brief Supported numerical schemes. */ -enum class SchemeKind { - DuFortFrankel - // Richardson, Laasonen, CrankNicolson (to be added later) +enum class SchemeKind +{ + DuFortFrankel, + CrankNicolson + // Richardson, Laasonen (to be added later) }; /** * @brief Abstract base class for time-stepping methods. */ -class Method { -public: +class Method +{ + public: virtual ~Method() = default; /// @return human-readable scheme name. @@ -37,7 +40,9 @@ class Method { * @param Tcurr Read-only temperature at the current time layer T^{n} * @param Tnext Output container that must be filled with T^{n+1} */ - virtual void step(const Grid& g, double D, double dt, + virtual void step(const Grid& g, + double D, + double dt, const std::vector& Tprev, const std::vector& Tcurr, std::vector& Tnext) = 0; @@ -48,7 +53,10 @@ class Method { * Three-level schemes (e.g. DuFort-Frankel) require both T^{n} and T^{n-1} * to produce T^{n+1}. Two-level schemes can ignore the returned value. */ - virtual bool uses_previous_step() const noexcept { return false; } + virtual bool uses_previous_step() const noexcept + { + return false; + } }; /** diff --git a/include/methods/crank_nicolson.hpp b/include/methods/crank_nicolson.hpp new file mode 100644 index 0000000..97892e5 --- /dev/null +++ b/include/methods/crank_nicolson.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include "grid.hpp" +#include + +/** + * @brief Crank-Nicolson scheme implementation. + */ +class CrankNicolson +{ + public: + /** + * @brief Advance one time step using Crank-Nicolson scheme. + * @param g Grid metadata + * @param D Diffusivity + * @param dt Time step + * @param T Output container for T^{n+1} + * @param Tn Input container for T^{n} + */ + void step(const Grid& g, + double D, + double dt, + std::vector& T, + const std::vector& Tn) const; +}; \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 08d4330..5793971 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,14 +9,18 @@ * 4. Prints grid information and first/last nodes */ -#include "user_input.h" -#include "types.hpp" #include "grid.hpp" -#include +#include "method.hpp" +#include "solver.hpp" +#include "types.hpp" +#include "user_input.h" #include +#include -int main() { - try { +int main() +{ + try + { // --- Physical constants (fixed from assignment) --- PhysParams phys; std::cout << "=== Grid initialization test ===\n"; @@ -24,10 +28,9 @@ int main() { // --- Ask user for dx --- double dx = askForDouble( - "Enter spatial step Δx [cm] (suggested 0.05): ", - [](double v){ return v > 0.0 && v <= 1.0; }, - "Δx must be positive and reasonably small (<= 1.0 cm)." - ); + "Enter spatial step Δx [cm] (suggested 0.05): ", + [](double v) { return v > 0.0 && v <= 1.0; }, + "Δx must be positive and reasonably small (<= 1.0 cm)."); // --- Create and validate grid --- Grid g(phys.L_cm, dx); @@ -40,20 +43,24 @@ int main() { std::cout << " Spatial step dx : " << g.dx << " cm\n"; std::cout << " Domain length L : " << g.L << " cm\n"; std::cout << " First node x[0] : " << g.x.front() << " cm\n"; - std::cout << " Last node x[N-1] : " << g.x.back() << " cm\n\n"; + std::cout << " Last node x[N-1] : " << g.x.back() << " cm\n\n"; + + // --- Run Crank-Nicolson Solver --- + std::cout << "=== Running Crank-Nicolson Solver ===\n"; + NumParams num; + num.dx = dx; + // Use default dt and tEnd from NumParams or ask user? + // For now use defaults: dt=0.01, tEnd=0.5 - // Optional: show first and last few nodes - std::cout << "First 5 nodes: "; - for (std::size_t i = 0; i < std::min(5, g.Nx); ++i) - std::cout << g.x[i] << " "; - std::cout << "\nLast 5 nodes: "; - for (std::size_t i = g.Nx - std::min(5, g.Nx); i < g.Nx; ++i) - std::cout << g.x[i] << " "; - std::cout << "\n"; + HeatSolver solver(phys, num, SchemeKind::CrankNicolson); + solver.run([](double t, const Grid& g, const std::vector& T) + { std::cout << "Time: " << t << " h, T[center] = " << T[g.Nx / 2] << "\n"; }); + std::cout << "Simulation finished.\n"; return 0; } - catch (const std::exception& e) { + catch (const std::exception& e) + { std::cerr << "Error: " << e.what() << "\n"; return 1; } diff --git a/src/method.cpp b/src/method.cpp index 91a6ae4..a97fa18 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -1,16 +1,85 @@ +/** + * @file method.cpp + * @brief Factory implementation. + */ + #include "method.hpp" -#include + +#include "methods/crank_nicolson.hpp" +// #include "methods/laasonen.hpp" #include +#include + +// --- Wrapper pour Laasonen --- +// class LaasonenMethod : public Method +// { +// public: +// std::string name() const override { return "Laasonen"; } +// bool uses_previous_step() const noexcept override +// { +// return false; +// } +// +// void step(const Grid& g, +// double D, +// double dt, +// const std::vector& /*Tnm1*/, +// const std::vector& Tn, +// std::vector& Tnp1) const override +// { +// // impl_.step(g, D, dt, Tnp1, Tn); +// } +// +// private: +// // Laasonen impl_; +// }; + +// --- Wrapper for Crank-Nicolson --- +class CrankNicolsonMethod : public Method +{ + public: + std::string name() const override + { + return "Crank-Nicolson"; + } + + bool uses_previous_step() const noexcept override + { + return false; + } + + void step(const Grid& g, + double D, + double dt, + const std::vector& /*Tnm1*/, + const std::vector& Tn, + std::vector& Tnp1) override + { + impl_.step(g, D, dt, Tnp1, Tn); + } + private: + CrankNicolson impl_; +}; + +// --- Factory --- std::unique_ptr make_method(SchemeKind scheme) { - // For now, no concrete method is available in this branch. - // The switch is kept so that the code remains aligned with the enum in the header. - switch (scheme) { + switch (scheme) + { + // case SchemeKind::Richardson: + // throw std::invalid_argument("make_method: Richardson not implemented yet"); + + // case SchemeKind::Laasonen: + // return std::make_unique(); + + case SchemeKind::CrankNicolson: + return std::make_unique(); + case SchemeKind::DuFortFrankel: - throw std::invalid_argument( - "make_method: DuFort-Frankel is not compiled in this branch"); - } + throw std::invalid_argument("make_method: DuFort-Frankel not implemented yet"); - throw std::invalid_argument("make_method: unsupported scheme"); -} + default: + throw std::invalid_argument("make_method: unknown scheme"); + } +} \ No newline at end of file diff --git a/src/methods/crank_nicolson.cpp b/src/methods/crank_nicolson.cpp new file mode 100644 index 0000000..8b15420 --- /dev/null +++ b/src/methods/crank_nicolson.cpp @@ -0,0 +1,73 @@ +/** + * @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, + std::vector& T, + const std::vector& Tn) const +{ + 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; + + // 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 * Tn[i - 1] + (1.0 - r) * Tn[i] + alpha * Tn[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 * T[0]; + d[M - 1] += alpha * T[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) { + T[j + 1] = x[j]; + } +} \ No newline at end of file From 7ef6fb1f87d978fd1b5535a298b4f8a62241fed6 Mon Sep 17 00:00:00 2001 From: Broky64 Date: Mon, 1 Dec 2025 14:42:47 +0100 Subject: [PATCH 2/2] feat: Refactor Crank-Nicolson to directly implement Method interface and add CSV output for simulation results. --- include/methods/crank_nicolson.hpp | 26 ++++++++++------ src/io.cpp | 7 +++-- src/main.cpp | 36 ++++++++++++++++++--- src/method.cpp | 36 +-------------------- src/methods/crank_nicolson.cpp | 50 ++++++++++++++++++------------ 5 files changed, 83 insertions(+), 72 deletions(-) diff --git a/include/methods/crank_nicolson.hpp b/include/methods/crank_nicolson.hpp index 97892e5..5129e90 100644 --- a/include/methods/crank_nicolson.hpp +++ b/include/methods/crank_nicolson.hpp @@ -1,25 +1,31 @@ #pragma once -#include "grid.hpp" -#include +#include "method.hpp" +#include /** * @brief Crank-Nicolson scheme implementation. */ -class CrankNicolson +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. - * @param g Grid metadata - * @param D Diffusivity - * @param dt Time step - * @param T Output container for T^{n+1} - * @param Tn Input container for T^{n} */ void step(const Grid& g, double D, double dt, - std::vector& T, - const std::vector& Tn) const; + const std::vector& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) override; }; \ No newline at end of file diff --git a/src/io.cpp b/src/io.cpp index 3b0dd9f..e52904a 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -1,8 +1,9 @@ #include "io.hpp" -#include + #include -#include +#include #include +#include namespace fs = std::filesystem; @@ -24,7 +25,7 @@ void save_profile_csv(const std::string& path, throw std::runtime_error("save_profile_csv: vector size mismatch"); // --- Ensure the result directory exists --- - fs::path fullPath = fs::path("result") / path; + fs::path fullPath = path; fs::create_directories(fullPath.parent_path()); // --- Open file for writing --- diff --git a/src/main.cpp b/src/main.cpp index 5793971..7e7a775 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -10,12 +10,15 @@ */ #include "grid.hpp" +#include "io.hpp" #include "method.hpp" #include "solver.hpp" #include "types.hpp" #include "user_input.h" +#include #include #include +#include int main() { @@ -49,13 +52,36 @@ int main() std::cout << "=== Running Crank-Nicolson Solver ===\n"; NumParams num; num.dx = dx; - // Use default dt and tEnd from NumParams or ask user? - // For now use defaults: dt=0.01, tEnd=0.5 + num.dt = 0.01; + num.tEnd = 0.5; + num.outEvery = 0.1; HeatSolver solver(phys, num, SchemeKind::CrankNicolson); - solver.run([](double t, const Grid& g, const std::vector& T) - { std::cout << "Time: " << t << " h, T[center] = " << T[g.Nx / 2] << "\n"; }); - std::cout << "Simulation finished.\n"; + + // Prepare output directory + // Note: io.cpp in this branch might not handle directory creation automatically if it's + // old. But let's assume standard behavior or check io.cpp. Actually, let's just use a + // simple callback that prints AND saves if possible. Since I haven't checked io.cpp in this + // branch, I'll rely on a lambda that constructs the path. + + std::string outDir = "results/Crank-Nicolson"; + // Create directory command (quick hack for cross-platform C++17/20) + std::filesystem::create_directories(outDir); + + solver.run( + [&](double t, const Grid& g, const std::vector& T) + { + std::cout << "Time: " << t << " h, T[center] = " << T[g.Nx / 2] << "\n"; + + std::ostringstream filename; + filename << outDir << "/profile_t_" << std::fixed << std::setprecision(2) << t + << "h.csv"; + + // Assuming save_profile_csv exists in io.hpp + save_profile_csv(filename.str(), g.x, T); + }); + + std::cout << "Simulation finished. Results saved to " << outDir << "\n"; return 0; } diff --git a/src/method.cpp b/src/method.cpp index a97fa18..4a264f2 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -34,47 +34,13 @@ // // Laasonen impl_; // }; -// --- Wrapper for Crank-Nicolson --- -class CrankNicolsonMethod : public Method -{ - public: - std::string name() const override - { - return "Crank-Nicolson"; - } - - bool uses_previous_step() const noexcept override - { - return false; - } - - void step(const Grid& g, - double D, - double dt, - const std::vector& /*Tnm1*/, - const std::vector& Tn, - std::vector& Tnp1) override - { - impl_.step(g, D, dt, Tnp1, Tn); - } - - private: - CrankNicolson impl_; -}; - // --- Factory --- std::unique_ptr make_method(SchemeKind scheme) { switch (scheme) { - // case SchemeKind::Richardson: - // throw std::invalid_argument("make_method: Richardson not implemented yet"); - - // case SchemeKind::Laasonen: - // return std::make_unique(); - case SchemeKind::CrankNicolson: - return std::make_unique(); + return std::make_unique(); case SchemeKind::DuFortFrankel: throw std::invalid_argument("make_method: DuFort-Frankel not implemented yet"); diff --git a/src/methods/crank_nicolson.cpp b/src/methods/crank_nicolson.cpp index 8b15420..7c0c72a 100644 --- a/src/methods/crank_nicolson.cpp +++ b/src/methods/crank_nicolson.cpp @@ -4,24 +4,32 @@ */ #include "methods/crank_nicolson.hpp" + #include #include void CrankNicolson::step(const Grid& g, double D, double dt, - std::vector& T, - const std::vector& Tn) const + 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; + const double r = D * dt / (dx * dx); + const double alpha = r / 2.0; // Safety check for grid size - if (N < 3) return; + 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; @@ -29,29 +37,31 @@ void CrankNicolson::step(const Grid& g, // --- 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 + 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) { + for (std::size_t j = 0; j < M; ++j) + { std::size_t i = j + 1; // Actual grid index - d[j] = alpha * Tn[i - 1] + (1.0 - r) * Tn[i] + alpha * Tn[i + 1]; + 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 * T[0]; - d[M - 1] += alpha * T[N - 1]; + 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) { + 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]; @@ -61,13 +71,15 @@ void CrankNicolson::step(const Grid& g, std::vector x(M); x[M - 1] = d[M - 1] / b[M - 1]; - for (std::size_t i = M - 1; i-- > 0; ) { + 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) { - T[j + 1] = x[j]; + for (std::size_t j = 0; j < M; ++j) + { + Tnext[j + 1] = x[j]; } } \ No newline at end of file