laasonen.cpp

Go to the documentation of this file.

#include "methods/laasonen.hpp"
 
#include <cstddef>
#include <vector>
 
void Laasonen::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;
    const double r = D * dt / (dx * dx);
 
    // Safety check
    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; // Number of internal nodes (excluding boundaries)
 
    // Tridiagonal coefficients
    // Equation: -r T_{i-1}^{n+1} + (1 + 2r) T_i^{n+1} - r T_{i+1}^{n+1} = T_i^n
    std::vector<double> a(M, -r);
    std::vector<double> b(M, 1.0 + 2.0 * r);
    std::vector<double> c(M, -r);
    std::vector<double> d(M, 0.0);
 
    // RHS construction (Right-Hand Side) = T^n
    for (std::size_t j = 0; j < M; ++j)
    {
        d[j] = Tcurr[j + 1];
    }
 
    // Adding boundary conditions to 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}
 
    // -----------------------
    // Thomas Algorithm
    // -----------------------
 
    // 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<double> x(M);
    x[M - 1] = d[M - 1];
 
    for (int i = static_cast<int>(M) - 2; i >= 0; --i)
    {
        x[i] = d[i] - c[i] * x[i + 1];
    }
 
    // Update T^{n+1} (internal nodes)
    for (std::size_t j = 0; j < M; ++j)
    {
        Tnext[j + 1] = x[j];
    }
}