From 3547ed1a67a6b9610f76c9f493a962209f9ad8b3 Mon Sep 17 00:00:00 2001 From: Corinne Wiesner-Friedman Date: Mon, 6 Apr 2026 12:41:51 -0400 Subject: [PATCH] Add AMM RDII integration --- src/engine/core/SWMMEngine.cpp | 14 + src/engine/core/SWMMEngine.hpp | 2 + src/engine/core/SimulationContext.hpp | 1 + src/engine/core/SimulationOptions.hpp | 14 + src/engine/data/InflowData.hpp | 40 +++ src/engine/hydrology/AMM.cpp | 285 +++++++++++++++++++ src/engine/hydrology/AMM.hpp | 189 ++++++++++++ src/engine/input/handlers/InflowsHandler.cpp | 90 +++++- src/engine/input/handlers/InflowsHandler.hpp | 3 + src/engine/input/handlers/OptionsHandler.cpp | 7 + src/engine/plugins/DefaultInputPlugin.cpp | 1 + 11 files changed, 645 insertions(+), 1 deletion(-) create mode 100644 src/engine/hydrology/AMM.cpp create mode 100644 src/engine/hydrology/AMM.hpp diff --git a/src/engine/core/SWMMEngine.cpp b/src/engine/core/SWMMEngine.cpp index 545906de8..b5a307843 100644 --- a/src/engine/core/SWMMEngine.cpp +++ b/src/engine/core/SWMMEngine.cpp @@ -1220,6 +1220,17 @@ void SWMMEngine::stepRouting(double dt_routing) noexcept { rdii_.applyRdiiInflows(ctx_); } + // B2a2. RDII inflows (antecedent moisture model) — AMM method + if (ctx_.options.rdii_method == RdiiMethod::AMM && amm_.componentCount() > 0) { + double avg_rainfall = 0.0; + for (int g = 0; g < ctx_.n_gages(); ++g) { + avg_rainfall += ctx_.gages.rainfall[static_cast(g)]; + } + if (ctx_.n_gages() > 0) avg_rainfall /= ctx_.n_gages(); + double air_temp = climate_.temperature; // current air temperature (deg F) + amm_.computeAll(ctx_, avg_rainfall, air_temp, dt_routing); + } + // B2b. Interface file inflows (from upstream model coupling) iface_.readInflows(ctx_, ctx_.current_date); @@ -2446,6 +2457,9 @@ void SWMMEngine::initHydraulics() noexcept { // 10a. RDII solver: initialize unit hydrograph groups rdii_.init(ctx_); + // 10a2. AMM solver: initialize antecedent moisture model groups + amm_.init(ctx_); + // 10b. Non-conduit hydraulic structures (pumps, orifices, weirs, outlets) hydstruct_.init(ctx_); diff --git a/src/engine/core/SWMMEngine.hpp b/src/engine/core/SWMMEngine.hpp index f2a60692d..14a373360 100644 --- a/src/engine/core/SWMMEngine.hpp +++ b/src/engine/core/SWMMEngine.hpp @@ -49,6 +49,7 @@ #include "../hydrology/LID.hpp" #include "../hydrology/Inflow.hpp" #include "../hydrology/RDII.hpp" +#include "../hydrology/AMM.hpp" #include "../quality/QualityRouting.hpp" #include "../quality/Landuse.hpp" #include "../controls/Controls.hpp" @@ -236,6 +237,7 @@ class SWMMEngine { controls::ControlEngine controls_; ///< Control rule evaluation inflow::InflowSolver inflow_; ///< External + DWF inflows rdii::RDIISolver rdii_; ///< RDII (unit hydrograph convolution) + amm::AMMSolver amm_; ///< AMM (antecedent moisture model) exfil::ExfilSolver exfil_; ///< Storage node exfiltration inlet::InletSolver inlet_; ///< Street inlet capture std::vector culvert_links_;///< Pre-built culvert link indices (avoid per-timestep alloc) diff --git a/src/engine/core/SimulationContext.hpp b/src/engine/core/SimulationContext.hpp index 73516707c..9c16f8234 100644 --- a/src/engine/core/SimulationContext.hpp +++ b/src/engine/core/SimulationContext.hpp @@ -313,6 +313,7 @@ struct SimulationContext { DwfData dwf_inflows; RDIIAssignData rdii_assigns; UnitHydData unit_hyds; ///< Parsed [HYDROGRAPHS] data + AMMAssignData amm_assigns; PatternData patterns; // ========================================================================= diff --git a/src/engine/core/SimulationOptions.hpp b/src/engine/core/SimulationOptions.hpp index c78e65f11..bc367437a 100644 --- a/src/engine/core/SimulationOptions.hpp +++ b/src/engine/core/SimulationOptions.hpp @@ -35,6 +35,17 @@ namespace openswmm { // Enumerators matching legacy SWMM enums.h // ============================================================================ +/** + * @brief RDII computation method. + * + * RTK uses the traditional 3-triangle unit hydrograph convolution. + * AMM uses the Antecedent Moisture Model (Edgren et al., JWMM C525). + */ +enum class RdiiMethod : int { + RTK = 0, ///< Traditional RTK unit hydrograph (default) + AMM = 1 ///< Antecedent Moisture Model +}; + /** * @brief Flow unit systems. * @see Legacy: FlowUnitsType in enums.h @@ -293,6 +304,9 @@ struct SimulationOptions { /** @brief Ignore RDII. */ bool ignore_rdii = false; + /** @brief RDII computation method (RTK or AMM). */ + RdiiMethod rdii_method = RdiiMethod::RTK; + /** @brief Ignore routing. */ bool ignore_routing = false; diff --git a/src/engine/data/InflowData.hpp b/src/engine/data/InflowData.hpp index 5b24b2444..a833eb9f2 100644 --- a/src/engine/data/InflowData.hpp +++ b/src/engine/data/InflowData.hpp @@ -17,6 +17,7 @@ #include #include +#include namespace openswmm { @@ -88,6 +89,7 @@ struct RDIIAssignData { } }; +// ============================================================================ // ============================================================================ // Unit Hydrograph data (from [HYDROGRAPHS] section) // ============================================================================ @@ -122,6 +124,44 @@ struct UnitHydData { void add(const UnitHydEntry& e) { entries.push_back(e); } }; +// ============================================================================ +// AMM assignments (from [AMM] section) +// ============================================================================ + +struct AMMAssignData { + int count() const { return static_cast(node_idx.size()); } + + std::vector node_idx; ///< Target node + std::vector amm_name; ///< AMM component group name + std::vector sewer_area; ///< Tributary sewer area + + void add(int ni, const std::string& name, double area) { + node_idx.push_back(ni); amm_name.push_back(name); + sewer_area.push_back(area); + } + + /// Raw component parameter definition (parsed from [AMM] section). + /// Transferred to AMMSolver::addComponentParams() during init(). + struct ComponentDef { + double area = 0.0; + double RD = 0.0; + double PAT = 0.0; + double HHL = 3600.0; + double AMHL = 28800.0; + double cold_SHCF = 0.0; + double hot_SHCF = 0.0; + double cold_temp = 30.0; + double hot_temp = 70.0; + double TAT = 0.0; + bool is_baseflow = false; + double cold_R = 0.0; + double hot_R = 0.0; + }; + + /// Map: component name → parsed parameters (populated by handle_amm). + std::unordered_map component_defs; +}; + // ============================================================================ // Time patterns (from [PATTERNS] section) // ============================================================================ diff --git a/src/engine/hydrology/AMM.cpp b/src/engine/hydrology/AMM.cpp new file mode 100644 index 000000000..cf4990d68 --- /dev/null +++ b/src/engine/hydrology/AMM.cpp @@ -0,0 +1,285 @@ +/** + * @file AMM.cpp + * @brief Antecedent Moisture Model — engine integration. + * + * @details Equations numbered per Edgren, Czachorski & Gonwa (2024), + * "Reparameterizing the Antecedent Moisture Model", + * J. Water Management Modeling, C525. + * + * @ingroup new_engine + * @license MIT License + */ + +#include "AMM.hpp" +#include "../core/SimulationContext.hpp" +#include +#include +#include + +namespace openswmm { +namespace amm { + +// ============================================================================ +// Public interface +// ============================================================================ + +int AMMSolver::addComponentParams(const std::string& name, const AMMParams& params) { + auto it = name_to_idx_.find(name); + if (it != name_to_idx_.end()) { + params_[static_cast(it->second)] = params; + return it->second; + } + int idx = static_cast(params_.size()); + params_.push_back(params); + name_to_idx_[name] = idx; + return idx; +} + +int AMMSolver::findComponent(const std::string& name) const { + auto it = name_to_idx_.find(name); + if (it == name_to_idx_.end()) return -1; + return it->second; +} + +void AMMSolver::init(SimulationContext& ctx) { + const auto& ad = ctx.amm_assigns; + + // Transfer parsed component definitions into the solver's parameter registry + for (const auto& [name, def] : ad.component_defs) { + AMMParams p; + p.area = def.area; + p.RD = def.RD; + p.PAT = def.PAT; + p.HHL = def.HHL; + p.AMHL = def.AMHL; + p.cold_SHCF = def.cold_SHCF; + p.hot_SHCF = def.hot_SHCF; + p.cold_temp = def.cold_temp; + p.hot_temp = def.hot_temp; + p.TAT = def.TAT; + p.is_baseflow = def.is_baseflow; + p.cold_R = def.cold_R; + p.hot_R = def.hot_R; + addComponentParams(name, p); + } + + int n = ad.count(); + if (n == 0) { + assigns_.clear(); + return; + } + + assigns_.resize(static_cast(n)); + states_.resize(static_cast(n)); + last_SHCF_.assign(static_cast(n), 0.0); + last_RW_.assign(static_cast(n), 0.0); + last_MAP_.assign(static_cast(n), 0.0); + + for (int i = 0; i < n; ++i) { + auto ui = static_cast(i); + assigns_[ui].node_idx = ad.node_idx[ui]; + assigns_[ui].area = ad.sewer_area[ui]; + + int pi = findComponent(ad.amm_name[ui]); + assigns_[ui].param_idx = pi; + + auto& s = states_[ui]; + s.Q_prev = 0.0; + s.RW_prev = 0.0; + + if (pi >= 0 && pi < static_cast(params_.size())) { + const auto& p = params_[static_cast(pi)]; + // Initial ring-buffer sizing (will be resized on first step if dt differs) + int precip_slots = 1; + if (p.PAT > 0.0 && p.HHL > 0.0) { + double est_dt = std::min(p.HHL, p.PAT); + precip_slots = static_cast(p.PAT / est_dt) + 2; + } + s.precip_buf.assign(static_cast(precip_slots), 0.0); + s.precip_head = 0; + + int temp_slots = 1; + if (p.TAT > 0.0 && p.HHL > 0.0) { + double est_dt = std::min(p.HHL, p.TAT); + temp_slots = static_cast(p.TAT / est_dt) + 2; + } + s.temp_buf.assign(static_cast(temp_slots), 0.0); + s.temp_head = 0; + } else { + s.precip_buf.assign(1, 0.0); + s.precip_head = 0; + s.temp_buf.assign(1, 0.0); + s.temp_head = 0; + } + } +} + +void AMMSolver::reset() { + for (auto& s : states_) { + s.Q_prev = 0.0; + s.RW_prev = 0.0; + std::fill(s.precip_buf.begin(), s.precip_buf.end(), 0.0); + s.precip_head = 0; + std::fill(s.temp_buf.begin(), s.temp_buf.end(), 0.0); + s.temp_head = 0; + } + std::fill(last_SHCF_.begin(), last_SHCF_.end(), 0.0); + std::fill(last_RW_.begin(), last_RW_.end(), 0.0); + std::fill(last_MAP_.begin(), last_MAP_.end(), 0.0); +} + +void AMMSolver::computeAll(SimulationContext& ctx, double rainfall, + double air_temp, double dt) { + if (dt <= 0.0) return; + + int n_nodes = ctx.n_nodes(); + + for (size_t c = 0; c < assigns_.size(); ++c) { + int pi = assigns_[c].param_idx; + if (pi < 0 || pi >= static_cast(params_.size())) continue; + + const auto& p = params_[static_cast(pi)]; + auto& s = states_[c]; + + // Use the assigned area, not the parameter area + double area = assigns_[c].area; + + // Ensure ring-buffer sizes match dt + int precip_slots = std::max(static_cast(p.PAT / dt) + 2, 1); + if (static_cast(s.precip_buf.size()) != precip_slots) { + s.precip_buf.assign(static_cast(precip_slots), 0.0); + s.precip_head = 0; + } + + int temp_slots = std::max(static_cast(p.TAT / dt) + 2, 1); + if (static_cast(s.temp_buf.size()) != temp_slots) { + s.temp_buf.assign(static_cast(temp_slots), 0.0); + s.temp_head = 0; + } + + // Eq 3: Moving Average Precipitation + double MAP = updateMovingAverage(s.precip_buf, s.precip_head, rainfall); + + // Eq 11: Moving Average Temperature + double MATemp = updateMovingAverage(s.temp_buf, s.temp_head, air_temp); + + // Eq 2: Shape factor + double SF = shapeFactor(dt, p.HHL); + + double Q_t; + double RW_t; + + if (p.is_baseflow) { + // Baseflow component (2-level, Eq 12-16) + double R_t = computeBaseflowR(MATemp, p.cold_R, p.hot_R, + p.cold_temp, p.hot_temp); + Q_t = computeQ(area, 0.0, R_t, s.RW_prev, + MAP, SF, dt, s.Q_prev, p.area_to_flow); + RW_t = R_t; + last_SHCF_[c] = 0.0; + } else { + // Standard 3-level component (Eq 1-11) + double SHCF = computeSHCF(MATemp, + p.cold_SHCF, p.hot_SHCF, + p.cold_temp, p.hot_temp); + double AMRF = antecedentMoistureRetentionFactor(dt, p.AMHL); + RW_t = computeRW(AMRF, SHCF, MAP, s.RW_prev); + Q_t = computeQ(area, p.RD, RW_t, s.RW_prev, + MAP, SF, dt, s.Q_prev, p.area_to_flow); + last_SHCF_[c] = SHCF; + } + + last_RW_[c] = RW_t; + last_MAP_[c] = MAP; + + // Scatter flow to node lateral inflow + int ni = assigns_[c].node_idx; + if (ni >= 0 && ni < n_nodes) { + ctx.nodes.lat_flow[static_cast(ni)] += Q_t; + } + + s.Q_prev = Q_t; + s.RW_prev = RW_t; + } +} + +// ============================================================================ +// Static helpers — paper equations +// ============================================================================ + +// Eq 2: SF = 0.5^(dt / HHL) +double AMMSolver::shapeFactor(double dt, double HHL) { + if (HHL <= 0.0) return 0.0; + return std::pow(0.5, dt / HHL); +} + +// Eq 6: AMRF = 0.5^(dt / AMHL) +double AMMSolver::antecedentMoistureRetentionFactor(double dt, double AMHL) { + if (AMHL <= 0.0) return 0.0; + return std::pow(0.5, dt / AMHL); +} + +// Eq 3 / Eq 11: Circular-buffer moving average +double AMMSolver::updateMovingAverage(std::vector& buf, int& head, + double new_value) { + if (buf.empty()) return new_value; + buf[static_cast(head % static_cast(buf.size()))] = new_value; + head = (head + 1) % static_cast(buf.size()); + + double sum = 0.0; + for (double v : buf) sum += v; + return sum / static_cast(buf.size()); +} + +// Eq 7-10: SHCF from temperature sigmoid +double AMMSolver::computeSHCF(double MATemp, + double cold_SHCF, double hot_SHCF, + double cold_temp, double hot_temp) { + double L = 1.2 * (cold_SHCF - hot_SHCF); + double temp_range = cold_temp - hot_temp; + if (std::abs(temp_range) < 1.0e-10) return cold_SHCF; + + double k = 4.7964 / temp_range; + double x0 = (cold_temp + hot_temp) / 2.0; + double SHCF = L / (1.0 + std::exp(-k * (MATemp - x0))) + + cold_SHCF - (11.0 / 12.0) * L; + return SHCF; +} + +// Eq 5: RW_t from antecedent moisture +double AMMSolver::computeRW(double AMRF, double SHCF, double MAP, + double RW_prev) { + if (AMRF <= 0.0 || AMRF >= 1.0 - 1.0e-15) { + return SHCF * MAP + RW_prev; + } + double correction = (AMRF - 1.0) / std::log(AMRF); + return std::max(correction * SHCF * MAP + AMRF * RW_prev, 0.0); +} + +// Eq 13: Baseflow capture fraction +double AMMSolver::computeBaseflowR(double MATemp, + double cold_R, double hot_R, + double cold_temp, double hot_temp) { + double L = 1.2 * (cold_R - hot_R); + double temp_range = cold_temp - hot_temp; + if (std::abs(temp_range) < 1.0e-10) return cold_R; + + double k = 4.7964 / temp_range; + double x0 = (cold_temp + hot_temp) / 2.0; + double R = L / (1.0 + std::exp(-k * (MATemp - x0))) + + cold_R - (11.0 / 12.0) * L; + return std::max(R, 0.0); +} + +// Eq 1: Q_t = A * (RD + (RW_t + RW_prev)/2) * MAP * (1-SF)/dt + SF * Q_prev +double AMMSolver::computeQ(double area, double RD, double RW_t, double RW_prev, + double MAP, double SF, double dt, + double Q_prev, double area_to_flow) { + if (dt <= 0.0) return 0.0; + double capture = std::max(RD + (RW_t + RW_prev) / 2.0, 0.0); + double Q_new = area * area_to_flow * capture * MAP * (1.0 - SF) / dt; + return std::max(Q_new + SF * Q_prev, 0.0); +} + +} // namespace amm +} // namespace openswmm diff --git a/src/engine/hydrology/AMM.hpp b/src/engine/hydrology/AMM.hpp new file mode 100644 index 000000000..7e88effc0 --- /dev/null +++ b/src/engine/hydrology/AMM.hpp @@ -0,0 +1,189 @@ +/** + * @file AMM.hpp + * @brief Antecedent Moisture Model — alternative RDII formulation. + * + * @details Implements the reparameterized AMM equations from: + * + * Edgren, D., Czachorski, R., & Gonwa, W. (2024). + * "Reparameterizing the Antecedent Moisture Model." + * Journal of Water Management Modeling, C525. + * https://doi.org/10.14796/JWMM.C525 + * + * The AMM replaces the RTK unit-hydrograph convolution with a three-level + * recursive model that accounts for antecedent soil moisture and seasonal + * temperature effects on rainfall capture fraction: + * + * Level 1 — Rainfall-runoff function (Eq 1-4) + * Level 2 — Antecedent moisture (Eq 5-6) + * Level 3 — Seasonal hydrologic cond. (Eq 7-11) + * + * Plus an optional baseflow component (Eq 12-16). + * + * @note Selected via RDII_METHOD AMM in [OPTIONS]. + * @ingroup new_engine + * + * @license MIT License + */ + +#ifndef OPENSWMM_AMM_HPP +#define OPENSWMM_AMM_HPP + +#include +#include +#include +#include + +namespace openswmm { + +struct SimulationContext; + +namespace amm { + +// ============================================================================ +// Parameter structure (one per component in [AMM] section) +// ============================================================================ + +/** + * @brief Parameters for one AMM component (standard 3-level or 2-level baseflow). + * + * ### Calibratable parameters (standard component) + * | Symbol | Field | Description | Units | + * |-----------|-------------|------------------------------------------|---------| + * | RD | RD | Dry-weather capture fraction | - | + * | PAT | PAT | Precipitation averaging time | seconds | + * | HHL | HHL | Hydrograph half-life | seconds | + * | AMHL | AMHL | Antecedent moisture half-life | seconds | + * | Cold SHCF | cold_SHCF | Cold-weather SHCF (high capture) | 1/L | + * | Hot SHCF | hot_SHCF | Hot-weather SHCF (low capture) | 1/L | + * | Cold Temp | cold_temp | Cold reference temperature | deg | + * | Hot Temp | hot_temp | Hot reference temperature | deg | + * | TAT | TAT | Temperature averaging time | seconds | + */ +struct AMMParams { + double area = 0.0; ///< Contributing area (project area units) + double RD = 0.0; ///< Min capture fraction (dry weather) + double PAT = 0.0; ///< Precipitation averaging time (sec) + double HHL = 3600.0; ///< Hydrograph half-life (sec) + double AMHL = 28800.0;///< Antecedent moisture half-life (sec) + double cold_SHCF = 0.0; ///< SHCF at cold temperature (1/length) + double hot_SHCF = 0.0; ///< SHCF at hot temperature (1/length) + double cold_temp = 30.0; ///< Cold reference temp + double hot_temp = 70.0; ///< Hot reference temp + double TAT = 0.0; ///< Temperature averaging time (sec) + + bool is_baseflow = false;///< True -> skip Level 2, use cold_R / hot_R + double cold_R = 0.0; ///< Cold capture fraction (baseflow) + double hot_R = 0.0; ///< Hot capture fraction (baseflow) + + double area_to_flow = 1.0 / 12.0; ///< area*depth -> volume conversion +}; + +// ============================================================================ +// Per-component runtime state +// ============================================================================ + +struct AMMState { + double Q_prev = 0.0; + double RW_prev = 0.0; + + std::vector precip_buf; + int precip_head = 0; + + std::vector temp_buf; + int temp_head = 0; +}; + +// ============================================================================ +// Solver +// ============================================================================ + +/** + * @brief Antecedent Moisture Model solver for engine integration. + * + * Manages multiple AMM components (each assigned to a node) and computes + * RDII inflow contributions at each timestep via the step() interface + * or the computeAll() interface that matches RDIISolver's pattern. + */ +class AMMSolver { +public: + /** + * @brief Add a named AMM component group. + * @param name Group name (matches [AMM] section 2nd field). + * @param params AMM parameters for this component. + * @returns Index of the registered parameter set. + */ + int addComponentParams(const std::string& name, const AMMParams& params); + + /** + * @brief Look up AMM component index by name. + * @returns Index, or -1 if not found. + */ + int findComponent(const std::string& name) const; + + /** + * @brief Initialize from SimulationContext AMM assignments. + * + * Resolves named AMM groups to parameter indices and allocates + * ring buffers for the moving-average computations. + */ + void init(SimulationContext& ctx); + + /** + * @brief Reset all state to zero (for re-running). + */ + void reset(); + + /** + * @brief Compute AMM inflows for all groups, scatter to node lat_flow. + * + * @param ctx SimulationContext (reads nodes.lat_flow, air_temp). + * @param rainfall Average rainfall intensity this step (depth/time). + * @param air_temp Air temperature this step (same units as params). + * @param dt Routing timestep (seconds). + */ + void computeAll(SimulationContext& ctx, double rainfall, + double air_temp, double dt); + + int componentCount() const { return static_cast(params_.size()); } + +private: + // Named parameter registry (from [AMM] component definitions) + std::vector params_; + std::unordered_map name_to_idx_; + + // Per-assignment runtime data + struct Assignment { + int node_idx = -1; + int param_idx = -1; + double area = 0.0; + }; + std::vector assigns_; + std::vector states_; + + // Diagnostic snapshots + std::vector last_SHCF_; + std::vector last_RW_; + std::vector last_MAP_; + + // ---- Helper functions matching paper equations ---- + static double shapeFactor(double dt, double HHL); + static double antecedentMoistureRetentionFactor(double dt, double AMHL); + static double updateMovingAverage(std::vector& buf, int& head, + double new_value); + static double computeSHCF(double MATemp, + double cold_SHCF, double hot_SHCF, + double cold_temp, double hot_temp); + static double computeRW(double AMRF, double SHCF, double MAP, + double RW_prev); + static double computeBaseflowR(double MATemp, + double cold_R, double hot_R, + double cold_temp, double hot_temp); + static double computeQ(double area, double RD, double RW_t, double RW_prev, + double MAP, double SF, double dt, + double Q_prev, double area_to_flow); +}; + +} // namespace amm +} // namespace openswmm + +#endif // OPENSWMM_AMM_HPP diff --git a/src/engine/input/handlers/InflowsHandler.cpp b/src/engine/input/handlers/InflowsHandler.cpp index 4f853edf8..2b77eb06f 100644 --- a/src/engine/input/handlers/InflowsHandler.cpp +++ b/src/engine/input/handlers/InflowsHandler.cpp @@ -46,6 +46,7 @@ #include #include #include +#include namespace openswmm::input { @@ -212,7 +213,7 @@ void handle_hydrographs(SimulationContext& ctx, const std::vector& // Otherwise: UHgroup Month Response R T K [Dmax Drecov Dinit] if (tok.size() < 6) continue; - // Parse month: "All" → -1, or numeric 1-12 → 0-based (0-11) + // Parse month: "All" -> -1, or numeric 1-12 -> 0-based (0-11) int month = -1; std::string month_upper = Tokenizer::to_upper(tok[1]); if (month_upper != "ALL") { @@ -243,4 +244,91 @@ void handle_hydrographs(SimulationContext& ctx, const std::vector& } } +// ============================================================================ +// handle_amm() +// ============================================================================ +/** + * @brief Parse the [AMM] section. + * + * Two line formats are accepted: + * + * ### Assignment line (3 tokens) — assigns an AMM component to a node + * ``` + * ;; Node AMMgroup SewerArea + * J1 Fast 1000.0 + * ``` + * + * ### Parameter line (first token == component name, second token == keyword) + * ``` + * ;; CompName Keyword Values... + * Fast PARAMS 0.10 1800 3600 28800 0.05 0.01 35.0 75.0 604800 + * Fast BASEFLOW 0.10 1800 3600 0.03 0.01 35.0 75.0 604800 + * Fast AREA 5.0 + * ``` + * + * PARAMS keyword fields (9): RD PAT HHL AMHL coldSHCF hotSHCF coldTemp hotTemp TAT + * BASEFLOW keyword fields (8): RD PAT HHL coldR hotR coldTemp hotTemp TAT + * AREA keyword fields (1): ContributingArea + */ +void handle_amm(SimulationContext& ctx, const std::vector& lines) { + // Two-pass: first collect component parameter definitions, then assignments. + + // Pass 1: scan all lines for PARAMS / BASEFLOW / AREA keywords + for (const auto& line : lines) { + auto tok = Tokenizer::tokenize(line); + if (tok.size() < 2) continue; + + std::string kw = tok[1]; + std::transform(kw.begin(), kw.end(), kw.begin(), ::toupper); + + if (kw == "PARAMS" && tok.size() >= 11) { + auto& p = ctx.amm_assigns.component_defs[tok[0]]; + p.RD = to_double(tok[2]); + p.PAT = to_double(tok[3]); + p.HHL = to_double(tok[4]); + p.AMHL = to_double(tok[5]); + p.cold_SHCF = to_double(tok[6]); + p.hot_SHCF = to_double(tok[7]); + p.cold_temp = to_double(tok[8]); + p.hot_temp = to_double(tok[9]); + p.TAT = to_double(tok[10]); + p.is_baseflow = false; + } + else if (kw == "BASEFLOW" && tok.size() >= 10) { + auto& p = ctx.amm_assigns.component_defs[tok[0]]; + p.RD = to_double(tok[2]); + p.PAT = to_double(tok[3]); + p.HHL = to_double(tok[4]); + p.cold_R = to_double(tok[5]); + p.hot_R = to_double(tok[6]); + p.cold_temp = to_double(tok[7]); + p.hot_temp = to_double(tok[8]); + p.TAT = to_double(tok[9]); + p.is_baseflow = true; + } + else if (kw == "AREA" && tok.size() >= 3) { + ctx.amm_assigns.component_defs[tok[0]].area = to_double(tok[2]); + } + } + + // Pass 2: scan for assignment lines (Node CompName SewerArea) + for (const auto& line : lines) { + auto tok = Tokenizer::tokenize(line); + if (tok.size() < 3) continue; + + // Skip parameter definition lines (second token is a keyword) + std::string kw = tok[1]; + std::transform(kw.begin(), kw.end(), kw.begin(), ::toupper); + if (kw == "PARAMS" || kw == "BASEFLOW" || kw == "AREA") continue; + + const int node_idx = ctx.node_names.find(tok[0]); + if (node_idx < 0) continue; + + const std::string& amm_name = tok[1]; + double sewer_area = to_double(tok[2]); + + ctx.amm_assigns.add(node_idx, amm_name, sewer_area); + } +} + } /* namespace openswmm::input */ diff --git a/src/engine/input/handlers/InflowsHandler.hpp b/src/engine/input/handlers/InflowsHandler.hpp index 35cf1a511..9ae490633 100644 --- a/src/engine/input/handlers/InflowsHandler.hpp +++ b/src/engine/input/handlers/InflowsHandler.hpp @@ -33,6 +33,9 @@ void handle_rdii(SimulationContext& ctx, const std::vector& lines); /** @brief Parse [HYDROGRAPHS] into UnitHydData. */ void handle_hydrographs(SimulationContext& ctx, const std::vector& lines); +/** @brief Parse [AMM] into AMMAssignData and register AMM component parameters. */ +void handle_amm(SimulationContext& ctx, const std::vector& lines); + } /* namespace openswmm::input */ #endif /* OPENSWMM_ENGINE_INFLOWS_HANDLER_HPP */ diff --git a/src/engine/input/handlers/OptionsHandler.cpp b/src/engine/input/handlers/OptionsHandler.cpp index a60ff147e..d01a87fd0 100644 --- a/src/engine/input/handlers/OptionsHandler.cpp +++ b/src/engine/input/handlers/OptionsHandler.cpp @@ -221,6 +221,13 @@ void handle_options(SimulationContext& ctx, const std::vector& line opt.ignore_groundwater = Tokenizer::parse_boolean(val); } else if (key == "IGNORE_RDII") { opt.ignore_rdii = Tokenizer::parse_boolean(val); + } else if (key == "RDII_METHOD") { + std::string uval = val; + std::transform(uval.begin(), uval.end(), uval.begin(), ::toupper); + if (uval == "AMM") + opt.rdii_method = RdiiMethod::AMM; + else + opt.rdii_method = RdiiMethod::RTK; } else if (key == "IGNORE_ROUTING") { opt.ignore_routing = Tokenizer::parse_boolean(val); } else if (key == "IGNORE_QUALITY") { diff --git a/src/engine/plugins/DefaultInputPlugin.cpp b/src/engine/plugins/DefaultInputPlugin.cpp index c7c51b03f..ccd288e0d 100644 --- a/src/engine/plugins/DefaultInputPlugin.cpp +++ b/src/engine/plugins/DefaultInputPlugin.cpp @@ -93,6 +93,7 @@ void DefaultInputPlugin::register_builtin_handlers() { registry_.register_builtin("DWF", input::handle_dwf); registry_.register_builtin("RDII", input::handle_rdii); registry_.register_builtin("HYDROGRAPHS", input::handle_hydrographs); + registry_.register_builtin("AMM", input::handle_amm); registry_.register_builtin("LOADINGS", input::handle_loadings); registry_.register_builtin("PATTERNS", input::handle_patterns);