LID.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_LID_HPP
#define OPENSWMM_LID_HPP
 
#include "Infiltration.hpp"
#include <utility>
#include <vector>
#include <cstddef>
#include <cstdint>
#include <string>
#include <fstream>
#include <map>
#include <memory>
 
namespace openswmm::detail {
class CompletedSourceArchive;
}
 
namespace openswmm {
 
struct SimulationContext;
 
namespace lid {
 
// ============================================================================
// LID type codes
// ============================================================================
 
enum class LIDType : int {
    BIO_CELL = 0,
    RAIN_GARDEN = 1,
    GREEN_ROOF = 2,
    INFIL_TRENCH = 3,
    PERM_PAVEMENT = 4,
    RAIN_BARREL = 5,
    VEG_SWALE = 6,
    ROOF_DISCON = 7
};

inline bool hasNativeBottom(LIDType type, double storage_thickness, double storage_ksat) {
    if (type == LIDType::VEG_SWALE)
        return true;
    if (type != LIDType::BIO_CELL && type != LIDType::RAIN_GARDEN &&
        type != LIDType::INFIL_TRENCH && type != LIDType::PERM_PAVEMENT)
        return false;
    return storage_thickness <= 0 || storage_ksat > 0;
}

inline bool hasNativeBottom(const std::string& type, double thickness, double ksat) {
    if (type == "VS")
        return true;
    return (type == "BC" || type == "RG" || type == "IT" || type == "PP") &&
           (thickness <= 0 || ksat > 0);
}
 
// ============================================================================
// Layer state (per LID unit)
// ============================================================================
 
constexpr int SURF = 0;
constexpr int SOIL = 1;
constexpr int STOR = 2;
constexpr int PAVE = 3;
constexpr int N_LAYERS = 4;
 
// ============================================================================
// LID unit parameters (SoA for batch processing)
// ============================================================================
 
struct LIDGroupSoA {
    LIDType type = LIDType::BIO_CELL;
    int count = 0;
 
    // Unit-conversion factors for the model's flow units (legacy UCF()).
    // Layer parameters are converted to internal ft / ft-per-sec at init();
    // these two factors also drive the underdrain head/rate conversion in
    // getDrainRate(). Defaults are the US-customary values (issue #102).
    double ucf_raindepth = 12.0;   
    double ucf_rainfall = 43200.0; 
 
    std::vector<int> subcatch_idx;   
    std::vector<int> control_idx;    
    std::vector<double> area;        
    std::vector<double> from_imperv; 
    std::vector<double> from_perv;   
    std::vector<int> to_perv;        
    std::vector<int> drain_node;     
    std::vector<int> drain_subcatch; 
    std::vector<double> inflow;      
    std::vector<double>
        evap_rate_unit; 
    std::vector<double>
        subcatch_rain; 
    std::vector<std::string>
        rpt_path; 
    // Per-step kernel rates staged for the LID report file (ft/s). Legacy
    // lidproc_saveResults reads these straight from lidproc.c's static flux
    // variables; the kernel here is per-call, so runUnitLegacy stages them.
    std::vector<double> rpt_surf_infil; 
    std::vector<double> rpt_pave_perc;  
    std::vector<double> rpt_soil_perc;  
 
    // Surface layer
    std::vector<double> surf_store; 
    std::vector<double> surf_rough; 
    std::vector<double> surf_slope; 
 
    // Soil layer
    std::vector<double> soil_thick;   
    std::vector<double> soil_poros;   
    std::vector<double> soil_fc;      
    std::vector<double> soil_wp;      
    std::vector<double> soil_ksat;    
    std::vector<double> soil_kslope;  
    std::vector<double> soil_suction; 
 
    // Storage layer
    std::vector<double> stor_thick; 
    std::vector<double> stor_void;  
    std::vector<double> stor_ksat;  
    std::vector<double> stor_clog;  
    std::vector<int> stor_covered;  
 
    // Drain
    std::vector<double> drain_coeff;  
    std::vector<double> drain_expon;  
    std::vector<double> drain_offset; 
    std::vector<double> drain_delay;  
    std::vector<double> drain_hopen;  
    std::vector<double> drain_hclose; 
    std::vector<int> drain_open;      
 
    // Pavement layer (PERM_PAVEMENT)
    std::vector<double> pave_thick;       
    std::vector<double> pave_void;        
    std::vector<double> pave_imperv_frac; 
    std::vector<double> pave_ksat;        
    std::vector<double> pave_clog_factor; 
    std::vector<double> pave_regen_days;  
    std::vector<double> pave_regen_deg;   
    std::vector<double> next_regen_day;   
 
    // Drainage mat layer (GREEN_ROOF)
    std::vector<double> drainmat_thick; 
    std::vector<double> drainmat_void;  
    std::vector<double> drainmat_rough; 
 
    // Surface geometry (for Manning's outflow)
    std::vector<double> surf_void_frac;  
    std::vector<double> surf_alpha;      
    std::vector<double> surf_side_slope; 
    std::vector<double> full_width;      
    std::vector<double> unit_area;       
    std::vector<double> unit_width;      
    std::vector<double> dry_time;        
 
    // State variables (updated each step)
    std::vector<double> surf_depth; 
    std::vector<double> soil_moist; 
    std::vector<double> stor_depth; 
    std::vector<double> pave_depth; 
 
    // Legacy lidproc.c per-unit state the kernel needs beyond the layer
    // depths (lid.c TLidUnit): the unit's own Green-Ampt infiltration
    // state for the soil layer (soilInfil), the drain flow the previous
    // step delivered (oldDrainFlow, the underdrain's hOpen/hClose
    // hysteresis memory), and the process's canOverflow flag
    // (validateLidProc: false for a roof disconnection and for a
    // bio-cell / rain garden / trench / pavement / green roof whose surface
    // has a Manning alpha).
    std::vector<GreenAmptState> soil_infil;
    std::vector<double> old_drain_flow;
    std::vector<uint8_t>
        is_wet; 
    std::vector<uint8_t> can_overflow;
    std::vector<double>
        drainmat_alpha; 
 
    // Outputs (per unit)
    std::vector<double> surface_runoff;
    std::vector<double> drain_flow;
    std::vector<double> evap_loss;  
    std::vector<double> infil_loss; 
 
    // Pollutant drain removal fractions: drain_rmvl[unit * n_pollutants + pollutant]
    std::vector<double> drain_rmvl; 
    int n_pollutants = 0;           
 
    // ---- A4: per-layer INFLOW rates (ft/sec of water per unit area) ----
    // Each layer receives from exactly one place: the layer above it in the
    // present stack, or externally for the topmost present layer. Every
    // batch*Flux routine already computes these as locals (soil_infil,
    // soil_perc, pavePerc, storageInflow, ...); they are published here
    // because a complete-mix age needs the INFLOW and nothing else — the
    // outflow leaves at the layer's own age.
    //
    // These are NOT the `f_old_*` below. Those are the Modified Puls
    // time-weighting term: `f_old_surf` is the NET dx/dt of the surface
    // layer, and the other three are never written at all. A net rate of
    // change is exactly the quantity that made phase A3 report elapsed time
    // instead of age.
    std::vector<double> in_surf; 
    std::vector<double> in_pave; 
    std::vector<double> in_soil; 
    std::vector<double> in_stor; 
    std::vector<double> evap_surf, evap_pave, evap_soil, evap_stor;
 
    // Previous flux rates (for Modified Puls time weighting)
    std::vector<double> f_old_surf; 
    std::vector<double> f_old_soil; 
    std::vector<double> f_old_stor; 
    std::vector<double> f_old_pave; 
 
    // Water balance tracking (cumulative per unit)
    std::vector<double> wb_inflow; 
    std::vector<double> wb_evap;   
    std::vector<double> spatial_infil_loss, wb_spatial_infil;
    std::vector<double> wb_infil;      
    std::vector<double> wb_surf_flow;  
    std::vector<double> wb_drain_flow; 
    std::vector<double> wb_init_vol;   
    std::vector<double> wb_final_vol;  
    std::vector<double> vol_treated;   

    void resize(int n);
};

double completedUnitStorage(const LIDGroupSoA& g, int unit);
 
// ============================================================================
// LID solver
// ============================================================================
 
class LIDSolver {
public:
    void init(SimulationContext& ctx);
 
    struct CompletedUnitInput {
        int type = -1, unit = -1;
        double inflow = 0.0, rain = 0.0, pet = 0.0;
        double native_infil = 0.0, max_native_infil = 1.0e10, infil_factor = 1.0;
        // Negative fraction retains ordinary completed-kernel arithmetic.
        double spatial_fraction = -1.0, spatial_max_infil = 1.0e10;
    };

    LIDSolver waterTrial() const;
    void executeCompleted(double dt, double start, double recovery_factor,
                          const std::vector<CompletedUnitInput>& inputs);

    LIDGroupSoA& group(int type_index) { return groups_[static_cast<size_t>(type_index)]; }
    const LIDGroupSoA& group(int type_index) const {
        return groups_[static_cast<size_t>(type_index)];
    }

    int numGroups() const { return static_cast<int>(groups_.size()); }

    const std::vector<std::pair<int, int>>& usageOrder() const { return usage_order_; }

    double storedVolume(const std::vector<int>* source_owners = nullptr) const;
    void checkpoint(detail::CompletedSourceArchive& a, bool model,
                    bool lazy_buffers = false); // internal accepted-boundary IO
    void validateCheckpoint(const LIDSolver& model, std::size_t runoff_count) const;

    double totalStoredVolume() const;
    double totalInitVolume() const;
    double totalInfilVolume() const;
    double totalEvapVolume() const;

    void execute(SimulationContext& ctx, double dt, double rainfall, double evap_rate,
                 const std::vector<int>* source_owners = nullptr);

    void setRunoffTime(double old_runoff_sec) { old_runoff_sec_ = old_runoff_sec; }

    void setNativeInfil(std::vector<double> native, std::vector<double> max_native) {
        native_infil_ = std::move(native);
        max_native_infil_ = std::move(max_native);
    }

    void setInfilFactors(std::vector<double> infil_factor, double recovery_factor) {
        infil_factor_ = std::move(infil_factor);
        recovery_factor_ = recovery_factor;
    }

    bool anyWet(const std::vector<int>* source_owners = nullptr) const;

    static void batchBioCellFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                 double dt);

    static void batchBarrelFlux(LIDGroupSoA& g, double rainfall, double dt);

    static void batchSwaleFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate, double dt);

    static void batchGreenRoofFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                   double dt);

    static void batchInfilTrenchFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                     double dt);

    static void batchPavementFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                  double dt);

    static void batchRoofDisconFlux(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                    double dt);

    static void batchSwaleModPuls(LIDGroupSoA& g, double rainfall, const double* evap_rate,
                                  double dt);
 
private:
    std::vector<LIDGroupSoA> groups_;
    std::vector<std::pair<int, int>> usage_order_; 
    double old_runoff_sec_ = 0.0;

    struct LidRptFile {
        std::ofstream f;
        int was_dry = 1;
        std::string held;
        bool failed = false;
    };
    std::map<std::pair<int, int>, std::unique_ptr<LidRptFile>> rpt_files_;
    void saveLidReport(SimulationContext& ctx, int gi, LIDGroupSoA& g, std::size_t u, double dt);
    double recovery_factor_ = 1.0;
    std::vector<double> native_infil_;
    std::vector<double> max_native_infil_;
    std::vector<double> infil_factor_;
};
 
} // namespace lid
} // namespace openswmm
 
#endif // OPENSWMM_LID_HPP