Snow.hpp

Go to the documentation of this file.

 
#ifndef OPENSWMM_SNOW_HPP
#define OPENSWMM_SNOW_HPP
 
#include <vector>
 
namespace openswmm {
 
struct SimulationContext;
 
namespace snow {
 
// ============================================================================
// Constants
// ============================================================================
 
constexpr int N_SUBAREAS = 3;    
constexpr int SNOW_PLOWABLE = 0;
constexpr int SNOW_IMPERV   = 1;
constexpr int SNOW_PERV     = 2;
 
// ============================================================================
// Per-subcatchment snowpack state (SoA for vectorization)
// ============================================================================
 
struct SnowSoA {
    int n_subcatch = 0;
 
    // Per subcatchment × 3 subareas = flat arrays [subcatch * 3 + subarea]
    std::vector<double> wsnow;    
    std::vector<double> fw;       
    std::vector<double> coldc;    
    std::vector<double> ati;      
    std::vector<double> awe;      
    std::vector<double> imelt;    
 
    // Per subcatchment × 3 subareas: parameters
    std::vector<double> tbase;    
    std::vector<double> dhm;      
    std::vector<double> dhmin;    
    std::vector<double> dhmax;    
    std::vector<double> fwfrac;   
 
    // Per subcatchment × 3 subareas: area fractions
    std::vector<double> fArea;    
 
    // Per subcatchment × 3 subareas: areal depletion state
    std::vector<double> si;       
    std::vector<double> sba;      
    std::vector<double> sbws;     
    std::vector<double> asc;      
 
    // Global ADC curves (shared by all subcatchments, 10 points each)
    double adc_imperv[10] = {1,1,1,1,1,1,1,1,1,1};
    double adc_perv[10]   = {1,1,1,1,1,1,1,1,1,1};
 
    // Per subcatchment: area fractions
    std::vector<double> snn;      
 
    // Per subcatchment: plowing parameters
    std::vector<double> weplow;   
    std::vector<double> sfrac;    
    std::vector<int>    to_subcatch; 
 
    // Global parameters
    double tipm    = 0.5;         
    double rnm     = 0.6;         
    double season  = 0.0;         
    double removed = 0.0;         
 
    void resize(int n);
};
 
// ============================================================================
// Snow solver
// ============================================================================
 
class SnowSolver {
public:
    void init(int n_subcatch);

    void execute(SimulationContext& ctx, double dt,
                 double temp, double wind, const double* rainfall,
                 const double* snowfall,
                 double gamma = 0.0, double ea = 0.0);

    void execute(SimulationContext& ctx, double dt,
                 double temp, double wind, double rainfall,
                 double snowfall = 0.0,
                 double gamma = 0.0, double ea = 0.0);

    void setMeltCoeffs(int day_of_year);

    void plowSnow(SimulationContext& ctx, double dt, const double* snowfall);

    void plowSnow(SimulationContext& ctx, double dt, double snowfall);
 
    SnowSoA& state() { return soa_; }
    const SnowSoA& state() const { return soa_; }
 
private:
    SnowSoA soa_;

    static double rainMeltRate(double temp, double wind, double gamma,
                               double ea, double rainfall);
};
 
} // namespace snow
} // namespace openswmm
 
#endif // OPENSWMM_SNOW_HPP