SubsurfaceData.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_2D_SUBSURFACE_DATA_HPP
#define OPENSWMM_ENGINE_2D_SUBSURFACE_DATA_HPP
 
#include <cstdint>
#include <string>
#include <vector>
#include "SurfaceOwnership.hpp"
 
namespace openswmm::twoD {

enum class SoilChar : int8_t {
    RUSSO = 0,        
    GARDNER = 1,      
    BROOKS_COREY = 2, 
    VAN_GENUCHTEN = 3 
};

enum class GwClosure : int8_t {
    AUTO = -1,       
    CLOSED_FORM = 0, 
    ENSLAVED = 1,    
    SIGMA = 2        
};

struct GwNodeBed {
    int node = -1;     
    int cell = -1;     
    double Kc = 0.0;   
    double dC = 0.0;   
    double area = 0.0; 
    bool locate = false;
    bool exchange = true;
    bool automatic = false;
};

enum class GwLinkMode : int8_t { NONE = 0, AUTO = 1, TWO_WAY = 2, DEFAULT = 3 };

struct GwLinkRow {
    int link = -1;
    double Kc = 0.0; 
    double dC = 0.0; 
    bool exchange = true; 
};

struct GwLinkShare {
    int link = -1;       
    int cell = -1;       
    double weight = 0.0; 
    struct Contact {
        double weight, begin, end;
    }; 
    std::vector<Contact> contacts;
};

struct GwSubcatchShare {
    int subcatch = -1;   
    int cell = -1;       
    double weight = 0.0; 
};

struct GwOptions {
    SoilChar soil_char = SoilChar::RUSSO;
    GwClosure closure = GwClosure::AUTO;
    int m_layers = 8;               
    bool capillary_diff = false;    
    double c_gw = 0.5;              
    double c_col = 0.9;             
    bool force_closed_form = false; 
    bool per_subcatch = false;
    bool dunne = true;
    std::string gw_et = "AUTO";
    double wilting_suction = 150.0; 
    bool wilting_suction_set = false;
    bool node_auto = true;
    GwLinkMode link_seepage = GwLinkMode::DEFAULT;
    bool authored = false;
};

struct GwAquiferRow {
    int scope = 0; 
    std::string tag;
    int cell = -1; 
 
    double Ks = 1.0e-5;    
    double zs = 5.0;       
    double theta_s = 0.45; 
    double theta_r = 0.10; 
    double alpha = 2.0;    
    double psi_b = 0.20;
    double lambda = 0.40;
    double vg_n = 1.6;
    double vg_L = 0.5;
    double c_loss = 0.0;
    double hg0 = -1.0;
 
    SoilChar soil_char = SoilChar::RUSSO;
    GwClosure closure = GwClosure::AUTO;
    int m_layers = -1; 
    bool soil_char_set = false, closure_set = false;
};

struct SubsurfaceState {
    int n_cells = 0;
    int m_layers = 0; 
 
    // ---- resolved parameters (per cell) ---------------------------------
    std::vector<double> Ks, zs, theta_s, theta_r, alpha;
    std::vector<double> psi_b, lambda, vg_n, vg_L, c_loss;
    std::vector<int8_t> soil_char; 
    std::vector<int8_t> closure;   
    std::vector<double> area;      
    std::vector<double> z_bed;     
 
    // ---- state ----------------------------------------------------------
    std::vector<double> hg; 
    std::vector<double> hu; 
    std::vector<double> theta_sigma;
 
    // ---- per-firing diagnostics + transport donors ------------------------
    std::vector<double> q0_last;    
    std::vector<double> qnode_last; 
    std::vector<double> qlat_last;  
    std::vector<double> qdeep_last; 
    std::vector<double> qet_last;   
    std::vector<double> dunne_last; 
    std::vector<double> qplus_last; 
    std::vector<double> infil_capacity, infil_remaining, infil_refresh, infil_interval;
    std::vector<double> wetting_front, reject_last, reject_cumulative;
    std::vector<double> dunne_cumulative;
    // R3: atmospheric demand is integrated at the surface cadence, in m3.
    std::vector<double> et_pending, et_potential_cumulative, et_surface_cumulative;
    std::vector<double> et_soil_cumulative, et_unused_cumulative;
    std::vector<double> et_surface_last, et_potential_last, et_stress, et_refresh;
    std::vector<double> qlink_last; 
 
    // ---- LTS ------------------------------------------------------------
    std::vector<double> dt_cell; 
    std::vector<uint8_t> tier;   
    std::vector<double> eacc_L, eacc_R;
    std::vector<double> xacc_from_surface;
    std::vector<double> xacc_to_surface;
    std::vector<double> nacc;
    std::vector<double> lacc;
 
    // ---- ledger (m³, cumulative) ----------------------------------------
    double led_recharge = 0.0; 
    double led_lateral = 0.0;  
    double led_deep = 0.0;     
    double led_node = 0.0;     
    double led_reject = 0.0;   
    double led_dunne = 0.0;    
    double led_caprise = 0.0;  
    double led_et = 0.0;       
    double led_infil_in = 0.0; 
    double led_link = 0.0;     
    double led_subcatch_in = 0.0;
    double led_init_storage = 0.0;
    double led_final_storage = 0.0;
    double led_source_in = 0.0,
           led_source_out = 0.0; 
 
    double led_boundary_in = 0,
           led_boundary_out = 0; 
    bool active = false;         
    void resize(int n, int m);
    double storage() const noexcept;
    double liveStorage() const noexcept;
    double ledgeredStorage() const noexcept;
    double continuityResidual() const noexcept;
};

struct GwUnitFactors {
    double length = 1.0;  
    double rate = 1.0;    
    double inv_len = 1.0; 
    double area = 1.0;    
};

struct GwBoundaryRow {
    int cell = -1, edge = -1;
    bool head = false;
    double value = 0.0;
    std::string series;
};

struct SubsurfaceConfig {
    GwOptions options;
    std::vector<GwAquiferRow> rows;
    std::vector<GwNodeBed> node_beds;
    std::vector<GwBoundaryRow> boundaries;
    std::vector<GwLinkRow> link_rows;               
    std::vector<SurfaceOwnerRecord> surface_owners; 
    bool empty() const noexcept {
        return !options.authored && rows.empty() && node_beds.empty() && link_rows.empty() &&
               boundaries.empty();
    }

    void clear() { *this = SubsurfaceConfig{}; }
};

const char* soilCharToken(SoilChar s) noexcept;
bool parseSoilChar(const std::string& t, SoilChar& s) noexcept;
const char* gwClosureToken(GwClosure c) noexcept;
bool parseGwClosure(const std::string& t, GwClosure& c) noexcept;
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SUBSURFACE_DATA_HPP