diff --git a/.gitignore b/.gitignore index 649aefa..154bef8 100644 --- a/.gitignore +++ b/.gitignore @@ -15,7 +15,6 @@ .cache/ # Outputs -/result/* /results/* !.gitkeep @@ -26,5 +25,4 @@ .temp/ # Logs -*.log - +*.log \ No newline at end of file diff --git a/include/method.hpp b/include/method.hpp index 42646f7..bce34ad 100644 --- a/include/method.hpp +++ b/include/method.hpp @@ -1,68 +1,67 @@ -#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 /** - * @brief Supported numerical schemes. + * @brief Enumeration of the numerical schemes available in the solver. */ enum class SchemeKind { - DuFortFrankel, - Richardson - // Laasonen, CrankNicolson (to be added later) + Richardson, ///< Central time, central space (explicit, unstable). + DuFortFrankel, ///< Modified Richardson (explicit, stable). + Laasonen, ///< Simple Implicit (forward time, central space). + CrankNicolson ///< Trapezoidal Implicit (second order accuracy). }; /** - * @brief Abstract base class for time-stepping methods. + * @class Method + * @brief Abstract interface for time-integration schemes solving the 1D heat equation. + * + * This interface abstracts away the difference between explicit/implicit schemes + * and 2-level/3-level time stepping schemes. */ class Method { - public: +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} + * @brief Computes the temperature field at the next time step (T^{n+1}). + * + * @param g Grid metadata (spacing dx and node count Nx). + * @param D Thermal diffusivity [cm^2/h]. + * @param dt Time step size [h]. + * @param Tprev Temperature at time step n-1 (Read-only). + * Required only for 3-level schemes (Richardson, DuFort-Frankel). + * For 2-level schemes, this vector may be ignored. + * @param Tcurr Temperature at time step n (Read-only). + * @param Tnext Output vector to be filled with temperature at time step n+1. */ virtual void step(const Grid& g, double D, double dt, const std::vector& Tprev, const std::vector& Tcurr, - std::vector& Tnext) = 0; + std::vector& Tnext) const = 0; /** - * @brief Indicate whether the scheme needs the previous time layer. + * @brief Indicates if the scheme requires data from time step n-1. * - * 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. + * @return true For 3-level schemes (e.g., Richardson, DuFort-Frankel). + * @return false For 2-level schemes (e.g., Laasonen, Crank-Nicolson). */ - virtual bool uses_previous_step() const noexcept - { - return false; - } + virtual bool uses_previous_step() 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 + * @brief Factory function to create a concrete method instance. + * + * @param scheme The type of numerical scheme requested. + * @return std::unique_ptr Pointer to the created solver method. */ -[[nodiscard]] std::unique_ptr make_method(SchemeKind scheme); +std::unique_ptr make_method(SchemeKind scheme); + +#endif // METHOD_HPP \ No newline at end of file diff --git a/include/methods/laasonen.hpp b/include/methods/laasonen.hpp new file mode 100644 index 0000000..fdebdff --- /dev/null +++ b/include/methods/laasonen.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include "method.hpp" +#include +#include + +/** + * @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, + const std::vector& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) const override; + + bool uses_previous_step() const noexcept override + { + return false; + } +}; diff --git a/src/io.cpp b/src/io.cpp index 1732cef..f55c6e1 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -10,21 +10,23 @@ namespace fs = std::filesystem; /** * @brief Save a temperature profile (x, T) to a CSV file. * - * This function automatically creates the "result/" directory if missing. + * This function automatically creates the directory structure if missing. * Example output path: "result/dufort_frankel/t_0.10.csv" * - * @param path Relative path to the file (inside result/) - * @param x Spatial coordinates [cm] - * @param T Temperature values [°C] + * @param path Relative or absolute path to the file. + * @param x Spatial coordinates [cm]. + * @param T Temperature values [°C]. */ void save_profile_csv(const std::string& path, const std::vector& x, const std::vector& T) { + // Validate input sizes if (x.size() != T.size()) throw std::runtime_error("save_profile_csv: vector size mismatch"); - // --- Ensure the directory exists --- + // --- Ensure the directory exists (Resolved Conflict) --- + // We use the safer check from the 'dev' branch. fs::path fullPath(path); if (fullPath.has_parent_path()) { @@ -41,4 +43,4 @@ void save_profile_csv(const std::string& path, file << std::fixed << std::setprecision(6); for (std::size_t i = 0; i < x.size(); ++i) file << x[i] << "," << T[i] << "\n"; -} +} \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 4987d7a..6318292 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 and Richardson schemes. + * Supports DuFort-Frankel, Richardson, and Laasonen schemes. */ #include "grid.hpp" @@ -15,16 +15,21 @@ #include #include #include +#include +#include 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 + * This helper function avoids code duplication by handling the solver setup, + * directory creation, and output callback for any given scheme. + * + * @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"). */ void run_simulation(const PhysParams& phys, const NumParams& num, @@ -37,6 +42,7 @@ void run_simulation(const PhysParams& phys, HeatSolver solver(phys, num, scheme); // --- Prepare output directory --- + // Results will be saved in: results// fs::path outDir = fs::path("results") / schemeName; if (!fs::exists(outDir)) { @@ -44,14 +50,15 @@ void run_simulation(const PhysParams& phys, } // --- Output callback --- + // Called by the solver at t=0 and every 'outEvery' hours. auto onDump = [&](double t, const Grid& g, const std::vector& T) { std::ostringstream name; + // Filename format: profile_t_0.10h.csv 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 --- @@ -67,23 +74,25 @@ int main() // --- Physical parameters --- PhysParams phys; - phys.L_cm = 31.0; - phys.Tin = 38.0; - phys.Tsur = 149.0; + 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 positive and <= 1.0 cm."); - num.dt = 0.01; - num.tEnd = 0.5; - num.outEvery = 0.1; + + num.dt = 0.01; // Time step [h] + num.tEnd = 0.5; // Total duration [h] + num.outEvery = 0.1; // Output interval [h] - // --- Grid validation --- + // --- Grid validation (Early check) --- Grid g(phys.L_cm, num.dx); validate_grid(g); std::cout << "Grid created ✅ Nx = " << g.Nx << ", dx = " << g.dx << " cm\n"; @@ -92,23 +101,35 @@ int main() std::cout << "\nSelect numerical method:\n"; std::cout << "1. DuFort-Frankel\n"; std::cout << "2. Richardson\n"; - std::cout << "3. Run All (Batch)\n"; + std::cout << "3. Laasonen\n"; + std::cout << "4. Run All (Batch)\n"; + // We use askForDouble (or askForInteger) to get a valid choice between 1 and 4 int choice = static_cast(askForDouble( - "Enter choice (1-3): ", - [](double v) { return v >= 1.0 && v <= 3.0; }, - "Please enter 1, 2, or 3.")); + "Enter choice (1-4): ", + [](double v) { return v >= 1.0 && v <= 4.0; }, + "Please enter a number between 1 and 4.")); + + // --- Execution Logic --- - if (choice == 1 || choice == 3) + // 1. DuFort-Frankel + if (choice == 1 || choice == 4) { run_simulation(phys, num, SchemeKind::DuFortFrankel, "DuFort-Frankel"); } - if (choice == 2 || choice == 3) + // 2. Richardson + if (choice == 2 || choice == 4) { run_simulation(phys, num, SchemeKind::Richardson, "Richardson"); } + // 3. Laasonen + if (choice == 3 || choice == 4) + { + run_simulation(phys, num, SchemeKind::Laasonen, "Laasonen"); + } + std::cout << "\nAll requested simulations completed successfully ✅\n"; return 0; } @@ -117,4 +138,4 @@ int main() std::cerr << "Fatal error: " << e.what() << "\n"; return 1; } -} +} \ No newline at end of file diff --git a/src/method.cpp b/src/method.cpp index fa37a27..59d6020 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -1,26 +1,109 @@ -#include "method.hpp" +/** + * @file method.cpp + * @brief Factory implementation for numerical schemes. + */ +#include "method.hpp" #include "methods/dufort_frankel.hpp" +#include "methods/laasonen.hpp" #include "methods/richardson.hpp" + #include #include +// ============================================================================= +// ADAPTER CLASSES (Wrappers) +// These classes adapt the specific scheme interfaces to the generic Method interface. +// ============================================================================= + +/** + * @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*/, + const std::vector& Tcurr, + std::vector& Tnext) const override + { + 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_; +}; + +// ============================================================================= +// 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. - * @throws std::invalid_argument if the scheme is not available in this build. + * @throws std::invalid_argument if the scheme is not supported. */ std::unique_ptr make_method(SchemeKind scheme) { switch (scheme) { + case SchemeKind::Laasonen: + 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(); + default: - throw std::invalid_argument("make_method: unsupported scheme"); + throw std::invalid_argument("make_method: unsupported or unknown scheme"); } -} +} \ No newline at end of file diff --git a/src/methods/laasonen.cpp b/src/methods/laasonen.cpp new file mode 100644 index 0000000..8022c5d --- /dev/null +++ b/src/methods/laasonen.cpp @@ -0,0 +1,102 @@ +#include "methods/laasonen.hpp" + +#include +#include + +void Laasonen::step(const Grid& g, + double D, + double dt, + 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); + + // Sécurité + 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) = T^n + for (std::size_t j = 0; j < M; ++j) + { + d[j] = Tcurr[j + 1]; + } + + // Ajout des conditions aux limites dans le RHS + // 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 + // 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 + // 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 (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) + { + Tnext[j + 1] = x[j]; + } +} diff --git a/src/solver.cpp b/src/solver.cpp index 472af48..13f6ff9 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -4,7 +4,6 @@ */ #include "solver.hpp" - #include // std::fill HeatSolver::HeatSolver(const PhysParams& phys, const NumParams& num, SchemeKind scheme) @@ -35,7 +34,7 @@ void HeatSolver::apply_dirichlet(std::vector& T) const { // Enforce surface temperature at both ends T.front() = phys_.Tsur; - T.back() = phys_.Tsur; + T.back() = phys_.Tsur; } void HeatSolver::init_field() @@ -52,18 +51,18 @@ void HeatSolver::bootstrap_if_needed() } // We have T_ = T^0 here - - // 1) Store T^0 in Tprev_ + // Tprev_ will store T^0 Tprev_ = T_; - // 2) Construct T^1 with a small FTCS step (1D diffusion) + // Construct T^1 with a small FTCS step (1D diffusion) const double dx = grid_.dx; - const double r = phys_.D_cm2h * num_.dt / (dx * dx); + const double r = phys_.D_cm2h * num_.dt / (dx * dx); const std::size_t N = grid_.Nx; next_ = T_; // FTCS scheme (conditional stability: r <= 0.5) + // If r > 0.5, we just clone T0 to T1 to avoid oscillations at the first step if (r <= 0.5) { for (std::size_t i = 1; i + 1 < N; ++i) @@ -73,12 +72,11 @@ void HeatSolver::bootstrap_if_needed() } else { - // If dt is too large, keep T^1 = T^0 to avoid oscillations next_ = T_; } apply_dirichlet(next_); - T_ = next_; // T^1 + T_ = next_; // Now T_ holds T^1, Tprev_ holds T^0 } void HeatSolver::run( @@ -93,6 +91,10 @@ void HeatSolver::run( // Initial output onDump(t, grid_, T_); + // IMPORTANT FIX: Ensure 'next_' has BCs before calculating the first step. + // This prevents implicit schemes from seeing 0°C at boundaries instead of Tsur. + apply_dirichlet(next_); + // Main time loop while (t < tEnd - 1e-12) { @@ -103,11 +105,13 @@ void HeatSolver::run( } else { - // Two-level method (e.g. FTCS, BTCS...) + // Two-level method (e.g. Laasonen, Crank-Nicolson) + // Note: For 2-level, the third argument (Tprev) is ignored by the method wrapper. + // We pass T_ (current) as Tprev just to satisfy the signature. method_->step(grid_, phys_.D_cm2h, dt, T_, T_, next_); } - // Enforce boundary conditions + // Enforce boundary conditions (overwrite boundaries just in case) apply_dirichlet(next_); // Shift time states @@ -127,4 +131,4 @@ void HeatSolver::run( nextDump += outEvery; } } -} +} \ No newline at end of file