richardson.hpp

Go to the documentation of this file.

#pragma once
/**
 * @file richardson.hpp
 * @brief Richardson explicit (central time - central space) scheme for the 1D heat equation.
 */
 
#include "../grid.hpp"
#include "../method.hpp"
#include <string>
#include <vector>
 
/**
 * @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<double>& Tprev,
              const std::vector<double>& Tcurr,
              std::vector<double>& Tnext) const override;
};