TransientPlasmaChannelSolver.cpp

Go to the documentation of this file.

#include "TransientPlasmaChannelSolver.h"
#include "PlasmaChannel1D.h"
 
#include "cantera/oneD/Sim1D.h"
#include "cantera/oneD/Boundary1D.h"
#include "cantera/oneD/MultiNewton.h"
#include "cantera/numerics/Integrator.h"
#include "cantera/numerics/FuncEval.h"
#include "cantera/base/ctexceptions.h"
 
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <sstream>
#include <functional>
#include <memory>
#include <stdexcept>
 
namespace rizer {
 
namespace {
 
// Method-of-lines RHS for the stiff (CVODE) backend:
// @f[
//     \dot{y} = f(y)
// @f]
// is the channel's spatial operator (transport + chemistry +
// Joule/elastic). Because PlasmaChannel1D::eval returns each node's
// @f$d\phi/dt@f$ (source/capacitance plus transport/capacitance) minus the
// transient term @f$r_{dt}(\phi-\phi_{prev})@f$, calling it with
// @f$r_{dt}=0@f$ yields exactly @f$\dot{y}@f$. The wall is a ghost-cell flux
// (no algebraic constraint), so the whole system is a clean ODE.
class ChannelRHS : public Cantera::FuncEval
{
public:
    ChannelRHS(PlasmaChannel1D* channel, Cantera::OneDim* sim, std::vector<double> y0,
               std::size_t n_vars, std::size_t channel_offset, std::size_t n_points,
               std::function<double(double)> fieldAt)
        : m_channel(channel), m_sim(sim), m_y0(std::move(y0)), m_n_vars(n_vars),
          m_channel_offset(channel_offset), m_n_points(n_points),
          m_field(std::move(fieldAt)) {}
 
    std::size_t neq() const override { return m_y0.size(); }

    void getState(Cantera::span<double> y) override {
        std::copy(m_y0.begin(), m_y0.end(), y.begin());
    }

    void eval(double t, Cantera::span<const double> y, Cantera::span<double> ydot,
              Cantera::span<const double> /*params*/) override {
        m_channel->setElectricField(m_field(t));
        // @f$r_{dt}=0@f$ -> the OneDim residual is exactly the
        // method-of-lines RHS @f$d\phi/dt@f$. Route through OneDim::eval
        // (the proven path that manages the mask workspace), not the
        // domain's eval directly.
        m_sim->eval(Cantera::npos, y, ydot, 0.0, 0);
    }

    void getConstraints(Cantera::span<double> constraints) override {
        std::fill(constraints.begin(), constraints.end(), 0.0);
        for (std::size_t j = 0; j < m_n_points; j++) {
            for (std::size_t comp = 2; comp < m_n_vars; comp++) {
                // 1.0 => y_i >= 0
                constraints[m_channel_offset + m_n_vars * j + comp] = 1.0;
            }
        }
    }
 
private:
    PlasmaChannel1D* m_channel;  // Non-owning; owned by solveChannelTransient's
                                  // `channel`
    Cantera::OneDim* m_sim;      // Non-owning; owned by solveChannelTransient's
                                  // `sim`
    std::vector<double> m_y0;    // Initial global state, units per-component
                                  // (Tg/Te [K], Y_k [-])
    std::size_t m_n_vars;          // Components per node, 2+K [-]
    std::size_t m_channel_offset;  // This domain's offset into the global
                                    // state [-]
    std::size_t m_n_points;        // Number of radial grid points [-]
    std::function<double(double)> m_field;  // E(t) [V/m], t in [s]
};
 
} // namespace
 
ChannelHistory solveChannelTransient(const ChannelOptions& opts)
{
    using std::make_shared;
 
    auto channel = make_shared<PlasmaChannel1D>(opts.cfg);
 
    auto left = make_shared<Cantera::Empty1D>();
    auto right = make_shared<Cantera::Empty1D>();
    std::vector<std::shared_ptr<Cantera::Domain1D>> domains{left, channel, right};
    Cantera::Sim1D sim(domains);
    sim.getInitialSoln();              // Apply PlasmaChannel1D::initialValue
 
    const std::size_t n_species = channel->nSpecies();
    const std::size_t n_points = channel->nPoints();
    const std::size_t n_vars = 2 + n_species;  // Per-node: Tg, Te, Y_1..Y_K
    const std::size_t channel_offset =
        channel->loc();  // This domain's offset in the global state
 
    std::vector<double> x(sim.size(), 0.0), xnew(sim.size());
 
    // Initialize x from the populated solution (getInitialSoln wrote it into
    // the internal m_state). NOTE: getState() copies the Newton work buffer
    // m_xnew (empty here), so read the real state via the domain's values()
    // (Domain1D::getValues reads m_state).
    for (std::size_t comp = 0; comp < n_vars; comp++) {
        std::vector<double> vals = channel->values(channel->componentName(comp));
        for (std::size_t j = 0; j < n_points; j++) {
            x[channel_offset + n_vars * j + comp] = vals[j];
        }
    }
 
    // Populate the transient mask before the first step. timeStep ->
    // initTimeInteg calls updateTransient()+factorize() BEFORE any residual
    // eval, so without this warm-up the mask is all-zero and the (zero) steady
    // Jacobian is singular on step 1. One eval(npos) fills m_mask (our diag=1).
    sim.eval(Cantera::npos, x, xnew, 1.0 / opts.dt, 0);
 
    // Static metadata (grid, species names, sizes) is fixed for the whole run
    // and only needs setting once, up front; only t/Tg/Te/ne/Y grow per frame.
    ChannelHistory history;
    history.npts = n_points;
    history.nsp = n_species;
    for (std::size_t j = 0; j < n_points; j++) history.r.push_back(channel->z(j));
    for (std::size_t k = 0; k < n_species; k++) {
        history.species.push_back(channel->componentName(2 + k));
    }
 
    // Appends one frame (all nodes) to `history` from the CURRENT flat state
    // `x`, in the flat layout x[channel_offset + n_vars*j + component] shared
    // by both backends below. Captured by reference so both the BDF and
    // Backward-Euler loops can call it identically at their own cadence.
    auto record = [&](double t) {
        history.t.push_back(t);
        for (std::size_t j = 0; j < n_points; j++) {
            const double Tg = x[channel_offset + n_vars * j + 0];
            const double Te = x[channel_offset + n_vars * j + 1];
            history.Tg.push_back(Tg);
            history.Te.push_back(Te);
            history.ne.push_back(
                channel->electronDensity(Tg, Te, &x[channel_offset + n_vars * j + 2]));
            for (std::size_t k = 0; k < n_species; k++) {
                history.Y.push_back(x[channel_offset + n_vars * j + 2 + k]);
            }
        }
    };
 
    // Field schedule E(t): an explicit table if provided, else the constant.
    PropertyTable Efield;
    if (opts.cfg.Et_t.size() >= 2
        && opts.cfg.Et_t.size() == opts.cfg.Et_v.size()) {
        Efield = PropertyTable(opts.cfg.Et_t, opts.cfg.Et_v);
    }
    auto fieldAt = [&](double t) {
        return Efield.empty() ? opts.cfg.electric_field : Efield.eval(t);
    };
 
    record(0.0);
 
    // Fixed output times, uniform over the total interval
    // @f$dt \cdot n_\text{steps}@f$.
    const double t_total =  // Total simulated time [s]
        opts.dt * static_cast<double>(opts.n_steps);
    const std::size_t n_intervals =  // Number of recorded frames after t=0 [-]
        std::max<std::size_t>(
            1, opts.n_steps / std::max<std::size_t>(1, opts.record_every));
 
    // ---------------- Stiff BDF backend (CVODE) ------------------------
    if (opts.integrator == "bdf") {
        ChannelRHS rhs(channel.get(), &sim, x, n_vars, channel_offset, n_points,
                       fieldAt);
        std::unique_ptr<Cantera::Integrator> cvode(Cantera::newIntegrator("CVODE"));
        cvode->setMethod(Cantera::BDF_Method);
        const char* lin_solver_env = std::getenv("RIZER_LINSOL");
        const std::string lin_solver_type =
            lin_solver_env ? std::string(lin_solver_env) : std::string("BAND");
        cvode->setLinearSolverType(lin_solver_type);
        if (lin_solver_type == "BAND") {
            // Each node couples to j +/- 1
            const int bandwidth = static_cast<int>(2 * n_vars);
            cvode->setBandwidth(bandwidth, bandwidth);
        }
        // Per-component absolute tolerances: K-scale for T, mass-fraction
        // floor for Y_k. atol[Y] must sit ABOVE the truly-trace species
        // level: a species at Y~1e-12 oscillating about zero cannot be
        // tracked to 1e-15, and CVODE then drives its internal step chasing
        // meaningless trace-species error (seen as a stall at the channel
        // edge, where chemistry is stiff but species are vanishing). 1e-12
        // ignores dynamically-irrelevant trace species; rtol 1e-7 governs the
        // O(1) carriers and the temperatures (this is a transport+chemistry
        // PDE solve, not a 0D point, so it does not need the 0D's rtol
        // 1e-10).
        std::vector<double> atol(x.size());
        for (std::size_t j = 0; j < n_points; j++) {
            atol[n_vars * j + 0] = 1.0e-3;              // Tg [K]
            atol[n_vars * j + 1] = 1.0e-3;              // Te [K]
            for (std::size_t k = 0; k < n_species; k++) {
                atol[n_vars * j + 2 + k] = 1.0e-9;      // Y_k [-]
            }
        }
        cvode->setTolerances(1.0e-7,
            Cantera::span<const double>(atol.data(), atol.size()));
        // Cap the internal step at the output cadence. This is essential for the
        // stiff plasma electron-energy equation (tiny electron heat capacity ->
        // sub-ns Te dynamics): without it CVODE attempts an over-large first step
        // and stalls. The 0D NRP reactor relies on the same max_step setting.
        cvode->setMaxStepSize(opts.dt);
        // A sane step budget: high enough for a long stiff run, low enough that a
        // genuine stall throws a clean (decoded) error instead of spinning at a
        // denormal step until SUNDIALS corrupts (the access-violation we saw).
        cvode->setMaxSteps(500000);
        cvode->initialize(0.0, rhs);
 
        // CVODE reports failures as bare flat-state indices ("275: 25.0 ..."),
        // which are meaningless to a caller. Translate each index into the node
        // (with its radius) and physical component (Tg / Te / species name) it
        // belongs to, so the message points at *where* and *what* went stiff.
        auto decodeError = [&](const std::string& message) -> std::string {
            // Failure time [s], if the message carries it
            // ("At t = <x> and h = ...").
            double t_fail = -1.0;
            // Character offset into `message` [-]
            std::size_t t_pos = message.find("At t = ");
            if (t_pos != std::string::npos) {
                try { t_fail = std::stod(message.substr(t_pos + 7)); } catch (...) {}
            }
            std::ostringstream out;
            out << "PlasmaChannel1D BDF solve failed to converge";
            if (t_fail >= 0.0) out << " near t=" << t_fail << " s";
            out << ".\nThe implicit corrector could not resolve the local "
                   "dynamics and the step size hit its floor. The dominant "
                   "error contributions are:\n";
            std::istringstream in(message);
            std::string line;
            bool in_error_block = false, found_any = false;
            while (std::getline(in, line)) {
                if (line.find("Components with largest") != std::string::npos) {
                    in_error_block = true;
                    continue;
                }
                if (!in_error_block) continue;
                std::size_t colon = line.find(':');
                if (colon == std::string::npos) continue;
                long global_index; double error_value;
                try {
                    global_index = std::stol(line.substr(0, colon));
                    error_value = std::stod(line.substr(colon + 1));
                } catch (...) { continue; }
                long local_index = global_index - static_cast<long>(channel_offset);
                if (local_index < 0) continue;
                std::size_t j = static_cast<std::size_t>(local_index) / n_vars;
                std::size_t comp = static_cast<std::size_t>(local_index) % n_vars;
                if (j >= n_points) continue;
                out << "  node " << j << " (r=" << channel->z(j) << " m)  "
                    << channel->componentName(comp)
                    << "  (weighted error " << error_value << ")\n";
                found_any = true;
            }
            if (!found_any) { out << message << "\n"; }
            out << "Remedies: reduce dt, loosen tolerances, soften the initial "
                   "profiles or the field ramp, or verify the initial state is "
                   "physically consistent at the flagged nodes.";
            return out.str();
        };
 
        // Graded startup. The 2T electron-energy + plasma-chemistry layer is
        // ferociously stiff at t=0 (tiny electron heat capacity ->
        // @f$dT_e/dt \sim 10^{15}@f$ K/s once the field is on). The proven 0D
        // reactor tames this with scipy vode's first_step=1e-15; Cantera's
        // CVODES exposes no initial-step setter, but CVODE caps its first
        // internal step at (tout - t0), so integrating to a sequence of tiny
        // times forces a small, safe first step that then grows under error
        // control -- the same effect, from the caller side.
        // Warm-up integration times [s]
        for (double t_warm : {1.0e-15, 1.0e-14, 1.0e-13, 1.0e-12, 1.0e-11, 1.0e-10}) {
            if (t_warm >= t_total) break;
            try {
                cvode->integrate(t_warm);
            } catch (const Cantera::CanteraError& err) {
                throw Cantera::CanteraError("solveChannelTransient",
                                            "{}", decodeError(err.what()));
            }
        }
 
        for (std::size_t frame = 1; frame <= n_intervals; frame++) {
            const double t_rec = t_total * static_cast<double>(frame)
                                           / static_cast<double>(n_intervals);
            std::string decoded;
            bool failed = false;
            try {
                cvode->integrate(t_rec);
            } catch (const Cantera::CanteraError& err) {
                decoded = decodeError(err.what());   // Format-free std::string
                failed = true;
            }
            if (failed) {
                // Re-throw OUTSIDE the catch handler, and pass the decoded text as
                // a format ARGUMENT ("{}") so CanteraError's fmt layer never parses
                // any stray braces in the CVODE message as format fields.
                throw Cantera::CanteraError("solveChannelTransient", "{}", decoded);
            }
            auto solution = cvode->solution();
            std::copy(solution.begin(), solution.end(), x.begin());
            record(t_rec);
            if (opts.loglevel > 0) {
                std::fprintf(stderr, "[bdf] frame %zu/%zu  t=%.3e s\n",
                             frame, n_intervals, t_rec);
                std::fflush(stderr);
            }
        }
        history.nframes = history.t.size();
        return history;
    }
 
    // ---------------- Adaptive Backward-Euler backend (MultiNewton) --------
    // Shrink dt on a failed Newton solve, grow it back on success. (A fixed step
    // -- m_tfactor=1 plus unbounded retries -- is what made stiff cases hang.
    // B-E is L-stable, so the constraint is the Newton basin, not stability;
    // letting dt adapt restores Cantera's own escape hatch.) Each sub-step
    // mirrors one SteadyStateSystem::timeStep iteration (initTimeInteg applies
    // the rdt transient term + factorizes, then a Newton solve) with explicit dt
    // and time control so frames land on the requested times.
    sim.newton().setOptions(/*maxJacAge=*/10);
    sim.setMinTimeStep(1.0e-16);
    sim.setMaxTimeStep(opts.dt);
    sim.setTimeStepFactor(0.5);
 
    const double dt_min = 1.0e-16;      // Smallest allowed adaptive step [s]
                                         // before giving up
    const long max_substeps = 2000000;  // Cumulative sub-step budget [-]
                                         // across the whole run
    auto& newton = sim.newton();
 
    double t = 0.0;            // Current simulated time [s]
    double dt = opts.dt;       // Current adaptive step size [s]
    int fails = 0;             // Consecutive Newton failures at the
                                // current dt [-]
    long total_substeps = 0;   // Cumulative sub-steps taken so far [-]
 
    for (std::size_t frame = 1; frame <= n_intervals; frame++) {
        const double t_rec =  // This frame's target recording time [s]
            t_total * static_cast<double>(frame) / static_cast<double>(n_intervals);
        while (t < t_rec - 1.0e-300) {
            // Clipped so frames land exactly on t_rec [s]
            const double step_size = std::min(dt, t_rec - t);
            channel->setElectricField(fieldAt(t));      // E(t) at the sub-step start
            sim.initTimeInteg(step_size, x);
            const int status = newton.solve(x, xnew, sim, opts.loglevel);
            if (status >= 0) {
                std::copy(xnew.begin(), xnew.end(), x.begin());
                t += step_size;
                fails = 0;
                dt = std::min(dt * 1.5, opts.dt);        // Ease the step size back up
            } else {
                if (++fails > 2) {
                    // Clamp @f$T>0@f$, @f$Y \in [0,1]@f$
                    sim.resetBadValues(x);
                    fails = 0;
                }
                dt *= 0.5;                               // Shrink and retry,
                                                           // even after a clamp
                if (dt < dt_min) {
                    throw Cantera::CanteraError("solveChannelTransient",
                        "Backward-Euler Newton failed to converge; dt fell below "
                        "{:g} s near t={:g} s. The electron-impact chemistry is too "
                        "stiff for this step -- use integrator='bdf'.", dt_min, t);
                }
            }
            if (++total_substeps > max_substeps) {
                throw Cantera::CanteraError("solveChannelTransient",
                    "Exceeded {} sub-steps; integration is not progressing.",
                    max_substeps);
            }
        }
        record(t);
    }
 
    // Reset every domain's rdt to 0 (Cantera's steady/transient mode flag),
    // matching Sim1D's own convention of leaving the object in steady mode
    // once time-stepping is done -- `sim` is local and destroyed right after,
    // so this has no observable effect here, but keeps the pattern consistent
    // with how Sim1D is used elsewhere (e.g. PlasmaColumnSolver).
    sim.setSteadyMode();
    history.nframes = history.t.size();
    return history;
}
 
} // namespace rizer