ThermalPlasmaColumn1D.cpp

Go to the documentation of this file.

#include "ThermalPlasmaColumn1D.h"
#include "numeric/DiffusiveFlux.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>
 
namespace rizer {
 
// ---------------------------------------------------------------------------------- //
//   Constructor                                                                      //
// ---------------------------------------------------------------------------------- //
 
ThermalPlasmaColumn1D::ThermalPlasmaColumn1D(double R, std::size_t npoints,
                                             double electric_field,
                                             double T_wall, PropertyTable sigma,
                                             PropertyTable kappa,
                                             PropertyTable p_rad,
                                             double T_center_guess, double rho_cp,
                                             PropertyTable init_profile)
    : Cantera::Domain1D(/*nv=*/1, /*points=*/npoints) // 1 solution component: T(r)
    , m_R(R)
    , m_E(electric_field)
    , m_Twall(T_wall)
    , m_Tcenter(T_center_guess)
    , m_rho_cp(rho_cp)
    , m_sigma(std::move(sigma))
    , m_kappa(std::move(kappa))
    , m_prad(std::move(p_rad))
    , m_init(std::move(init_profile))
{
    if (R <= 0.0 || npoints < 3) {
        throw std::invalid_argument(
            "ThermalPlasmaColumn1D: require R > 0 and at least 3 grid points "
            "(got R=" + std::to_string(R) + ", n=" + std::to_string(npoints)
            + ").");
    }
    // Uniform initial radial grid r in [0, R]; PlasmaColumnSolver will call
    // sim.solve() which refines it adaptively via the Refiner criteria.
    setupUniformGrid(npoints, R, 0.0);
 
    setComponentName(0, "T");
    setBounds(/*n=*/0, /*lower=*/0.0, /*upper=*/1.0e5);  // [K]
    setSteadyTolerances(/*rtol=*/1.0e-4, /*atol=*/1.0e-9);
    setTransientTolerances(/*rtol=*/1.0e-4, /*atol=*/1.0e-11);
}
 
// ---------------------------------------------------------------------------------- //
//   Public methods                                                                   //
// ---------------------------------------------------------------------------------- //
 
void ThermalPlasmaColumn1D::eval(std::size_t jg, Cantera::span<const double> xg,
                                 Cantera::span<double> rg, Cantera::span<int> maskg,
                                 double rdt)
{
    // Cantera calls eval() both for a full residual sweep (jg == npos) and for
    // Jacobian columns (jg = global index of the perturbed variable), once per
    // column during finite-difference Jacobian assembly. In the latter case,
    // perturbing node jg only changes the residuals at jg itself and its two
    // immediate neighbours (the 3-point conduction stencil below), so a
    // cheap early-out here skips domains/points the perturbation cannot
    // possibly affect, saving redundant PropertyTable evaluations.
    if (jg != Cantera::npos
        && (jg + 1 < firstPoint() || jg > lastPoint() + 1)) {
        return;
    }
 
    // Slice the global arrays down to this domain's local range.
    auto x    = xg.subspan(loc(), size());    // local solution vector
                                               // (T at each grid point)
    auto rsd  = rg.subspan(loc(), size());    // local residual vector
                                               // (dT/dt at each grid point)
    auto diag = maskg.subspan(loc(), size()); // local mask vector (1 if
                                               // dT/dt exists, 0 if algebraic)
 
    // Determine the range of local indices to update.
    std::size_t jmin, jmax;
    if (jg == Cantera::npos) {             // full sweep
        jmin = 0;
        jmax = m_points - 1;
    } else {                               // single Jacobian column — 3-point stencil
        std::size_t jpt = (jg == 0) ? 0 : jg - firstPoint();
        jmin = std::max<std::size_t>(jpt, 1) - 1;
        jmax = std::min(jpt + 1, m_points - 1);
    }
 
    const double E2 = m_E * m_E;
 
    for (std::size_t j = jmin; j <= jmax; j++) {
        const double Tj = x[index(0, j)];
 
        // ---- Dirichlet wall BC (algebraic row, diag=0 tells the solver it
        //      has no time derivative) ------------------------------------
        //  The residual is simply @f$ T(r=R) - T_\text{wall} @f$, which
        //  Newton drives to zero by forcing the last node's temperature to
        //  the prescribed wall value; there is no flux balance at this row
        //  (the last interior cell's east flux uses this Dirichlet value as
        //  its neighbour, see flux_e above for j = m_points - 2).
        if (j == m_points - 1) {
            rsd[index(0, j)] = Tj - m_Twall;
            diag[index(0, j)] = 0;
            continue;
        }
 
        // ---- Conductive flux on the east face (between j and j+1) -------
        //  The finite-volume divergence is
        //  @f$ \frac{1}{r}\frac{d}{dr}\!\left(r \kappa \frac{dT}{dr}\right) @f$.
        //  The east face is at @f$ r_e = 0.5(r_j + r_{j+1}) @f$, with
        //  spacing @f$ dr_e = r_{j+1} - r_j @f$.
        //  The thermal conductivity at the face is the arithmetic average
        //  of the two nodes.
        //
        //  zj     : Radius of node j (i.e. z(0)=0, z(npts-1)=R)
        //  r_e    : Mid-point radius of the east face
        //           (i.e. r_e(0)=0.5*(z(0)+z(1)=0.5*z(1)))
        //  dr_e   : Node spacing to the east
        //  Tjp1   : Temperature at node j+1
        //  kap_e  : Arithmetic-average thermal conductivity at the face
        //  flux_e : @f$ r_e \kappa_e \left.\frac{dT}{dr}\right|_e
        //           = r_e \kappa_e (T_{j+1} - T_j) / dr_e @f$ (units: W/m)
        const double zj     = z(j);
        const double re     = 0.5 * (zj + z(j + 1));
        const double dre    = z(j + 1) - zj;
        const double Tjp1   = x[index(0, j + 1)];
        const double kap_e  = 0.5 * (m_kappa.eval(Tj) + m_kappa.eval(Tjp1));
        const double flux_e = numeric::faceFlux(re, kap_e, Tj, Tjp1, dre);
 
        // ---- Conductive flux on the west face (between j-1 and j) -------
        //  At the axis (j=0) there is no j-1 neighbour; the r=0 symmetry BC
        //  (@f$ dT/dr = 0 @f$ at @f$ r = 0 @f$, i.e. no heat flows "through"
        //  the centreline) is imposed geometrically rather than via an
        //  explicit equation: the west face radius r_w is 0 there, and
        //  since flux ~ r_w * ..., the west flux is forced to exactly zero
        //  regardless of dT/dr.
        //  No ghost node or extra BC row is needed for this case.
        double flux_w, rw;
        if (j == 0) {
            rw     = 0.0;
            flux_w = 0.0;
        } else {
            rw = 0.5 * (z(j - 1) + zj);
            const double Tjm1 = x[index(0, j - 1)];
            const double drw   = zj - z(j - 1);
            const double kap_w = 0.5 * (m_kappa.eval(Tjm1) + m_kappa.eval(Tj));
            flux_w = numeric::faceFlux(rw, kap_w, Tjm1, Tj, drw);
        }
 
        // ---- Finite-volume divergence ------------------------------------
        //  @f$ \text{vol} = \int r\,dr @f$ over the cell
        //  @f$ = 0.5(r_e^2 - r_w^2) @f$.
        //  div_flux approximates
        //  @f$ \frac{1}{r}\frac{d}{dr}\!\left(r \kappa \frac{dT}{dr}\right) @f$
        //  at node j.
        const double vol      = numeric::annularVolume(re, rw);
        const double div_flux = (flux_e - flux_w) / vol;
 
        // ---- Source term: Joule heating minus radiative loss -------------
        //  Unlike kappa (looked up at the two face-adjacent nodes and
        //  averaged, since it multiplies a flux between nodes), sigma and
        //  P_rad are cell-centred quantities and are looked up once at this
        //  node's own temperature Tj via PropertyTable::eval() (linear
        //  interpolation on the tabulated grid; m_prad.eval() returns 0 if
        //  no radiation table was supplied).
        const double source = m_sigma.eval(Tj) * E2 - m_prad.eval(Tj);
 
        // ---- Residual (scaled by rho*cp so Newton updates are in K) ------
        //  The pseudo-transient term rdt*(T - T_prev) adds a diagonal shift
        //  during Cantera's time-stepping fallback and vanishes at steady
        //  state (rdt -> 0 once Newton converges).
        rsd[index(0, j)] = (div_flux + source) / m_rho_cp
                           - rdt * (Tj - prevSoln(0, j));
        diag[index(0, j)] = 1;
    }
}
 
 
std::size_t ThermalPlasmaColumn1D::componentIndex(const std::string& name,
                                                   bool /*checkAlias*/) const
{
    if (name == "T") {
        return 0;
    }
    throw std::invalid_argument(
        "ThermalPlasmaColumn1D: no component named '" + name + "'.");
}
 
 
bool ThermalPlasmaColumn1D::hasComponent(const std::string& name,
                                          bool /*checkAlias*/) const
{
    return name == "T";
}
 
 
void ThermalPlasmaColumn1D::getValues(const std::string& component,
                                       Cantera::span<double> values) const
{
    if (component != "T") {
        throw std::invalid_argument(
            "ThermalPlasmaColumn1D: no component named '" + component + "'.");
    }
    if (!m_state) {
        throw std::runtime_error(
            "ThermalPlasmaColumn1D::getValues: domain not installed in a container.");
    }
    const double* soln = m_state->data() + loc();
    for (std::size_t j = 0; j < m_points; j++) {
        values[j] = soln[index(0, j)];
    }
}
 
 
double ThermalPlasmaColumn1D::initialValue(std::size_t /*n*/, std::size_t j)
{
    double r = z(j);
    if (!m_init.empty()) {
        // Seeded guess (e.g. the analytical EH profile), interpolated in r.
        return m_init.eval(r);
    }
    // Default parabola: T_center at the axis, T_wall at @f$ r = R @f$.
    double s = r / m_R;
    return m_Twall + (m_Tcenter - m_Twall) * (1.0 - s * s);
}
 
 
void ThermalPlasmaColumn1D::resetBadValues(Cantera::span<double> xg)
{
    auto x = xg.subspan(loc(), size());
    double lo = lowerBound(0);
    double hi = upperBound(0);
    for (std::size_t j = 0; j < m_points; j++) {
        double T = x[index(0, j)];
        x[index(0, j)] = std::min(std::max(T, lo), hi);
    }
}
 
} // namespace rizer