PlasmaColumnSolver.cpp

Go to the documentation of this file.

#include "PlasmaColumnSolver.h"
#include "ThermalPlasmaColumn1D.h"
 
#include "cantera/oneD/Sim1D.h"
#include "cantera/oneD/Boundary1D.h"
#include "cantera/oneD/Domain1D.h"
 
#include <memory>
 
namespace rizer {
 
namespace { constexpr double kPi = 3.14159265358979323846; }
 
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.
    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);
    }
 
    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>();
 
    std::vector<std::shared_ptr<Cantera::Domain1D>> domains{left, column, right};
    Cantera::Sim1D sim(domains);
 
    // Refinement applies to the column (domain index 1).
    sim.setRefineCriteria(1, opts.refine_ratio, opts.refine_slope,
                          opts.refine_curve, opts.refine_prune);
 
    sim.solve(opts.loglevel, opts.refine_grid);
 
    // Extract the converged profile from the (possibly refined) 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);
 
    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: I = int_0^R sigma(T) * E * 2*pi*r dr  (trapezoidal).
    double I = 0.0;
    for (std::size_t j = 0; j + 1 < np; j++) {
        // Trapezoidal rule:  int_a^b f(r) dr ~ 0.5*(f(a)+f(b))*(b-a)
        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;
    }
    res.current = 2.0 * kPi * electric_field * I;
 
    return res;
}
 
} // namespace rizer