HeatData.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_ENGINE_DATA_HEAT_DATA_HPP
#define OPENSWMM_ENGINE_DATA_HEAT_DATA_HPP
 
#include <limits>
#include <vector>
 
#include "BedZoneData.hpp"      // SedimentConfig — H6b's bed zone configuration
#include "HeatOverrideData.hpp" // HeatElement / HeatAttr / HeatScope — PE
 
namespace openswmm {

enum class HeatSource : int {
    RAINFALL        = 0,  
    DWF             = 1,
    GW              = 2,
    RDII            = 3,
    EXTERNAL_INFLOW = 4,
    IFACE           = 5,
    INITIAL_STATE   = 6,  
    COUNT_          = 7
};

enum class HeatSubArea : int {
    IMPERV0 = 0,
    IMPERV1 = 1,
    PERV    = 2,
    COUNT_  = 3
};

enum class DryTempPolicy : int {
    HOLD    = 0,
    AIR     = 1,
    DEFAULT = 2
};
 

struct ConductionConfig {
    double water_conductivity = 0.606;   
    double sed_conductivity   = 2.6;     
    double sed_density        = 1970.0;  
    double sed_specific_heat  = 2758.0;  
};

enum class ShortwaveMode : int {
    CONSTANT   = 0,  
    TIMESERIES = 1,  
    COMPUTED   = 2   
};

struct SolarConfig {
    double latitude_deg   = 0.0;    
    double longitude_deg  = 0.0;    
    double timezone_hours = 0.0;    
    bool   has_latitude   = false;  
    bool   has_longitude  = false;
    bool   has_timezone   = false;

    double elevation_m    = 0.0;
    bool   has_elevation  = false;
 
    // ---- Bird & Hulstrom (1981) atmosphere. Defaults are the paper's
    //      standard atmosphere.
    double aod380        = 0.30;   
    double aod500        = 0.20;   
    double precip_water_cm = 1.42; 
    double ozone_cm        = 0.34; 
    double ground_albedo   = 0.20; 
};

struct CloudConfig {
    bool   configured = false;  
    bool   use_timeseries = false;
    int    ts_index   = -1;     
    double fraction   = 0.0;    
 
    double sw_atten_k = 0.75;   
    double sw_atten_n = 3.4;    
    double lw_cloud_k = 0.17;   
};

struct RadiativeConfig {
    double shortwave_wm2   = 0.0;   
    double albedo          = 0.0;   
    double shade_factor    = 0.0;   
    double sky_view        = 1.0;   
    double emiss_water     = 0.97;  
    double emiss_landcover = 0.97;  
    double atm_emiss_coeff = 0.5;   
    double lw_reflection   = 0.03;  

    double landcover_temp = std::numeric_limits<double>::quiet_NaN();

    ShortwaveMode sw_mode = ShortwaveMode::CONSTANT;
    int sw_ts_index = -1;
};

struct HeatOverrideData {
    std::vector<HeatOverrideRow> rows;   
 
    std::vector<RadiativeConfig> rad_link;   
    std::vector<RadiativeConfig> rad_node;   
    std::vector<SedimentConfig>  sed_link;   
 
    bool resolved = false;
 
    bool empty() const noexcept { return rows.empty(); }
    void clear() { *this = HeatOverrideData{}; }
};

struct HeatConfigData {
    bool configured = false;

    bool radiative_exchange = false;

    RadiativeConfig radiative;

    bool surface_exchange = false;

    DryTempPolicy dry_temp_policy = DryTempPolicy::HOLD;

    bool layer_conduction = false;

    ConductionConfig conduction;

    bool sediment_exchange = false;

    SedimentConfig sediment;

    HeatOverrideData overrides;

    SolarConfig solar;

    CloudConfig cloud;

    static constexpr double kDefaultTemp = 20.0;

    double global_temp[static_cast<int>(HeatSource::COUNT_)] = {
        kDefaultTemp, kDefaultTemp, kDefaultTemp, kDefaultTemp,
        kDefaultTemp, kDefaultTemp, kDefaultTemp};

    bool configured_source[static_cast<int>(HeatSource::COUNT_)] = {};
 
    // NODE-scope overrides (DWF / EXTERNAL_INFLOW), parallel arrays.
    std::vector<int>    node_over_source;  
    std::vector<int>    node_over_node;    
    std::vector<double> node_over_temp;    

    double source_temp(HeatSource s, int node) const noexcept {
        for (std::size_t i = 0; i < node_over_source.size(); ++i)
            if (node_over_source[i] == static_cast<int>(s) &&
                node_over_node[i] == node)
                return node_over_temp[i];
        return global_temp[static_cast<int>(s)];
    }
};

struct HeatState {
    std::vector<double> node_temp_vol_in;  
    std::vector<double> node_lid_drain_temp_vol_in;
    std::vector<double> node_temp;         
    std::vector<double> link_temp;         

    bool legacy_seeded = false;
 
    // ---- H6a per-step solar forcing (plan §2.5). RESOLVED ONCE PER STEP by
    //      `updateSolarForcing`, then read const by every flux call.
    //
    //      This is state, not config, which is why it lives here: `Jin`
    //      under TIMESERIES or COMPUTED changes every step, while
    //      `RadiativeConfig::shortwave_wm2` is what the deck wrote and must
    //      not be overwritten (a hot-started or re-opened model would
    //      otherwise resume from a stale interpolation rather than from its
    //      own configuration).
    //
    //      It is also what keeps the SPA cost argument honest: shortwave is
    //      GLOBAL scope, so the position is computed once per step and every
    //      element reads the same cached number. A per-element solve is
    //      plan §7 work and would not use this field.

    double shortwave_now = -1.0;
    double cloud_now = 0.0;
 
    // ---- H5a watershed rows. Sized by `resizeWatershed`, NOT by `resize`.
    //      Kept a separate call deliberately: A3 widened `WaterAgeState::
    //      resize` with a defaulted third parameter, four existing two-
    //      argument call sites then silently emptied the new arrays at
    //      runtime, and the validator had to remove the default to move the
    //      failure back to the compiler. A second function cannot be called
    //      short, so the hazard is unrepresentable rather than observed.

    std::vector<double> subarea_temp;
    std::vector<double> subarea_vol_prev;
    std::vector<double> subcatch_runoff_temp;
    std::vector<double> subcatch_runon_temp_vol_in;

    std::vector<double> subcatch_runon_temp_rate;

    std::vector<double> subcatch_outfall_temp_vol;
 
    static constexpr int kNSubArea = static_cast<int>(HeatSubArea::COUNT_);
 
    void resize(int n_nodes, int n_links, double initial_temp) {
        node_temp_vol_in.assign(static_cast<std::size_t>(n_nodes), 0.0);
        node_lid_drain_temp_vol_in.assign(static_cast<std::size_t>(n_nodes), 0.0);
        node_temp.assign(static_cast<std::size_t>(n_nodes), initial_temp);
        link_temp.assign(static_cast<std::size_t>(n_links), initial_temp);
        legacy_seeded = false;
        // H6a. `clear()` resets these via whole-struct assignment, but
        // `resize()` is what runs at INITIALIZE — so a re-initialize on an
        // already-run context would otherwise leave the previous run's
        // forcing readable, and `swmm_heat_get_current_shortwave` documents
        // itself as unresolved before the first step.
        shortwave_now = -1.0;
        cloud_now     = 0.0;
    }
 
    void resizeWatershed(int n_subcatch, double initial_temp) {
        const auto n = static_cast<std::size_t>(n_subcatch);
        subarea_temp.assign(n * static_cast<std::size_t>(kNSubArea),
                            initial_temp);
        subarea_vol_prev.assign(n * static_cast<std::size_t>(kNSubArea), 0.0);
        subcatch_runoff_temp.assign(n, initial_temp);
        subcatch_runon_temp_vol_in.assign(n, 0.0);
        subcatch_runon_temp_rate.assign(n, 0.0);
        subcatch_outfall_temp_vol.assign(n, 0.0);
    }
 
    bool watershedSized(int n_subcatch) const noexcept {
        return subarea_temp.size() == static_cast<std::size_t>(n_subcatch) *
                                          static_cast<std::size_t>(kNSubArea);
    }
 
    void clear() { *this = HeatState{}; }
};
 
}  // namespace openswmm
 
#endif  // OPENSWMM_ENGINE_DATA_HEAT_DATA_HPP