ThermalPlasmaColumn1D.cpp

Go to the documentation of this file.

#include "ThermalPlasmaColumn1D.h"
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <string>
 
namespace rizer {
 
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(0, 0.0, 1.0e5);             // physical temperature bounds [K]
    setSteadyTolerances(1.0e-4, 1.0e-9);  // relative, absolute
    setTransientTolerances(1.0e-4, 1.0e-11);
}
 
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 r = R.
    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);
    }
}
 
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).  In the
    // latter case we only need to update the three points that the stencil of
    // node jg touches; skip everything else to avoid redundant work.
    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) ------------------------------------
        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 (1/r) d/dr(r * kappa * dT/dr).
        //  The east face is at r_e = 0.5*(r_j + r_{j+1}), with spacing dr_e = r_{j+1} - r_j.
        //  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 : r_e * kap_e * dT/dr|_e = r_e * kap_e * (T(j+1) - Tj) / dr_e   (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 = re * kap_e * (Tjp1 - Tj) / dre;
 
        // ---- Conductive flux on the west face (between j-1 and j) -------
        //  At the axis (j=0) the face radius r_w = 0, so the flux vanishes
        //  naturally — this is the symmetry BC (dT/dr = 0 at r = 0).
        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 = rw * kap_w * (Tj - Tjm1) / drw;
        }
 
        // ---- Finite-volume divergence ------------------------------------
        //  vol = integral of r dr over the cell = 0.5*(r_e^2 - r_w^2)
        //  div_flux approximates (1/r) d/dr(r * kappa * dT/dr) at node j.
        const double vol      = 0.5 * (re * re - rw * rw);
        const double div_flux = (flux_e - flux_w) / vol;
 
        // ---- Source term: Joule heating minus radiative loss -------------
        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;
    }
}
 
} // namespace rizer