Groundwater.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_GROUNDWATER_HPP
#define OPENSWMM_GROUNDWATER_HPP
 
#ifndef OPENSWMM_RESTRICT
#  if defined(_MSC_VER)
#    define OPENSWMM_RESTRICT __restrict
#  else
#    define OPENSWMM_RESTRICT __restrict__
#  endif
#endif
 
#include <vector>
#include "../math/MathExpr.hpp"
 
namespace openswmm {
 
struct SimulationContext;
 
namespace groundwater {
 
enum GWVar : int {
    GWV_HGW   = 0,  
    GWV_HSW   = 1,  
    GWV_HCB   = 2,  
    GWV_HGS   = 3,  
    GWV_KS    = 4,  
    GWV_K     = 5,  
    GWV_THETA = 6,  
    GWV_PHI   = 7,  
    GWV_FI    = 8,  
    GWV_FU    = 9,  
    GWV_A     = 10, 
    GWV_MAX   = 11
};
 
static const char* GW_VAR_NAMES[] = {
    "HGW", "HSW", "HCB", "HGS", "KS", "K",
    "THETA", "PHI", "FI", "FU", "A"
};
 
// ============================================================================
// Per-subcatchment GW state (SoA for vectorization)
// ============================================================================
 
struct GWSoA {
    int n_subcatch = 0;
 
    // Aquifer properties (set at init, constant during sim)
    std::vector<double> porosity;
    std::vector<double> field_cap;
    std::vector<double> wilt_point;
    std::vector<double> k_sat;        
    std::vector<double> k_slope;      
    std::vector<double> tension_slope;
    std::vector<double> upper_evap_frac;
    std::vector<int>    upper_evap_pat;   
    std::vector<double> lower_evap_depth;
    std::vector<double> lower_loss_coeff;   
    std::vector<double> total_depth;        
 
    // Lateral flow coefficients
    std::vector<double> a1, b1;       
    std::vector<double> a2, b2;       
    std::vector<double> a3;           
    std::vector<double> h_star;       
 
    // State variables (updated each timestep)
    std::vector<double> theta;        
    std::vector<double> lower_depth;  
 
    // State
    std::vector<double> old_flow;     
 
    // Outputs
    std::vector<double> gw_flow;      
    std::vector<double> upper_evap;   
    std::vector<double> lower_evap;   
    std::vector<double> deep_loss;    
    std::vector<double> max_infil_vol;
 
    // Custom flow expressions (from [GWF] section)
    std::vector<mathexpr::Expression> lateral_expr;
    std::vector<mathexpr::Expression> deep_expr;
 
    void resize(int n);
};
 
// ============================================================================
// GW solver
// ============================================================================
 
class GWSolver {
public:
    void init(int n_subcatch);
 
    void execute(SimulationContext& ctx, double dt, double max_evap,
                 const double* infil_rate, const double* sw_head,
                 const double* frac_perv, const double* perv_evap_rate);
 
    GWSoA&       state()       { return soa_; }
    const GWSoA& state() const { return soa_; }
 
private:
    GWSoA soa_;
 
    static void batchUpperPerc(
        const double* OPENSWMM_RESTRICT theta,
        const double* OPENSWMM_RESTRICT field_cap,
        const double* OPENSWMM_RESTRICT k_sat,
        const double* OPENSWMM_RESTRICT k_slope,
        double*       OPENSWMM_RESTRICT perc,
        int count
    );
 
    static void batchGWFlow(
        const double* OPENSWMM_RESTRICT lower_depth,
        const double* OPENSWMM_RESTRICT h_star,
        const double* OPENSWMM_RESTRICT a1, const double* OPENSWMM_RESTRICT b1,
        const double* OPENSWMM_RESTRICT a2, const double* OPENSWMM_RESTRICT b2,
        const double* OPENSWMM_RESTRICT a3,
        const double* OPENSWMM_RESTRICT sw_head,
        double*       OPENSWMM_RESTRICT gw_flow,
        int count
    );
};
 
} // namespace groundwater
} // namespace openswmm
 
#endif // OPENSWMM_GROUNDWATER_HPP