SolverOptions2D.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_SOLVER_OPTIONS_HPP
#define OPENSWMM_ENGINE_2D_SOLVER_OPTIONS_HPP
 
#include <cstdint>
#include <cmath>
#include <string>
#include <vector>
 
#include "../../core/FilePathPair.hpp"
 
namespace openswmm::twoD {

enum class CellClosure2D : int8_t {
    FLAT = 0, 
    VFR = 1   
};

enum class FaceDepth2D : int8_t {
    MEAN = 0,    
    VFR_FACE = 1 
};

enum class RainfallMode : int8_t {
    NATURAL_NEIGHBOUR = 0, 
    SYSTEM = 1,            
    NONE = 2,              
    NEAREST_NEIGHBOUR = 3  
};

enum class CouplingLaw2D : int8_t { ORIFICE = 0, WEIR_ORIFICE = 1 };

enum class CouplingHead2D : int8_t { CELL = 0, VERTEX = 1, STENCIL_WET_MEAN = 2 };

enum class CouplingPondMode2D : int8_t { FOOTPRINT = 0, CHAMBER = 1, TRACKING = 2 };

enum class Backend2D : int8_t {
    CPU = 0,  
    AUTO = 1, 
    OMP = 2,  
    CUDA = 3, 
    HIP = 4,  
    SYCL = 5  
};

enum class Momentum2D : int8_t { LOCAL_INERTIAL = 0, FULL_SWE = 1, DIFFUSIVE_WAVE = 2 };

enum class OutputPrecision2D : int8_t {
    FLOAT32 = 0, 
    FLOAT64 = 1  
};

namespace report2d {
enum Var : unsigned {
    DEPTH = 1u << 0,        
    VELOCITY = 1u << 1,     
    EDGE_FLUX = 1u << 2,    
    NODE_HEAD = 1u << 3,    
    SPECIES = 1u << 4,      
    RAINFALL = 1u << 5,     
    INFILTRATION = 1u << 6, 
    COUPLING = 1u << 7,     
    GRADIENTS = 1u << 8,    
    CONTINUITY = 1u << 9,   
    ENVELOPES = 1u << 10,   
    BUILDUP = 1u << 11,     
    GROUNDWATER =
        1u
        << 12, 
    GW_DETAILED =
        1u << 13,    
    HEAT = 1u << 14, 
    TURBULENCE = 1u << 15, 
    ALL_MASK = (1u << 16) - 1u,
    DEFAULT_MASK = DEPTH | VELOCITY | EDGE_FLUX | NODE_HEAD | SPECIES | RAINFALL | INFILTRATION |
                   ENVELOPES | BUILDUP | GROUNDWATER,
    MINIMAL_MASK = DEPTH | NODE_HEAD | ENVELOPES,
};
} // namespace report2d

enum class Coriolis2D : int8_t { OFF = 0, LATITUDE = 1, CONSTANT_F = 2 };
enum class Turbulence2D : int8_t { NONE = 0, SMAGORINSKY = 1, MIXING_LENGTH = 2 };
inline const char* turbulenceName(Turbulence2D t) noexcept {
    return t == Turbulence2D::SMAGORINSKY ? "SMAGORINSKY" : t == Turbulence2D::MIXING_LENGTH ? "MIXING_LENGTH" : "NONE";
}
 
struct SolverOptions2D {
    // Explicit opt-in. Unauthored fields are omitted by the writer.
    bool coupled_heat = false, heat_model_authored = false;
    int heat_surface = -1, heat_radiative = -1; // AUTO inherits network module
    bool heat_surface_authored = false, heat_radiative_authored = false;
    int bed_heat_model = 0; // NONE, LOCAL_STORAGE, GROUND_COLUMN
    bool bed_heat_authored = false;
    // Derived during initialization: legacy runtime-water conversion stays
    // 4.186 MJ/(m3 K); opt-in coupled heat uses configured water rho*cp.
    double water_heat_capacity = 4186000.0;
    Momentum2D momentum = Momentum2D::LOCAL_INERTIAL;
    Coriolis2D coriolis = Coriolis2D::OFF;
    Turbulence2D turbulence = Turbulence2D::NONE;
    double coriolis_latitude = 0.0, coriolis_f = 0.0, smagorinsky_cs = 0.10;
    double eddy_viscosity_background = 0.0;
    double mixing_length_coefficient = 0.5;
    bool coriolis_latitude_set = false, coriolis_f_set = false;
    bool wall_drag_enabled = false;
    bool momentumPhysics() const noexcept {
        return wall_drag_enabled || coriolis != Coriolis2D::OFF || turbulence != Turbulence2D::NONE;
    }
    double coriolisParameter() const noexcept {
        return coriolis == Coriolis2D::CONSTANT_F ? coriolis_f
             : coriolis == Coriolis2D::LATITUDE
                 ? 1.4584231711061317e-4 * std::sin(coriolis_latitude * 0.017453292519943295)
                 : 0.0;
    }

    int front_rebuild = -1; 
    int reconstruction_order = 1;
    double max_timestep = 10.0;
    double dry_depth = 0.001;        
    double limiter_epsilon = 1.0e-6; 
    double flux_dh_eps = 0.004; 
    double coupling_sync = 0.0;
 
    double coupling_cd = 0.65; 
    CouplingLaw2D coupling_law = CouplingLaw2D::ORIFICE;
    CouplingHead2D coupling_head = CouplingHead2D::CELL;
    int coupling_coeffs = 0;
    CouplingPondMode2D coupling_pond_mode = CouplingPondMode2D::FOOTPRINT;
    double coupling_pond_relax = 0.0;
    double coupling_rim_band = -1.0;
    int8_t coupling_equalize = -1;
    bool report_2d = true; 
 
    // ---- Results-file size controls (plans/2D_OUTPUT_FLOAT32_AND_VARIABLE_SELECTION_PLAN) ----
    OutputPrecision2D output_precision = OutputPrecision2D::FLOAT32;
    int output_compression = 4;
    unsigned report_2d_vars = report2d::DEFAULT_MASK;
    std::vector<std::string> report_2d_species;
    double report_2d_step = 0.0;
 
    // -----------------------------------------------------------------------
    // E2 (2026-09-07) — process enables. Every default reproduces the
    // behaviour before the keys existed, so a deck without them is
    // bit-identical. Consumed in SurfaceRouter2D::initialize() (infiltration
    // rows, transport row layout via transport::surface2DEnables) and the
    // per-step forcing refresh (evaporation).
    // -----------------------------------------------------------------------
    int8_t infiltration = -1;
    double infil_step = 0.0;
    std::string infil_default_method;
    std::string infil_destination;
    int8_t evaporation = 1;
    bool transport_pollutants = true;
    bool transport_msx = true;
    bool transport_age = true;
    bool transport_temperature = true;
 
    // -----------------------------------------------------------------------
    // U5 (2026-09-07), rewired 2026-09-07 once the G1 two-zone kernel landed.
    // These are the PROCESS ENABLES for the integrated 2D subsurface; the
    // parameters themselves live in [2D_AQUIFER_OPTIONS] / [2D_AQUIFER] /
    // [2D_AQUIFER_NODE]. Same shape as INFILTRATION / INFIL_STEP above: a
    // tri-state enable whose AUTO reproduces the pre-existing behaviour, and
    // an alias that folds into the authoritative struct at open.
    // -----------------------------------------------------------------------
    int8_t groundwater = -1;
    std::string gw_et;
 
    // Rainfall→mesh mapping. Default NATURAL_NEIGHBOUR (spatial interpolation
    // across all located gages); SYSTEM applies the uniform all-gage mean.
    // Parsed from [2D_OPTIONS] RAINFALL_MODE; env OPENSWMM_2D_RAINFALL_MODE
    // (natural|nearest|system|none) overrides.
    RainfallMode rainfall_mode = RainfallMode::NATURAL_NEIGHBOUR;
 
    // Volume → free-surface cell closure. Default FLAT (legacy η = tri_cz + V/A):
    // fast and the right choice for typical deep-water urban flooding, where a
    // partially wet cell is the exception. VFR (planar-bed volume/free-surface,
    // fully implemented on all backends) restores the C-property at shorelines
    // and removes the water-climbs-uphill artifact, but resolves the shoreline
    // wetting/drying FLAT freezes out — ~3–8× more marcher substeps — so it is
    // OPT-IN (best on shallow water / gentle slopes). Parsed from [2D_OPTIONS]
    // CELL_CLOSURE (FLAT|VFR). See plans/2d/2D_VFR_SOLVER_CLOSURE_PLAN.md.
    CellClosure2D cell_closure = CellClosure2D::FLAT;
 
    // Effective conveyance depth at shared edges. Default MEAN (legacy upwind
    // cell-mean depth). VFR_FACE (B&S Eq. 14 face depth + wetting gate) pairs
    // with CELL_CLOSURE=VFR to complete the artifact fix; opt-in for the same
    // reason. Parsed from [2D_OPTIONS] FACE_RECONSTRUCTION (MEAN|VFR_FACE).
    FaceDepth2D face_reconstruction = FaceDepth2D::MEAN;

    double vfr_min_wet_frac = 0.01;
 
    // -----------------------------------------------------------------------
    // Explicit local-inertial marcher (INTEGRATOR EXPLICIT — the only 2D
    // integrator since the D2 CVODE/ARKODE retirement) options. Defaults per
    // the 2026-07-29 reimplementation plan.
    // -----------------------------------------------------------------------
    double theta = 0.8;      
    double theta_reference_step = 0.0; 
    double cfl_number = 0.7; 
    double h_move = 0.003;   
    int lts_tiers = 4;       
    double froude_max = 1.5; 
    bool advection = false; 
    double exchange_beta = 0.8;

    double dispersion = 0.0;
    double exchange_relax = 1.0;
    bool coupling_area_auto = false;

    bool coupling_in_flooding = false;

    Backend2D backend = Backend2D::AUTO;

    bool auto_backend_cpu = false;

    openswmm::FilePathPair mesh_file;

    openswmm::FilePathPair output_file;
 
    // -----------------------------------------------------------------------
    // Unit-system coupling factors — NOT parsed from input. Computed once in
    // SurfaceRouter2D::initialize() from the project FLOW_UNITS.
    //
    // The 2D solver runs internally in SI (metres, m², m³, m³/s, g=9.80665).
    // The 1D SWMM engine ALWAYS computes internally in FEET (g=32.2, PHI=1.486)
    // — for EVERY project, US or SI: its reader converts metric inputs to feet
    // on load and only converts back at the display/output boundary. So these
    // coupling factors are ALWAYS the feet⇄metres conversion, independent of
    // FLOW_UNITS. SurfaceRouter2D::initialize() overwrites the 1.0 defaults
    // with the real ft⇄m factors; the defaults only stand when 2D is inactive
    // (no coupling occurs). The MESH scaling factor is separate and IS driven
    // by FLOW_UNITS (the mesh is authored in project units) — see initialize().
    // -----------------------------------------------------------------------
    double len_1d_to_2d = 1.0;  
    double len_2d_to_1d = 1.0;  
    double vol_1d_to_2d = 1.0;  
    double flow_1d_to_2d = 1.0; 
    double flow_2d_to_1d = 1.0; 

    int num_threads = 1;

    int requested_threads = 0;

    bool mesh_units_si = false;

    bool mesh_scaled_to_si = false;

    double mesh_to_si_factor = 1.0;

    bool pending_rows_drained = false;

    double bc_flow_to_si_applied = 1.0;
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_SOLVER_OPTIONS_HPP