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;  
    std::vector<RadiativeConfig> rad_cell; 
    std::vector<SedimentConfig> sed_cell;  
 
    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<unsigned char> mesh_contact_link; 
    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;

    std::vector<double> subcatch_lid_drain_temp_cfs;
 
    static constexpr int kNSubArea = static_cast<int>(HeatSubArea::COUNT_);

    void resize(int n_nodes, int n_links, double initial_temp) {
        // Lazy 1D backend seeding must retain the router's already-bound
        // contact ownership. Model/router initialization clears that mask.
        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);
        subcatch_lid_drain_temp_cfs.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