crank_nicolson.cpp

Go to the documentation of this file.

/**
 * @file crank_nicolson.cpp
 * @brief Implementation of the Crank-Nicolson scheme using the Thomas Algorithm.
 */
 
#include "methods/crank_nicolson.hpp"
 
#include <cstddef>
#include <vector>
 
void CrankNicolson::step(const Grid& g,
                         double D,
                         double dt,
                         const std::vector<double>& /*Tprev*/,
                         const std::vector<double>& Tcurr,
                         std::vector<double>& Tnext) 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;
 
    // 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;
 
    // --- 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<double> a(M, -alpha);  // Lower diagonal
    std::vector<double> b(M, 1.0 + r); // Main diagonal
    std::vector<double> c(M, -alpha);  // Upper diagonal
    std::vector<double> 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 * 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 * 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)
    {
        double m = a[i] / b[i - 1];
        b[i] -= m * c[i - 1];
        d[i] -= m * d[i - 1];
    }
 
    // Backward substitution phase
    std::vector<double> 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)
    {
        Tnext[j + 1] = x[j];
    }
}