PlasmaColumnSolver.cpp

Go to the documentation of this file.

#include "PlasmaColumnSolver.h"
#include "ThermalPlasmaColumn1D.h"
#include "units.h"
#include "numeric/Format.h"
 
#include "cantera/oneD/Sim1D.h"
#include "cantera/oneD/Boundary1D.h"
#include "cantera/oneD/Domain1D.h"
 
#include <algorithm>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <string>
 
namespace rizer {
 
// ---------------------------------------------------------------------------------- //
//   Free functions                                                                   //
// ---------------------------------------------------------------------------------- //
 
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.
    // Both empty leaves init_profile default-constructed, and the domain builds
    // its own parabolic seed from T_wall and T_center_guess. Any other pair that
    // is not >= 2 matching points is refused: falling back to the parabolic seed
    // would discard the caller's guess without a word.
    PropertyTable init_profile;
    if (!opts.init_r.empty() || !opts.init_T.empty()) {
        if (opts.init_r.size() < 2 || opts.init_r.size() != opts.init_T.size()) {
            throw std::invalid_argument(
                "solveColumn: init_r and init_T must be given together with the "
                "same length of at least 2 points, or both omitted for the "
                "parabolic seed; got " + std::to_string(opts.init_r.size()) +
                " and " + std::to_string(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(/*dom=*/1, /*ratio=*/opts.refine_ratio,
                          /*slope=*/opts.refine_slope,
                          /*curve=*/opts.refine_curve,
                          /*prune=*/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(/*loglevel=*/opts.loglevel, /*refine_grid=*/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");
    // Newton iterates may leave the tables mid-solve (resetBadValues brings them
    // back), so the tables are checked against the converged column instead.
    // The tolerance absorbs the round-off of a wall held at a table's first
    // point.
    auto requireCovered = [&](const PropertyTable &table, const char *name)
    {
        for (double T : Tvals)
        {
            const double tolerance = 1.0e-6 * std::max(1.0, std::abs(T));
            if (!table.empty() &&
                !(T >= table.first() - tolerance && T <= table.last() + tolerance))
            {
                throw std::runtime_error(
                    std::string("solveColumn: the converged column reaches T = ") +
                    significant(T) + " K, outside the " + name + " table [" +
                    significant(table.first()) + ", " +
                    significant(table.last()) +
                    "] K, where the property is held at its end value.");
            }
        }
    };
    requireCovered(sigma, "sigma(T)");
    requireCovered(kappa, "kappa(T)");
    requireCovered(p_rad, "P_rad(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