PropertyTable.h

Go to the documentation of this file.

// PropertyTable: an immutable 1D piecewise-linear lookup with endpoint
// clamping.
//
// The abscissa is whatever the caller tabulated against, not necessarily a
// temperature. The same container holds, in this extension:
//   - sigma(T), kappa(T), P_rad(T) and the optional seed profile T_seed(r)
//     of ThermalPlasmaColumn1D (steady LTE column);
//   - the sigma(Te) and nu_m(Te) fallbacks of PlasmaCell, and CollisionModel's
//     Maxwellian kernel Qbar(Te) together with the momentum-transfer cross
//     section sigma(energy) it is built from;
//   - kappa(T), kappa_e(T) and the initial radial profiles Tg(r), Te(r),
//     Y_s(r), rho(r), u_r(r) of PlasmaChannel1D;
//   - the applied field E(t) of the transient solver.
// eval()'s argument and return carry the units of whatever was tabulated;
// the class itself is unit-agnostic.
//
// Three shapes, decided in one place. fromTableOrScalar() maps a caller's
// (grid, values) pair onto exactly one of
//   - a real table, from >= 2 matching points;
//   - a constant, from one value against an empty or one-point grid;
//   - an empty table, when neither is, whose eval() returns 0.0 -- which is
//     how "this property is not configured" is spelled.
// Any other pair throws, naming the pair. It has to: an empty table is
// indistinguishable from a property that is genuinely zero everywhere, so
// returning one for an off-by-one in the caller's arrays would turn a setup
// error into a silently different physical problem rather than a failure.
//
// Ordering is checked at construction, and is non-decreasing rather than
// strictly increasing. A repeated knot is legal and is how a step profile is
// written (two values at one abscissa), which NUMERICS.md section "Validation ladder"'s P3
// initial condition needs. A decreasing knot throws, because eval() would not
// fail on one: std::upper_bound returns an arbitrary bracket there, and the
// interpolation weight leaves [0, 1] for the rest of the run.
//
// Outside [front, back] eval() clamps to the endpoint value: no
// extrapolation, no NaN, and no error either. The transient closures evaluate
// their tables through numeric/TableLookup.h, which refuses a state outside
// them, and PlasmaChannel1D checks its initial profiles at construction. Three
// lookups still clamp: ThermalPlasmaColumn1D's Newton residual (its converged
// column is checked after the solve instead), CollisionModel::meanCrossSection's
// sigma(energy), whose integral runs past any finite table by design, and the
// LTE column's seed profile.
//
// Nothing here reads a file. Every table is built by the caller, either from
// numpy arrays through bindings.cpp's makeTable() or from a Config pair
// through fromTableOrScalar().
 
#ifndef RIZER_PROPERTY_TABLE_H
#define RIZER_PROPERTY_TABLE_H
 
#include <vector>
#include <stdexcept>
#include <algorithm>
#include <string>
 
namespace rizer
{
 
    class PropertyTable
    {
    public:
        PropertyTable() = default;

        PropertyTable(std::vector<double> T, std::vector<double> values)
            : m_T(std::move(T)), m_v(std::move(values))
        {
            if (m_T.size() != m_v.size() || m_T.size() < 2)
            {
                throw std::invalid_argument(
                    "PropertyTable: T and values must have equal length >= 2");
            }
            // eval() brackets with std::upper_bound and clamps on front()/back(),
            // both of which need the knots ordered. An unsorted grid does not
            // fail there -- it returns a value off an arbitrary bracket, with the
            // interpolation weight outside [0, 1], for the rest of the run.
            //
            // Non-decreasing, not strictly increasing: a REPEATED knot is how a
            // step profile is written (two values at the same radius), which the
            // step-function initial condition of NUMERICS.md section "Validation ladder" needs.
            // Repeats cannot divide by zero below, because upper_bound returns the
            // first knot strictly greater than T, so m_T[hi] > T >= m_T[hi-1]
            // forces m_T[hi] != m_T[hi-1].
            for (std::size_t i = 1; i < m_T.size(); i++)
            {
                if (m_T[i] < m_T[i - 1])
                {
                    throw std::invalid_argument(
                        "PropertyTable: T must be non-decreasing, but T[" +
                        std::to_string(i) + "] = " + std::to_string(m_T[i]) +
                        " is below T[" + std::to_string(i - 1) +
                        "] = " + std::to_string(m_T[i - 1]) + ".");
                }
            }
        }

        static PropertyTable constant(double value)
        {
            return PropertyTable({0.0, 1.0e9}, {value, value});
        }

        static PropertyTable fromTableOrScalar(const std::vector<double> &x,
                                               const std::vector<double> &v,
                                               const char *name = "property table")
        {
            if (x.size() >= 2 && x.size() == v.size())
                return PropertyTable(x, v);
            // A grid of two or more points against one value is the mismatched
            // pair below, not a constant: dropping the grid would hand back a
            // different property than the one tabulated.
            if (v.size() == 1 && x.size() <= 1)
                return PropertyTable::constant(v[0]);
            if (x.empty() && v.empty())
                return PropertyTable();
            throw std::invalid_argument(
                std::string("PropertyTable: ") +
                name +
                " has a grid of " +
                std::to_string(x.size()) +
                " point(s) against " +
                std::to_string(v.size()) +
                " value(s). They must be the same length and at least 2 points, or a "
                "single value for a constant, or both empty to leave the property unset.");
        }

        bool empty() const { return m_T.empty(); }

        bool covers(double x) const
        {
            return !m_T.empty() && x >= m_T.front() && x <= m_T.back();
        }

        double first() const { return m_T.front(); }

        double last() const { return m_T.back(); }

        double eval(double T) const
        {
            if (m_T.empty())
            {
                return 0.0;
            }
            if (T <= m_T.front())
            {
                return m_v.front();
            }
            if (T >= m_T.back())
            {
                return m_v.back();
            }
            // upper_bound gives the first node strictly greater than T; subtract 1
            // to get the left bracket of the containing interval.
            auto it = std::upper_bound(m_T.begin(), m_T.end(), T);
            size_t hi = static_cast<size_t>(it - m_T.begin());
            size_t lo = hi - 1;
            double w = (T - m_T[lo]) / (m_T[hi] - m_T[lo]);
            return m_v[lo] + w * (m_v[hi] - m_v[lo]);
        }
 
    private:
        std::vector<double> m_T; 
        std::vector<double> m_v; 
    };
 
} // namespace rizer
 
#endif