From 33aae72c56ea411f4c13d9c76d91f83e9e914117 Mon Sep 17 00:00:00 2001 From: Broky64 Date: Wed, 5 Nov 2025 17:09:07 +0100 Subject: [PATCH 1/3] Implement Laasonen implicit scheme for 1D heat equation solver - Refactor method interface to include Laasonen scheme - Add Laasonen class with time-stepping logic - Update factory method to create LaasonenMethod instance - Modify HeatSolver to support Laasonen scheme execution - Enhance main program to test Laasonen scheme and save results as CSV --- include/method.hpp | 73 ++++++++-------------- include/methods/laasonen.hpp | 26 ++++++++ src/main.cpp | 115 +++++++++++++++++++++++------------ src/method.cpp | 56 ++++++++++++++--- src/methods/laasonen.cpp | 58 ++++++++++++++++++ src/solver.cpp | 58 +++++++----------- 6 files changed, 257 insertions(+), 129 deletions(-) create mode 100644 include/methods/laasonen.hpp create mode 100644 src/methods/laasonen.cpp diff --git a/include/method.hpp b/include/method.hpp index 5488170..2633172 100644 --- a/include/method.hpp +++ b/include/method.hpp @@ -1,60 +1,39 @@ -#pragma once -/** - * @file method.hpp - * @brief Numerical scheme interface for the 1D heat equation. - */ +#ifndef METHOD_HPP +#define METHOD_HPP -#include "grid.hpp" -#include #include -#include +#include +#include "grid.hpp" -/** - * @brief Supported numerical schemes. - */ +// Tous les schémas que tu gères dans ton prog enum class SchemeKind { - DuFortFrankel - // Richardson, Laasonen, CrankNicolson (to be added later) + Richardson, + DuFortFrankel, + Laasonen }; -/** - * @brief Abstract base class for time-stepping methods. - */ +// Interface commune à tous les schémas class Method { public: virtual ~Method() = default; - /// @return human-readable scheme name. - virtual std::string name() const = 0; - - /** - * @brief Advance one time step, computing T^{n+1} from T^n (and possibly T^{n-1}). - * @param g Grid metadata (spacing and node count) - * @param D Diffusivity [cm^2/h] - * @param dt Time step [h] - * @param Tprev Read-only temperature at the previous time layer T^{n-1} - * (may be empty when the scheme is two-level) - * @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, - const std::vector& Tprev, - const std::vector& Tcurr, - std::vector& Tnext) = 0; + // true → schéma à 3 niveaux (a besoin de T^{n-1}) + // false → schéma à 2 niveaux (juste T^{n}) + virtual bool uses_previous_step() const = 0; - /** - * @brief Indicate whether the scheme needs the previous time layer. - * - * 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; } + // Appelé par le solver à chaque pas de temps + // Tnm1 : T^{n-1} (inutile pour les schémas à 2 niveaux) + // Tn : T^{n} + // Tnp1 : T^{n+1} (à remplir) + virtual void step(const Grid& g, + double D, + double dt, + const std::vector& Tnm1, + const std::vector& Tn, + std::vector& Tnp1) const = 0; }; -/** - * @brief Factory helper that builds a Method instance for the requested scheme. - * @param scheme Scheme identifier (e.g. SchemeKind::DuFortFrankel) - * @return Unique pointer owning the concrete method implementation - * @throws std::invalid_argument if the scheme is not supported - */ -[[nodiscard]] std::unique_ptr make_method(SchemeKind scheme); +// Fabrique qui crée le bon schéma +std::unique_ptr make_method(SchemeKind scheme); + +#endif // METHOD_HPP diff --git a/include/methods/laasonen.hpp b/include/methods/laasonen.hpp new file mode 100644 index 0000000..5e5f683 --- /dev/null +++ b/include/methods/laasonen.hpp @@ -0,0 +1,26 @@ +#ifndef LAASONEN_HPP +#define LAASONEN_HPP + +#include +#include "grid.hpp" + +// ----------------------------------------------------------------------------- +// Méthode implicite de Laasonen (Backward Euler en temps + central en espace) +// Inconditionnellement stable pour l'équation de diffusion 1D +// ----------------------------------------------------------------------------- +class Laasonen { +public: + // Avance d’un pas de temps : + // - g : maillage (Nx, dx) + // - D : coefficient de diffusion + // - dt : pas de temps + // - T : vecteur de température (sera modifié → T^{n+1}) + // - Tprev : température à l’étape précédente T^n + void step(const Grid& g, + double D, + double dt, + std::vector& T, + const std::vector& Tprev) const; +}; + +#endif // LAASONEN_HPP diff --git a/src/main.cpp b/src/main.cpp index 08d4330..c65e2d3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,60 +1,95 @@ /** * @file main.cpp - * @brief Simple test program for the Grid and user input utilities. - * - * This program: - * 1. Asks the user for Δx (spatial step) - * 2. Builds a uniform grid from 0 to L = 31 cm - * 3. Validates its consistency (assertions in Debug mode) - * 4. Prints grid information and first/last nodes + * @brief Test driver for the 1D heat conduction solver (Laasonen implicit scheme), + * saving temperature profiles as CSV files. */ -#include "user_input.h" #include "types.hpp" #include "grid.hpp" +#include "solver.hpp" +#include "method.hpp" +#include "user_input.h" +#include "io.hpp" // <-- pour save_profile_csv() + #include #include +#include +#include +#include + +/** + * @brief Output callback: save temperature profiles to CSV at each dump time. + * + * The output folder is automatically created under `result/laasonen/`. + * Example files: `result/laasonen/t_0.00.csv`, `t_0.10.csv`, etc. + * + * @param t Current time [h] + * @param g Grid object (for x-locations) + * @param T Temperature field [°C] + */ +void save_to_csv(double t, const Grid& g, const std::vector& T) +{ + namespace fs = std::filesystem; + + const std::string folder = "result/laasonen"; + if (!fs::exists(folder)) { + fs::create_directories(folder); + } -int main() { + // Build file name like "t_0.10.csv" + std::ostringstream name; + name << folder << "/t_" << std::fixed << std::setprecision(2) << t << ".csv"; + + // Save profile to CSV using helper from io.cpp + save_profile_csv(name.str(), g.x, T); + + std::cout << "Saved: " << name.str() << "\n"; +} + +/** + * @brief Entry point for testing the Laasonen implicit scheme. + */ +int main() +{ try { - // --- Physical constants (fixed from assignment) --- + std::cout << "=== 1D Heat Equation Solver (Laasonen test) ===\n"; + + // --- Physical parameters (assignment constants) --- PhysParams phys; - std::cout << "=== Grid initialization test ===\n"; - std::cout << "Domain length L = " << phys.L_cm << " cm\n"; - - // --- 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)." - ); - - // --- Create and validate grid --- - Grid g(phys.L_cm, dx); + phys.L_cm = 31.0; ///< Wall thickness [cm] + phys.Tin = 38.0; ///< Initial temperature [°C] + phys.Tsur = 149.0; ///< Boundary temperature [°C] + phys.D_cm2h = 93.0; ///< Diffusivity [cm²/h] + + // --- Numerical parameters --- + NumParams num; + 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; ///< Simulation duration [h] + num.outEvery = 0.1; ///< Output interval [h] + + // --- Build grid --- + Grid g(phys.L_cm, num.dx); validate_grid(g); - // --- Output basic info --- - std::cout << std::fixed << std::setprecision(4); - std::cout << "\nGrid successfully created ✅\n"; - std::cout << " Number of nodes Nx : " << g.Nx << "\n"; - 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"; - - // 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"; + std::cout << "\nGrid created ✅ Nx = " << g.Nx + << ", dx = " << g.dx << " cm\n"; + + // --- Create solver with Laasonen scheme --- + HeatSolver solver(phys, num, SchemeKind::Laasonen); + std::cout << "Running Laasonen implicit scheme...\n"; + + // --- Run simulation and save results --- + solver.run(save_to_csv); + std::cout << "\nSimulation completed successfully ✅\n"; + std::cout << "Results available under ./result/laasonen/\n"; return 0; } catch (const std::exception& e) { - std::cerr << "Error: " << e.what() << "\n"; + std::cerr << "Fatal error: " << e.what() << "\n"; return 1; } } diff --git a/src/method.cpp b/src/method.cpp index 91a6ae4..145b3d2 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -1,16 +1,58 @@ +/** + * @file method.cpp + * @brief Factory for numerical schemes. + */ + #include "method.hpp" +#include "methods/laasonen.hpp" // <-- ce qu'on vient d’ajouter #include #include +// ----------------------------------------------------------------------------- +// Adaptateur pour utiliser ta classe Laasonen (qui a step(..., T, Tprev)) +// dans l’interface commune Method +// ----------------------------------------------------------------------------- +class LaasonenMethod : public Method { +public: + bool uses_previous_step() const override { + // Laasonen est implicite à 2 niveaux → pas besoin de T^{n-1} + return false; + } + + void step(const Grid& g, + double D, + double dt, + const std::vector& /*Tnm1*/, + const std::vector& Tn, + std::vector& Tnp1) const override + { + // on délègue au Laasonen "pur" qu’on a écrit + impl_.step(g, D, dt, Tnp1, Tn); + } + +private: + Laasonen impl_; +}; + +/** + * @brief Build a concrete numerical method from the requested scheme. + */ 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) { - case SchemeKind::DuFortFrankel: - throw std::invalid_argument( - "make_method: DuFort-Frankel is not compiled in this branch"); - } + case SchemeKind::Richardson: + throw std::invalid_argument( + "make_method: Richardson is not compiled in this branch"); + + case SchemeKind::Laasonen: + return std::make_unique(); - throw std::invalid_argument("make_method: unsupported scheme"); + case SchemeKind::DuFortFrankel: + // tu pourras remplacer ça quand ta classe DuFort-Frankel sera prête + throw std::invalid_argument( + "make_method: DuFort-Frankel is not compiled in this branch"); + + default: + throw std::invalid_argument("make_method: unknown scheme"); + } } diff --git a/src/methods/laasonen.cpp b/src/methods/laasonen.cpp new file mode 100644 index 0000000..f90919a --- /dev/null +++ b/src/methods/laasonen.cpp @@ -0,0 +1,58 @@ +#include "methods/laasonen.hpp" +#include +#include + +void Laasonen::step(const Grid& g, + double D, + double dt, + std::vector& T, + const std::vector& Tprev) const +{ + const std::size_t N = g.Nx; + const double dx = g.dx; + const double r = D * dt / (dx * dx); + + // Sécurité + if (N < 3) return; + + const std::size_t M = N - 2; // nombre de points internes (sans BC) + + // Coefficients tridiagonaux + std::vector a(M, -r); + std::vector b(M, 1.0 + 2.0 * r); + std::vector c(M, -r); + std::vector d(M, 0.0); + + // Construction du RHS (Right-Hand Side) + for (std::size_t j = 0; j < M; ++j) { + d[j] = Tprev[j + 1]; + } + + // Ajout des conditions aux limites dans le RHS + d[0] += r * T[0]; // dépend de T[0] (bord gauche) + d[M-1] += r * T[N - 1]; // dépend de T[N-1] (bord droit) + + // ----------------------- + // Algorithme de Thomas + // ----------------------- + + // Forward elimination + 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 + std::vector x(M, 0.0); + 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]; + } + + // Mise à jour de T^{n+1} (les points internes) + for (std::size_t j = 0; j < M; ++j) { + T[j + 1] = x[j]; + } +} diff --git a/src/solver.cpp b/src/solver.cpp index f0e5f86..3489525 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -1,81 +1,71 @@ +/** + * @file solver.cpp + * @brief Time-marching driver for the 1D heat equation. + */ + #include "solver.hpp" -#include +#include // 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 computation method (factory according to your code) + // on récupère la bonne méthode via la factory (définie dans method.cpp) method_ = make_method(scheme); uses_prev_step_ = method_->uses_previous_step(); - // Allocate vectors + // allocation T_.resize(grid_.Nx); next_.resize(grid_.Nx); if (uses_prev_step_) { Tprev_.resize(grid_.Nx); } - // Initialize the field + // init init_field(); - - // Bootstrapping for three-level methods (such as DuFort–Frankel) bootstrap_if_needed(); } void HeatSolver::apply_dirichlet(std::vector& T) const { - // Impose surface temperature at both ends T.front() = phys_.Tsur; T.back() = phys_.Tsur; } void HeatSolver::init_field() { - // Uniform initial field at Tin std::fill(T_.begin(), T_.end(), phys_.Tin); - - // Boundary conditions apply_dirichlet(T_); } void HeatSolver::bootstrap_if_needed() { - if (!uses_prev_step_) { + if (!uses_prev_step_) return; - } - // Here, T_ = T^0 - - // 1) Store T^0 in Tprev_ + // T_ = T^0 Tprev_ = T_; - // 2) Compute T^1 using a small FTCS step (1D diffusion) const double dx = grid_.dx; 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) + // petit FTCS pour construire T^1 if (r <= 0.5) { for (std::size_t i = 1; i + 1 < N; ++i) { - next_[i] = T_[i] + r * (T_[i + 1] - 2.0 * T_[i] + T_[i - 1]); + next_[i] = T_[i] + r * (T_[i+1] - 2.0 * T_[i] + T_[i-1]); } } else { - // If dt is too large, keep T^1 = T^0 to avoid oscillations next_ = T_; } - // Boundary conditions apply_dirichlet(next_); - - // 3) Update T_ to contain T^1 - T_ = next_; + T_ = next_; // T^1 } void HeatSolver::run( @@ -85,35 +75,33 @@ void HeatSolver::run( const double dt = num_.dt; const double tEnd = num_.tEnd; const double outEvery = num_.outEvery; + double nextDump = 0.0; - double nextDump = 0.0; - - // Initial dump + // sortie initiale onDump(t, grid_, T_); - // Main time loop while (t < tEnd - 1e-12) { if (uses_prev_step_) { - // Three-level method (e.g., DuFort–Frankel) + // schémas à 3 niveaux (ex. Richardson) method_->step(grid_, phys_.D_cm2h, dt, Tprev_, T_, next_); } else { - // Two-level method (e.g., FTCS, BTCS…) + // schémas à 2 niveaux (ex. Laasonen implicite) method_->step(grid_, phys_.D_cm2h, dt, T_, T_, next_); } - // Apply boundary conditions + // BC apply_dirichlet(next_); - // Shift time states + // shift if (uses_prev_step_) { Tprev_ = T_; } T_ = next_; - // Advance time + // temps t += dt; - // Periodic output + // sorties périodiques if (t + 1e-12 >= nextDump) { onDump(t, grid_, T_); nextDump += outEvery; From 88b7a51a602ed3eb9063ec854d87f6bebd80476c Mon Sep 17 00:00:00 2001 From: Broky64 Date: Mon, 1 Dec 2025 08:32:53 +0100 Subject: [PATCH 2/3] fix: init for laasonen was false for N+1 ans Nmax-1 --- .gitignore | 2 +- src/solver.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 2787c44..c014479 100644 --- a/.gitignore +++ b/.gitignore @@ -14,7 +14,7 @@ .idea/ # Outputs -/results/* +/result/* !.gitkeep # MacOS diff --git a/src/solver.cpp b/src/solver.cpp index 3489525..3dd30b5 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -80,12 +80,12 @@ void HeatSolver::run( // sortie initiale onDump(t, grid_, T_); + apply_dirichlet(next_); + while (t < tEnd - 1e-12) { if (uses_prev_step_) { - // schémas à 3 niveaux (ex. Richardson) method_->step(grid_, phys_.D_cm2h, dt, Tprev_, T_, next_); } else { - // schémas à 2 niveaux (ex. Laasonen implicite) method_->step(grid_, phys_.D_cm2h, dt, T_, T_, next_); } From b9fc9365ec0b2f1235eed2d714786d6c6d4e5aeb Mon Sep 17 00:00:00 2001 From: Broky64 Date: Mon, 1 Dec 2025 12:02:27 +0100 Subject: [PATCH 3/3] feat: Refactor Laasonen to directly implement Method interface, update its step signature and Thomas algorithm, and add method naming. --- .gitignore | 6 ++- include/method.hpp | 15 +++++-- include/methods/laasonen.hpp | 42 +++++++++--------- src/io.cpp | 7 +-- src/main.cpp | 57 +++++++++++++------------ src/method.cpp | 56 +++++++----------------- src/methods/laasonen.cpp | 82 +++++++++++++++++++++++++++--------- 7 files changed, 153 insertions(+), 112 deletions(-) diff --git a/.gitignore b/.gitignore index c014479..6809f88 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ # IDE .vscode/ .idea/ +.cache/ # Outputs /result/* @@ -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 2633172..e80f76e 100644 --- a/include/method.hpp +++ b/include/method.hpp @@ -1,22 +1,29 @@ #ifndef METHOD_HPP #define METHOD_HPP +#include "grid.hpp" #include +#include #include -#include "grid.hpp" // Tous les schémas que tu gères dans ton prog -enum class SchemeKind { +enum class SchemeKind +{ Richardson, DuFortFrankel, Laasonen + // CrankNicolson (to be added later) }; // Interface commune à tous les schémas -class Method { -public: +class Method +{ + public: virtual ~Method() = default; + /// @return human-readable scheme name. + virtual std::string name() const = 0; + // true → schéma à 3 niveaux (a besoin de T^{n-1}) // false → schéma à 2 niveaux (juste T^{n}) virtual bool uses_previous_step() const = 0; diff --git a/include/methods/laasonen.hpp b/include/methods/laasonen.hpp index 5e5f683..fdebdff 100644 --- a/include/methods/laasonen.hpp +++ b/include/methods/laasonen.hpp @@ -1,26 +1,30 @@ -#ifndef LAASONEN_HPP -#define LAASONEN_HPP +#pragma once +#include "method.hpp" +#include #include -#include "grid.hpp" -// ----------------------------------------------------------------------------- -// Méthode implicite de Laasonen (Backward Euler en temps + central en espace) -// Inconditionnellement stable pour l'équation de diffusion 1D -// ----------------------------------------------------------------------------- -class Laasonen { -public: - // Avance d’un pas de temps : - // - g : maillage (Nx, dx) - // - D : coefficient de diffusion - // - dt : pas de temps - // - T : vecteur de température (sera modifié → T^{n+1}) - // - Tprev : température à l’étape précédente T^n +/** + * @brief Laasonen implicit method (Backward Euler in time + Central in space). + * Unconditionally stable for the 1D heat equation. + */ +class Laasonen : public Method +{ + public: + std::string name() const override + { + return "Laasonen"; + } + void step(const Grid& g, double D, double dt, - std::vector& T, - const std::vector& Tprev) const; -}; + const std::vector& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) const override; -#endif // LAASONEN_HPP + bool uses_previous_step() const noexcept override + { + return false; + } +}; 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 c65e2d3..3dc8b88 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,18 +4,17 @@ * saving temperature profiles as CSV files. */ -#include "types.hpp" #include "grid.hpp" -#include "solver.hpp" +#include "io.hpp" // <-- pour save_profile_csv() #include "method.hpp" +#include "solver.hpp" +#include "types.hpp" #include "user_input.h" -#include "io.hpp" // <-- pour save_profile_csv() - -#include +#include #include +#include #include #include -#include /** * @brief Output callback: save temperature profiles to CSV at each dump time. @@ -32,7 +31,8 @@ void save_to_csv(double t, const Grid& g, const std::vector& T) namespace fs = std::filesystem; const std::string folder = "result/laasonen"; - if (!fs::exists(folder)) { + if (!fs::exists(folder)) + { fs::create_directories(folder); } @@ -51,31 +51,35 @@ void save_to_csv(double t, const Grid& g, const std::vector& T) */ int main() { - try { + try + { std::cout << "=== 1D Heat Equation Solver (Laasonen test) ===\n"; // --- Physical parameters (assignment constants) --- PhysParams phys; - phys.L_cm = 31.0; ///< Wall thickness [cm] - phys.Tin = 38.0; ///< Initial temperature [°C] - phys.Tsur = 149.0; ///< Boundary temperature [°C] - phys.D_cm2h = 93.0; ///< Diffusivity [cm²/h] + phys.L_cm = 31.0; ///< Wall thickness [cm] + phys.Tin = 38.0; ///< Initial temperature [°C] + phys.Tsur = 149.0; ///< Boundary temperature [°C] + phys.D_cm2h = 93.0; ///< Thermal diffusivity [cm^2/h] - // --- Numerical parameters --- - NumParams num; - 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; ///< Simulation duration [h] - num.outEvery = 0.1; ///< Output interval [h] + // --- Method Selection --- + std::cout << "\nSelect numerical method:\n"; + std::cout << "1. Laasonen\n"; + // std::cout << "2. Richardson (Not available)\n"; + // std::cout << "3. DuFort-Frankel (Not available)\n"; - // --- Build grid --- - Grid g(phys.L_cm, num.dx); - validate_grid(g); + int choice = static_cast(askForDouble( + "Enter choice (1): ", [](double v) { return v == 1.0; }, "Please enter 1.")); - std::cout << "\nGrid created ✅ Nx = " << g.Nx - << ", dx = " << g.dx << " cm\n"; + // --- Numerical parameters --- + NumParams num; + 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; ///< Simulation duration [h] + num.outEvery = 0.1; ///< Output interval [h] // --- Create solver with Laasonen scheme --- HeatSolver solver(phys, num, SchemeKind::Laasonen); @@ -88,7 +92,8 @@ int main() std::cout << "Results available under ./result/laasonen/\n"; return 0; } - catch (const std::exception& e) { + catch (const std::exception& e) + { std::cerr << "Fatal error: " << e.what() << "\n"; return 1; } diff --git a/src/method.cpp b/src/method.cpp index 145b3d2..79d79d2 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -4,55 +4,31 @@ */ #include "method.hpp" -#include "methods/laasonen.hpp" // <-- ce qu'on vient d’ajouter -#include -#include - -// ----------------------------------------------------------------------------- -// Adaptateur pour utiliser ta classe Laasonen (qui a step(..., T, Tprev)) -// dans l’interface commune Method -// ----------------------------------------------------------------------------- -class LaasonenMethod : public Method { -public: - bool uses_previous_step() const override { - // Laasonen est implicite à 2 niveaux → pas besoin de T^{n-1} - return false; - } - - void step(const Grid& g, - double D, - double dt, - const std::vector& /*Tnm1*/, - const std::vector& Tn, - std::vector& Tnp1) const override - { - // on délègue au Laasonen "pur" qu’on a écrit - impl_.step(g, D, dt, Tnp1, Tn); - } -private: - Laasonen impl_; -}; +// #include "methods/dufort_frankel.hpp" // Not available +#include "methods/laasonen.hpp" +// #include "methods/richardson.hpp" // Not available +#include +#include /** * @brief Build a concrete numerical method from the requested scheme. */ std::unique_ptr make_method(SchemeKind scheme) { - switch (scheme) { - case SchemeKind::Richardson: - throw std::invalid_argument( - "make_method: Richardson is not compiled in this branch"); + switch (scheme) + { + case SchemeKind::Laasonen: + return std::make_unique(); - case SchemeKind::Laasonen: - return std::make_unique(); + case SchemeKind::Richardson: + throw std::invalid_argument("make_method: Richardson is not compiled in this branch"); - case SchemeKind::DuFortFrankel: - // tu pourras remplacer ça quand ta classe DuFort-Frankel sera prête - throw std::invalid_argument( - "make_method: DuFort-Frankel is not compiled in this branch"); + case SchemeKind::DuFortFrankel: + throw std::invalid_argument( + "make_method: DuFort-Frankel is not compiled in this branch"); - default: - throw std::invalid_argument("make_method: unknown scheme"); + default: + throw std::invalid_argument("make_method: unknown scheme"); } } diff --git a/src/methods/laasonen.cpp b/src/methods/laasonen.cpp index f90919a..8022c5d 100644 --- a/src/methods/laasonen.cpp +++ b/src/methods/laasonen.cpp @@ -1,58 +1,102 @@ #include "methods/laasonen.hpp" + #include #include void Laasonen::step(const Grid& g, double D, double dt, - std::vector& T, - const std::vector& Tprev) const + const std::vector& /*Tprev*/, + const std::vector& Tcurr, + std::vector& Tnext) const { const std::size_t N = g.Nx; const double dx = g.dx; - const double r = D * dt / (dx * dx); + const double r = D * dt / (dx * dx); // Sécurité - if (N < 3) return; + if (N < 3) + return; + + // Ensure output size + if (Tnext.size() != N) + Tnext.resize(N); + + // Copy boundary conditions from current (or apply later, but solver handles BC usually) + // The solver (HeatSolver) applies BCs after step(), but we need them for the system if we solve + // for internal nodes. Actually, Tnext will be overwritten for internal nodes. We can just copy + // Tcurr to Tnext to initialize (and keep boundaries if they don't change in this step, though + // solver overrides them). + Tnext = Tcurr; const std::size_t M = N - 2; // nombre de points internes (sans BC) // Coefficients tridiagonaux + // Equation: -r T_{i-1}^{n+1} + (1 + 2r) T_i^{n+1} - r T_{i+1}^{n+1} = T_i^n std::vector a(M, -r); std::vector b(M, 1.0 + 2.0 * r); std::vector c(M, -r); std::vector d(M, 0.0); - // Construction du RHS (Right-Hand Side) - for (std::size_t j = 0; j < M; ++j) { - d[j] = Tprev[j + 1]; + // Construction du RHS (Right-Hand Side) = T^n + for (std::size_t j = 0; j < M; ++j) + { + d[j] = Tcurr[j + 1]; } // Ajout des conditions aux limites dans le RHS - d[0] += r * T[0]; // dépend de T[0] (bord gauche) - d[M-1] += r * T[N - 1]; // dépend de T[N-1] (bord droit) + // T_{i-1}^{n+1} for i=1 is T_0^{n+1} (Boundary Left) + // T_{i+1}^{n+1} for i=M is T_{N-1}^{n+1} (Boundary Right) + // We assume Dirichlet BCs are fixed or handled by the solver. + // The solver calls apply_dirichlet(next_) AFTER step(). + // But for the implicit solve, we need the boundary values at n+1. + // Since we don't have them passed explicitly as "next boundaries", we usually assume they are + // same as Tcurr or fixed. In this problem, Tsur is constant. So Tcurr[0] and Tcurr[N-1] are + // correct for T^{n+1} boundaries too. + + d[0] += r * Tcurr[0]; // T_0^{n+1} assumed approx T_0^n or fixed + d[M - 1] += r * Tcurr[N - 1]; // T_{N-1}^{n+1} // ----------------------- // Algorithme de Thomas // ----------------------- // Forward elimination - 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]; + // c[0] /= b[0]; + // d[0] /= b[0]; + // for i = 1 to M-1 + // temp = 1 / (b[i] - a[i] * c[i-1]) + // c[i] *= temp + // d[i] = (d[i] - a[i] * d[i-1]) * temp + + // Standard Thomas implementation + // Modify coefficients in place + c[0] /= b[0]; + d[0] /= b[0]; + + for (std::size_t i = 1; i < M; ++i) + { + double temp = 1.0 / (b[i] - a[i] * c[i - 1]); + c[i] *= temp; + d[i] = (d[i] - a[i] * d[i - 1]) * temp; } // Backward substitution - std::vector x(M, 0.0); - x[M - 1] = d[M - 1] / b[M - 1]; + // x[M-1] = d[M-1] + // for i = M-2 down to 0 + // x[i] = d[i] - c[i] * x[i+1] + + std::vector x(M); + x[M - 1] = d[M - 1]; - for (std::size_t i = M - 1; i-- > 0; ) { - x[i] = (d[i] - c[i] * x[i + 1]) / b[i]; + for (int i = static_cast(M) - 2; i >= 0; --i) + { + x[i] = d[i] - c[i] * x[i + 1]; } // Mise à jour de T^{n+1} (les points internes) - 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]; } }