Reach.h

#ifndef HYDROBRICKS_REACH_H
#define HYDROBRICKS_REACH_H
 
#include "Includes.h"
 
class Subbasin;
struct ChannelRoutingSettings;
 
class Reach {
  public:
    enum class Scheme {
        None,           // instantaneous pass-through
        Lag,            // pure translation by length / celerity
        Muskingum,      // Muskingum method with K = length / celerity and the weighting factor X
        MuskingumCunge  // Muskingum with K and X derived from the channel geometry and the discharge
    };
 
    explicit Reach(Subbasin* subbasin);
 
    virtual ~Reach() = default;
 
    static Scheme SchemeFromString(const string& name);
 
    static string SchemeToString(Scheme scheme);
 
    void Initialize();
 
    void SetChannelRouting(const ChannelRoutingSettings& settings);
 
    double RouteUpstream(double inflowVolume, double timeStepInDays);
 
    double RouteLocal(double localVolume, double timeStepInDays);
 
    [[nodiscard]] bool RoutesLocalRunoff() const {
        return _routeLocal && _scheme != Scheme::None;
    }
 
    void Reset();
 
    void SaveAsInitialState();
 
    [[nodiscard]] Subbasin* GetSubbasin() const {
        return _subbasin;
    }
 
    [[nodiscard]] Scheme GetScheme() const {
        return _scheme;
    }
 
    [[nodiscard]] double GetLength() const {
        return _length;
    }
 
    [[nodiscard]] double GetSlope() const {
        return _slope;
    }
 
    [[nodiscard]] double GetWidth() const;
 
    [[nodiscard]] double GetManning() const;
 
    [[nodiscard]] double GetCelerity() const;
 
    [[nodiscard]] double GetCelerityForDischarge(double dischargeM3s) const;
 
    [[nodiscard]] double GetMuskingumX() const;
 
    [[nodiscard]] double GetTravelTimeInDays() const;
 
    [[nodiscard]] int GetSubreachCount() const {
        return _main.subreaches;
    }
 
    [[nodiscard]] double GetInflow() const {
        return _inflow;
    }
 
    [[nodiscard]] double GetOutflow() const {
        return _outflow;
    }
 
    [[nodiscard]] double GetStorage() const {
        return _main.storage + _local.storage;
    }
 
  protected:
    struct RoutingState {
        // Lag scheme: the delivery schedule (slot j: water due j time steps from now) and its ordinates.
        vecDouble schedule;
        vecDouble ordinates;
        vecDouble initialSchedule;
 
        // Muskingum schemes: the coefficients of the sub-step, the sub-step count and the previous values.
        int substeps = 1;
        double c0 = 1, c1 = 0, c2 = 0;
        bool instantaneous = false;  // the travel time is too short for the step: pass-through
 
        // The sub reaches routed in series (one, unless the Muskingum-Cunge criterion asks for more), each
        // with its own memory of the previous sub-step and its own share of the volume in transit.
        int subreaches = 1;
        vecDouble previousInflow;   // m3 per sub-step, per sub reach
        vecDouble previousOutflow;  // m3 per sub-step, per sub reach
        vecDouble subStorage;       // m3 in transit, per sub reach
        vecDouble initialSubStorage;
        double subreachTimeStep = -1;  // the time step the sub reach count was computed for
 
        double storage = 0;  // m3 in transit, over all the sub reaches
        double initialStorage = 0;
 
        // What the ordinates and coefficients were computed for, to recompute them when it changes.
        double lastTravelTime = -1;
        double lastX = -1;
        double lastTimeStep = -1;
    };
 
    Subbasin* _subbasin;  // non-owning
    Scheme _scheme = Scheme::None;
    bool _routeLocal = false;
    double _length = 0;  // m
    double _slope = 0;   // m/m
 
    // Model-wide parameters (non-owning pointers into the channel routing settings).
    const float* _celerity = nullptr;
    const float* _celerityExponent = nullptr;
    const float* _referenceDischarge = nullptr;
    const float* _x = nullptr;
    const float* _width = nullptr;
    const float* _manning = nullptr;
 
    // Reach-specific values from the subbasin properties (NaN when absent).
    double _celerityOverride = NAN_D;
    double _xOverride = NAN_D;
    double _widthOverride = NAN_D;
    double _manningOverride = NAN_D;
 
    double _inflow = 0;  // m3 per time step
    double _outflow = 0;
 
    RoutingState _main;   // the upstream inflow, over the whole reach
    RoutingState _local;  // the subbasin's own runoff, over half the reach
 
  private:
    double RouteBranch(RoutingState& state, double volume, double length, double timeStepInDays);
 
    double RouteLag(RoutingState& state, double volume, double travelTime, double timeStepInDays);
 
    double RouteMuskingum(RoutingState& state, double volume, double travelTime, double x, double timeStepInDays);
 
    void ComputeCungeParameters(double dischargeM3s, double length, double& travelTime, double& x) const;
 
    [[nodiscard]] int ComputeSubreachCount(double length, double timeStepInDays) const;
 
    static void SetSubreachCount(RoutingState& state, int count);
 
    static void ComputeOrdinates(RoutingState& state, double travelTime, double timeStepInDays);
 
    void ComputeCoefficients(RoutingState& state, double travelTime, double x, double timeStepInDays) const;
 
    static void ResetState(RoutingState& state);
};
 
#endif  // HYDROBRICKS_REACH_H