PlasmaColumnSolver.cpp

Go to the documentation of this file.

#include "PlasmaColumnSolver.h"
#include "ThermalPlasmaColumn1D.h"
#include "units.h"
 
#include "cantera/oneD/Sim1D.h"
#include "cantera/oneD/Boundary1D.h"
#include "cantera/oneD/Domain1D.h"
 
#include <memory>
 
namespace rizer {
 
ColumnResult solveColumn(double R, double electric_field, double T_wall,
                          const PropertyTable& sigma, const PropertyTable& kappa,
                          const PropertyTable& p_rad, const ColumnOptions& opts)
{
    using std::make_shared;
 
    // Optional initial guess T(r), folded into the domain so it is applied at
    // construction (via initialValue) without the deprecated Sim1D::setProfile.
    // Requires at least two points and matching lengths for init_r/init_T;
    // otherwise init_profile is left default-constructed and the domain falls
    // back to its own parabolic seed built from T_wall and T_center_guess.
    PropertyTable init_profile;
    if (opts.init_r.size() >= 2 && opts.init_r.size() == opts.init_T.size()) {
        init_profile = PropertyTable(opts.init_r, opts.init_T);
    }
 
    // The interior domain owns the actual Elenbaas-Heller physics: the radial
    // conduction/Joule-heating/radiation balance is assembled as residuals on
    // this domain's grid points; sigma/kappa/p_rad are consulted at each
    // Newton iterate via their tabulated PropertyTable lookups.
    auto column = make_shared<ThermalPlasmaColumn1D>(
        R, opts.n_points, electric_field, T_wall,
        sigma, kappa, p_rad, opts.T_center_guess, opts.rho_cp, init_profile);
 
    // Dummy no-op terminators (no Solution required), matching the pattern of
    // Cantera's BoundaryValueProblem sample. The column owns its own BCs.
    auto left = make_shared<Cantera::Empty1D>();
    auto right = make_shared<Cantera::Empty1D>();
 
    // Assemble the three-domain Sim1D stack [Empty1D | ThermalPlasmaColumn1D |
    // Empty1D]. Sim1D indexes domains in this vector order, so the column is
    // always domain index 1 regardless of how many points it has.
    std::vector<std::shared_ptr<Cantera::Domain1D>> domains{left, column, right};
    Cantera::Sim1D sim(domains);
 
    // Refinement applies to the column (domain index 1). These thresholds
    // (ratio/slope/curve/prune) come straight from ColumnOptions and control
    // how aggressively Sim1D inserts/removes grid points between Newton
    // solves when refine_grid is enabled below.
    sim.setRefineCriteria(1, opts.refine_ratio, opts.refine_slope,
                          opts.refine_curve, opts.refine_prune);
 
    // Drive the steady solve. Internally, Sim1D::solve() first attempts a
    // direct Newton solve on the current grid; if Newton fails to converge it
    // automatically falls back to pseudo-transient time-stepping (using
    // rho_cp baked into the column's residuals as the effective transient
    // time constant) to move the iterate closer to the solution before
    // retrying Newton. When refine_grid is true, this Newton/time-stepping
    // cycle is repeated on progressively refined grids (per the criteria set
    // above) until no further refinement is needed.
    sim.solve(opts.loglevel, opts.refine_grid);
 
    // Extract the converged profile from the (possibly refined) grid. Note
    // nPoints()/z(j) reflect the final refined grid, which may differ in size
    // and spacing from the initial opts.n_points uniform grid.
    const std::size_t np = column->nPoints();
    ColumnResult res;
    res.n_points = np;
    res.electric_field = electric_field;
    res.r.reserve(np);
    res.T.reserve(np);
    res.sigma.reserve(np);
    res.kappa.reserve(np);
 
    // Re-evaluate sigma(T)/kappa(T) at the converged temperatures (rather than
    // caching values from the last Newton iterate) so the returned arrays are
    // exactly consistent with the final T profile.
    std::vector<double> Tvals = column->values("T");
    for (std::size_t j = 0; j < np; j++) {
        double r = column->z(j);
        double T = Tvals[j];
        res.r.push_back(r);
        res.T.push_back(T);
        res.sigma.push_back(column->sigmaAt(T));
        res.kappa.push_back(column->kappaAt(T));
    }
 
    // Integrated current:
    // @f$I = \int_0^R \sigma(T)\, E\, 2\pi r\, dr@f$ (trapezoidal).
    double I = 0.0;
    for (std::size_t j = 0; j + 1 < np; j++) {
        // Trapezoidal rule:
        // @f$\int_a^b f(r)\, dr \approx \tfrac{1}{2}(f(a)+f(b))(b-a)@f$
        double f0 = res.sigma[j] * res.r[j];
        double f1 = res.sigma[j + 1] * res.r[j + 1];
        double dr = res.r[j + 1] - res.r[j];
        I += 0.5 * (f0 + f1) * dr;
    }
    // E is uniform along the column (fixed axial field), so it can be pulled
    // out of the radial integral: @f$I = 2\pi E \int \sigma(T)\, r\, dr@f$.
    res.current = 2.0 * units::pi * electric_field * I;
 
    return res;
}
 
} // namespace rizer