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 <vector>
#include <cstddef>
 
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
};
 
// ============================================================================
// 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; 
 
    // 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> 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;    
 
    // 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;      
 
    // 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> 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);
};
 
// ============================================================================
// LID solver
// ============================================================================
 
class LIDSolver {
public:
    void init(SimulationContext& ctx);

    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()); }

    double storedVolume() const;

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

    void execute(SimulationContext& ctx, double dt,
                 double rainfall, double evap_rate);

    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_;
};
 
} // namespace lid
} // namespace openswmm
 
#endif // OPENSWMM_LID_HPP