diff --git a/.gitignore b/.gitignore index e0a43c8..a448d57 100644 --- a/.gitignore +++ b/.gitignore @@ -15,7 +15,7 @@ .cache/ # Outputs -/results/* +/result/* !.gitkeep #MacOs Files diff --git a/include/method.hpp b/include/method.hpp index 80f27f0..42646f7 100644 --- a/include/method.hpp +++ b/include/method.hpp @@ -14,8 +14,9 @@ */ enum class SchemeKind { - DuFortFrankel - // Richardson, Laasonen, CrankNicolson (to be added later) + DuFortFrankel, + Richardson + // Laasonen, CrankNicolson (to be added later) }; /** diff --git a/include/methods/richardson.hpp b/include/methods/richardson.hpp new file mode 100644 index 0000000..80f23da --- /dev/null +++ b/include/methods/richardson.hpp @@ -0,0 +1,60 @@ +#pragma once +/** + * @file richardson.hpp + * @brief Richardson explicit (central time - central space) scheme for the 1D heat equation. + */ + +#include +#include +#include "../grid.hpp" +#include "../method.hpp" + +/** + * @class Richardson + * @brief Implementation of the explicit Richardson scheme (CTCS). + * + * The scheme reads: + * \f[ + * T_i^{n+1} = T_i^{n-1} + 2 r \left( T_{i+1}^n - 2 T_i^n + T_{i-1}^n \right), + * \quad r = \frac{D \, \Delta t}{\Delta x^2} + * \f] + * + * It is a three-level scheme, therefore it requires both \f$T^{n-1}\f$ and \f$T^n\f$ + * to produce \f$T^{n+1}\f$. Boundary conditions are not enforced here, but left + * to the solver. + */ +class Richardson : public Method { +public: + /** + * @brief Indicate that this is a three-level scheme. + * @return true always, since Richardson needs \f$T^{n-1}\f$. + */ + bool uses_previous_step() const noexcept override { return true; } + + /** + * @brief Get human-readable name of the scheme. + * @return "Richardson" + */ + std::string name() const override { return "Richardson"; } + + /** + * @brief Advance one time step using the Richardson explicit scheme. + * + * This computes \f$T^{n+1}\f$ from \f$T^{n}\f$ and \f$T^{n-1}\f$. + * Interior nodes \f$i = 1, \dots, N-2\f$ are updated, while boundary + * nodes are left untouched (the solver will impose Dirichlet BCs). + * + * @param g Grid metadata (spacing and number of nodes) + * @param D Diffusivity [cm^2/h] + * @param dt Time step [h] + * @param Tprev Temperature field at previous time layer \f$T^{n-1}\f$ + * @param Tcurr Temperature field at current time layer \f$T^{n}\f$ + * @param Tnext Output temperature field to be filled with \f$T^{n+1}\f$ + */ + void step(const Grid& g, + double D, + double dt, + const std::vector& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) override; +}; diff --git a/src/io.cpp b/src/io.cpp index 3b0dd9f..1732cef 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; @@ -23,9 +24,12 @@ void save_profile_csv(const std::string& path, if (x.size() != T.size()) throw std::runtime_error("save_profile_csv: vector size mismatch"); - // --- Ensure the result directory exists --- - fs::path fullPath = fs::path("result") / path; - fs::create_directories(fullPath.parent_path()); + // --- Ensure the directory exists --- + fs::path fullPath(path); + if (fullPath.has_parent_path()) + { + fs::create_directories(fullPath.parent_path()); + } // --- Open file for writing --- std::ofstream file(fullPath); diff --git a/src/main.cpp b/src/main.cpp index d398532..4987d7a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,9 +1,12 @@ /** * @file main.cpp - * @brief Run DuFort-Frankel simulation with user-chosen dx, dt; dump CSV every 0.1 h. + * @brief Test driver for the 1D heat conduction solver. + * Supports DuFort-Frankel and Richardson schemes. */ + #include "grid.hpp" #include "io.hpp" +#include "method.hpp" #include "solver.hpp" #include "types.hpp" #include "user_input.h" @@ -11,57 +14,107 @@ #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 + */ +void run_simulation(const PhysParams& phys, + const NumParams& num, + SchemeKind scheme, + const std::string& schemeName) +{ + std::cout << "\n=== Running " << schemeName << " Solver ===\n"; + + // --- Create solver --- + HeatSolver solver(phys, num, scheme); + + // --- Prepare output directory --- + fs::path outDir = fs::path("results") / schemeName; + if (!fs::exists(outDir)) + { + fs::create_directories(outDir); + } + + // --- Output callback --- + auto onDump = [&](double t, const Grid& g, const std::vector& T) + { + std::ostringstream name; + 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 --- + solver.run(onDump); + std::cout << "Results saved to " << outDir << "\n"; +} int main() { try { - PhysParams phys; // L=31 cm, D=93 cm^2/h, Tin=38°C, Tsur=149°C - NumParams num; // defaults: dx=0.05, dt=0.01, tEnd=0.5, outEvery=0.1 + std::cout << "=== 1D Heat Equation Solver ===\n"; - // Ask user for dx, dt (dt choices per assignment) - std::cout << "=== DuFort-Frankel simulation ===\n"; + // --- Physical parameters --- + PhysParams phys; + 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 in (0, 1]."); - // --- Run DuFort-Frankel Solver --- - std::cout << "=== Running DuFort-Frankel Solver ===\n"; - // Use default dt and tEnd from NumParams or ask user? - // For now use defaults: dt=0.01, tEnd=0.5 - - HeatSolver solver(phys, num, SchemeKind::DuFortFrankel); - - // --- Prepare output directory --- - namespace fs = std::filesystem; - std::string schemeName = - "DuFort-Frankel"; // Could get this from solver.method_name() if exposed, but - // hardcoding for now is fine or better yet, expose it. - // Actually, let's just hardcode "DuFort-Frankel" as requested or use the scheme enum. - // The user said: results/"method name"/... - - fs::path outDir = fs::path("results") / schemeName; - if (!fs::exists(outDir)) + "Δx must be positive and <= 1.0 cm."); + num.dt = 0.01; + num.tEnd = 0.5; + num.outEvery = 0.1; + + // --- Grid validation --- + Grid g(phys.L_cm, num.dx); + validate_grid(g); + std::cout << "Grid created ✅ Nx = " << g.Nx << ", dx = " << g.dx << " cm\n"; + + // --- Method Selection --- + std::cout << "\nSelect numerical method:\n"; + std::cout << "1. DuFort-Frankel\n"; + std::cout << "2. Richardson\n"; + std::cout << "3. Run All (Batch)\n"; + + int choice = static_cast(askForDouble( + "Enter choice (1-3): ", + [](double v) { return v >= 1.0 && v <= 3.0; }, + "Please enter 1, 2, or 3.")); + + if (choice == 1 || choice == 3) { - fs::create_directories(outDir); + run_simulation(phys, num, SchemeKind::DuFortFrankel, "DuFort-Frankel"); } - auto onDump = [&](double t, const Grid& g, const std::vector& T) + if (choice == 2 || choice == 3) { - std::ostringstream name; - 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); - std::cout << "Saved " << fullPath.string() << "\n"; - }; - - solver.run(onDump); - std::cout << "Simulation complete.\n"; + run_simulation(phys, num, SchemeKind::Richardson, "Richardson"); + } + + std::cout << "\nAll requested simulations completed successfully ✅\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 2028f61..fa37a27 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -1,16 +1,25 @@ #include "method.hpp" #include "methods/dufort_frankel.hpp" +#include "methods/richardson.hpp" #include #include +/** + * @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. + */ std::unique_ptr make_method(SchemeKind scheme) { switch (scheme) { - // @brief Factory implementation. case SchemeKind::DuFortFrankel: return std::make_unique(); + case SchemeKind::Richardson: + return std::make_unique(); default: throw std::invalid_argument("make_method: unsupported scheme"); } diff --git a/src/methods/richardson.cpp b/src/methods/richardson.cpp new file mode 100644 index 0000000..a79ab77 --- /dev/null +++ b/src/methods/richardson.cpp @@ -0,0 +1,25 @@ +/** + * @file richardson.cpp + * @brief Richardson explicit scheme implementation. + */ + +#include "methods/richardson.hpp" +#include // for std::size_t + +void Richardson::step(const Grid& g, + double D, + double dt, + const std::vector& Tprev, + const std::vector& Tcurr, + std::vector& Tnext) +{ + const double r = D * dt / (g.dx * g.dx); + const std::size_t N = g.Nx; + + // Update interior nodes only; boundaries are imposed by the solver. + for (std::size_t i = 1; i + 1 < N; ++i) { + Tnext[i] = Tprev[i] + + 2.0 * r * ( Tcurr[i + 1] - 2.0 * Tcurr[i] + Tcurr[i - 1] ); + } + // Tnext[0] and Tnext[N-1] left as-is +} diff --git a/src/solver.cpp b/src/solver.cpp index 3f92418..472af48 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -1,11 +1,16 @@ +/** + * @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 calculation method (factory based on scheme) method_ = make_method(scheme); @@ -35,10 +40,7 @@ void HeatSolver::apply_dirichlet(std::vector& T) const void HeatSolver::init_field() { - // Uniform initial field at Tin std::fill(T_.begin(), T_.end(), phys_.Tin); - - // Boundary conditions apply_dirichlet(T_); } @@ -59,7 +61,6 @@ void HeatSolver::bootstrap_if_needed() 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) @@ -76,11 +77,8 @@ void HeatSolver::bootstrap_if_needed() next_ = T_; } - // Boundary conditions apply_dirichlet(next_); - - // 3) Update T_ to contain T^1 - T_ = next_; + T_ = next_; // T^1 } void HeatSolver::run( @@ -90,10 +88,9 @@ void HeatSolver::run( const double dt = num_.dt; const double tEnd = num_.tEnd; const double outEvery = num_.outEvery; - double nextDump = 0.0; - // Initial dump + // Initial output onDump(t, grid_, T_); // Main time loop @@ -101,7 +98,7 @@ void HeatSolver::run( { if (uses_prev_step_) { - // Three-level method (e.g. DuFort-Frankel) + // Three-level method (e.g. DuFort-Frankel, Richardson) method_->step(grid_, phys_.D_cm2h, dt, Tprev_, T_, next_); } else