Groundwater.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_GROUNDWATER_HPP
#define OPENSWMM_GROUNDWATER_HPP
 
#ifndef OPENSWMM_RESTRICT
#if defined(_MSC_VER)
#define OPENSWMM_RESTRICT __restrict
#else
#define OPENSWMM_RESTRICT __restrict__
#endif
#endif
 
#include <vector>
#include "../math/MathExpr.hpp"
 
namespace openswmm {
 
struct SimulationContext;
 
namespace groundwater {

enum GWVar : int {
    GWV_HGW = 0,   
    GWV_HSW = 1,   
    GWV_HCB = 2,   
    GWV_HGS = 3,   
    GWV_KS = 4,    
    GWV_K = 5,     
    GWV_THETA = 6, 
    GWV_PHI = 7,   
    GWV_FI = 8,    
    GWV_FU = 9,    
    GWV_A = 10,    
    GWV_MAX = 11
};

static const char* GW_VAR_NAMES[] = {"HGW",   "HSW", "HCB", "HGS", "KS", "K",
                                     "THETA", "PHI", "FI",  "FU",  "A"};

static const char* GW_VAR_DESCRIPTIONS[] = {
    "Water table height above aquifer bottom (ft or m)",
    "Surface-water head at the receiving node above aquifer bottom (ft or m)",
    "Channel bottom height above aquifer bottom (Hstar) (ft or m)",
    "Ground surface height above aquifer bottom (ft or m)",
    "Saturated hydraulic conductivity (in/hr or mm/hr)",
    "Unsaturated hydraulic conductivity (in/hr or mm/hr)",
    "Upper-zone moisture content (fraction)",
    "Soil porosity (fraction)",
    "Surface infiltration rate (in/hr or mm/hr)",
    "Upper-zone percolation rate (in/hr or mm/hr)",
    "Subcatchment area (ac or ha)"};

int gwf_validate(const std::string& expr, std::string& msg, int& col);
 
// ============================================================================
// Per-subcatchment GW state (SoA for vectorization)
// ============================================================================
 
struct GWSoA {
    int n_subcatch = 0;
 
    // Aquifer properties (set at init, constant during sim)
    std::vector<double> porosity;
    std::vector<double> field_cap;
    std::vector<double> wilt_point;
    std::vector<double> k_sat;   
    std::vector<double> k_slope; 
    std::vector<double> tension_slope;
    std::vector<double> upper_evap_frac;
    std::vector<int> upper_evap_pat; 
    std::vector<double> lower_evap_depth;
    std::vector<double> lower_loss_coeff; 
    std::vector<double> total_depth;      
    std::vector<double>
        bottom_elev; 
 
    // Lateral flow coefficients
    std::vector<double> a1, b1; 
    std::vector<double> a2, b2; 
    std::vector<double> a3;     
    std::vector<double> h_star; 
 
    // State variables (updated each timestep)
    std::vector<double> theta;       
    std::vector<double> lower_depth; 
 
    // State
    std::vector<double> old_flow; 
 
    // Outputs
    std::vector<double> gw_flow;    
    std::vector<double> upper_evap; 
    std::vector<double> lower_evap; 
    std::vector<double> deep_loss;  
    std::vector<double> max_infil_vol;
 
    // Custom flow expressions (from [GWF] section)
    std::vector<mathexpr::Expression> lateral_expr;
    std::vector<mathexpr::Expression> deep_expr;

    void resize(int n);
};
 
// ============================================================================
// GW solver
// ============================================================================

class GWSolver {
public:
    void init(int n_subcatch);

    void execute(SimulationContext& ctx, double dt, double max_evap, const double* infil_rate,
                 const double* sw_head, const double* frac_perv, const double* perv_evap_rate,
                 int evap_pattern_month = 0);

    GWSoA& state() { return soa_; }
    const GWSoA& state() const { return soa_; }
 
private:
    GWSoA soa_;

    static void batchUpperPerc(const double* OPENSWMM_RESTRICT theta,
                               const double* OPENSWMM_RESTRICT field_cap,
                               const double* OPENSWMM_RESTRICT k_sat,
                               const double* OPENSWMM_RESTRICT k_slope,
                               double* OPENSWMM_RESTRICT perc, int count);

    static void batchGWFlow(const double* OPENSWMM_RESTRICT lower_depth,
                            const double* OPENSWMM_RESTRICT h_star,
                            const double* OPENSWMM_RESTRICT a1, const double* OPENSWMM_RESTRICT b1,
                            const double* OPENSWMM_RESTRICT a2, const double* OPENSWMM_RESTRICT b2,
                            const double* OPENSWMM_RESTRICT a3,
                            const double* OPENSWMM_RESTRICT sw_head,
                            double* OPENSWMM_RESTRICT gw_flow, int count);
};
 
} // namespace groundwater
} // namespace openswmm
 
#endif // OPENSWMM_GROUNDWATER_HPP