DegreeDaySnow.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_SNOW_HPP
#define OPENSWMM_SNOW_HPP
 
#include <cstddef>
#include <vector>
#include "SnowModel.hpp"
 
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;
    // Native UEB columns publish into the reporting arrays below, but their
    // water must never also be advanced by the degree-day solver.
    std::vector<bool> external;
    std::vector<double> external_net; 
 
    // 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; 

    std::vector<double> age;
    std::vector<double>
        ueb_transfer_energy;            // pending DD-to-UEB plow heat, kJ/m2 of receiving subarea
    std::vector<unsigned char> enabled; // driver ownership mask, true by default
    bool conservative_plowing =
        false; // completed source footprints use actual donor/receiver areas
    std::vector<double> out_age;
 
    // 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; 

    double precip_age = 0.0;
    bool track_age = false;

    void resize(int n);
};
 
// ============================================================================
// Snow solver
// ============================================================================
 
class DegreeDaySnow : public surface::SnowModel {
public:
    void init(int n_subcatch);
    void initializeElements(int n) override { init((n + N_SUBAREAS - 1) / N_SUBAREAS); }
    void step(int e, const surface::SnowInputs& input, surface::SnowOutputs& output) override;
    int stateSize() const override { return 12; }
    void pack(int e, double* output) const override;
    void unpack(int e, const double* input) override;

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

    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) override;

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

    SnowSoA& state() { return soa_; }
    const SnowSoA& state() const { return soa_; }
 
private:
    SnowSoA soa_;
    void meltElement(std::size_t ui, int subarea, double dt, double temp, double rain, double snow,
                     double rmelt);

    static double rainMeltRate(double temp, double wind, double gamma, double ea, double rainfall);
};
 
using SnowSolver = DegreeDaySnow;
} // namespace snow
} // namespace openswmm
 
namespace openswmm::surface {
using DegreeDaySnow = snow::DegreeDaySnow;
}
 
#endif // OPENSWMM_SNOW_HPP