solver.cpp

Go to the documentation of this file.

/**
 * @file solver.cpp
 * @brief Time-marching driver for the 1D heat equation.
 */
 
#include "solver.hpp"
 
#include <algorithm> // std::fill
 
HeatSolver::HeatSolver(const PhysParams& phys, const NumParams& num, SchemeKind scheme)
    : phys_(phys)
    , num_(num)
    , grid_(phys.L_cm, num.dx)
{
    // Create the calculation method (factory based on scheme)
    method_ = make_method(scheme);
    uses_prev_step_ = method_->uses_previous_step();
 
    // Vector allocation
    T_.resize(grid_.Nx);
    next_.resize(grid_.Nx);
    if (uses_prev_step_)
    {
        Tprev_.resize(grid_.Nx);
    }
 
    // Field initialization
    init_field();
 
    // Bootstrap for 3-level methods (like DuFort-Frankel)
    bootstrap_if_needed();
}
 
void HeatSolver::apply_dirichlet(std::vector<double>& T) const
{
    // Enforce surface temperature at both ends
    T.front() = phys_.Tsur;
    T.back() = phys_.Tsur;
}
 
void HeatSolver::init_field()
{
    std::fill(T_.begin(), T_.end(), phys_.Tin);
    apply_dirichlet(T_);
}
 
void HeatSolver::bootstrap_if_needed()
{
    if (!uses_prev_step_)
    {
        return;
    }
 
    // We have T_ = T^0 here
    // Tprev_ will store T^0
    Tprev_ = T_;
 
    // 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 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)
        {
            next_[i] = T_[i] + r * (T_[i + 1] - 2.0 * T_[i] + T_[i - 1]);
        }
    }
    else
    {
        next_ = T_;
    }
 
    apply_dirichlet(next_);
    T_ = next_; // Now T_ holds T^1, Tprev_ holds T^0
}
 
void HeatSolver::run(
        const std::function<void(double, const Grid&, const std::vector<double>&)>& onDump)
{
    double t = 0.0;
    const double dt = num_.dt;
    const double tEnd = num_.tEnd;
    const double outEvery = num_.outEvery;
    double nextDump = 0.0;
 
    // -------------------------------------------------------------------------
    // 1. Initial State 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_);
 
    // -------------------------------------------------------------------------
    // 2. Main Time Stepping Loop
    // -------------------------------------------------------------------------
    while (t < tEnd - 1e-12)
    {
        // --- A. Compute Next Time Step (T^{n+1}) ---
        if (uses_prev_step_)
        {
            // Three-level method (e.g. DuFort-Frankel, Richardson)
            method_->step(grid_, phys_.D_cm2h, dt, Tprev_, T_, next_);
        }
        else
        {
            // 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_);
        }
 
        // --- B. Enforce Boundary Conditions ---
        // Overwrite boundaries just in case the method touched them (though it shouldn't)
        apply_dirichlet(next_);
 
        // --- C. Shift Time States (Rotation) ---
        if (uses_prev_step_)
        {
            Tprev_ = T_;
        }
        T_ = next_;
 
        // --- D. Advance Time ---
        t += dt;
 
        // --- E. Periodic Output ---
        if (t + 1e-12 >= nextDump)
        {
            onDump(t, grid_, T_);
            nextDump += outEvery;
        }
    }
}