CircuitCoupling.cpp

Go to the documentation of this file.

#include "CircuitCoupling.h"
 
#include "cantera/base/ct_defs.h"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
 
namespace rizer
{
 
    // ------------------------------------------------------------------------------ //
    //   Constructor                                                                  //
    // ------------------------------------------------------------------------------ //
 
    ChannelCircuit::ChannelCircuit(
        const PlasmaChannel1D &channel, double gap,
        const std::function<double(double, double)> &plasma_voltage,
        std::size_t channel_offset, std::size_t n_vars)
        : m_channel(&channel), m_gap(gap),
          m_plasma_voltage(&plasma_voltage),
          m_channel_offset(channel_offset), m_n_vars(n_vars),
          m_n_points(channel.nPoints()),
          m_Rp(std::numeric_limits<double>::infinity())
    {
        if (!active())
        {
            return;
        }
        if (!(m_gap > 0.0))
        {
            throw std::runtime_error(
                "ChannelCircuit: the circuit closes on "
                "R_p = gap / integral(sigma dA) and sets the field to "
                "V_p / gap, so gap must be positive; got " +
                std::to_string(m_gap) + " m.");
        }
        // A single node has no spacing to mirror, so its outer face is a
        // stand-in built from the domain radius rather than a cell edge, and a
        // conductance over it would describe no cross-section in particular.
        if (m_n_points < 2)
        {
            throw std::runtime_error(
                "ChannelCircuit: a circuit closes on integral(sigma dA) over "
                "the cross-section, which needs at least 2 radial nodes to "
                "have an area; got " +
                std::to_string(m_n_points) + ".");
        }
 
        // The same cells the hydro sub-step, the conduction stencil and the
        // Joule source act on, out to the wall face half a spacing beyond the
        // last node (NUMERICS.md section 2). Any other cells would put the
        // circuit's power V^2/R_p and the deposited Joule power out of balance.
        m_cell_area.resize(m_n_points);
        for (std::size_t j = 0; j < m_n_points; j++)
        {
            const double inner = channel.faceRadius(j);
            const double outer = channel.faceRadius(j + 1);
            m_cell_area[j] = Cantera::Pi * (outer * outer - inner * inner);
        }
    }
 
    // ------------------------------------------------------------------------------ //
    //   Public methods                                                               //
    // ------------------------------------------------------------------------------ //
 
    void ChannelCircuit::seedResistance(const std::vector<double> &x)
    {
        if (!active())
        {
            return;
        }
        m_Rp = resistanceFrom(conductance(x));
    }
 
    void ChannelCircuit::closeOver(double t0, double t1,
                                   const std::vector<double> &x)
    {
        if (!active())
        {
            return;
        }
        m_Rp = resistanceFrom(conductance(x));
        m_t0 = t0;
        m_dt = t1 - t0;
        m_Vp0 = (*m_plasma_voltage)(t0, m_Rp);
        m_Vp1 = (*m_plasma_voltage)(t1, m_Rp);
    }
 
    double ChannelCircuit::electricField(double t) const
    {
        // Clamped to the step: CVODE can evaluate past t1 (its internal step
        // overshoots the output time), and an extrapolated V_p there is a
        // voltage the circuit never produced.
        const double fraction =
            (m_dt > 0.0) ? std::min(std::max((t - m_t0) / m_dt, 0.0), 1.0) : 0.0;
        return (m_Vp0 + (m_Vp1 - m_Vp0) * fraction) / m_gap;
    }
 
    CircuitSample ChannelCircuit::sample(const std::vector<double> &x) const
    {
        CircuitSample out;
        out.Rp = m_Rp;
        out.Rp_state = resistanceFrom(conductance(x));
        out.Vp = m_Vp1;
        out.current = (m_Rp > 0.0 && std::isfinite(m_Rp)) ? m_Vp1 / m_Rp : 0.0;
        out.electric_field = m_Vp1 / m_gap;
        return out;
    }
 
    // ------------------------------------------------------------------------------ //
    //   Private methods                                                              //
    // ------------------------------------------------------------------------------ //
 
    double ChannelCircuit::conductance(const std::vector<double> &x) const
    {
        double total = 0.0;
        for (std::size_t j = 0; j < m_n_points; j++)
        {
            total += m_channel->nodeElectricalConductivity(
                         x[m_channel_offset + m_n_vars * j + 0],
                         x[m_channel_offset + m_n_vars * j + 1],
                         &x[m_channel_offset + m_n_vars * j + 2], j) *
                     m_cell_area[j];
        }
        return total;
    }
 
    double ChannelCircuit::resistanceFrom(double conductance_value) const
    {
        return (conductance_value > 0.0)
                   ? m_gap / conductance_value
                   : std::numeric_limits<double>::infinity();
    }
 
} // namespace rizer