PySWMM Documentation which is a Python package the allows seamless interaction with the USEPA-SWMM5

Source of this file is https://pypi.python.org/pypi/pyswmm/

Start here to begin working with pyswmm.

The pyswmm package allows seamless interaction with the USEPA-SWMM5
data model.  Parameter getters/setters and results getters have been
exposed, allowing the user to see results while a simulation is
running as well as update link settings.

Loading a Model
---------------

There are three options to load a model. If there is no desire to
interact with the simulation then the simplest way to run the
model is the following:

.. code-block:: python

>>> from pyswmm import Simulation
>>>
>>> sim = Simulation('./testmodel.inp')
>>> sim.execute()


The following method allows the user to read in a model and
manually step through the simulation and get/set parameters and
results.  This scenario is the cleanest solution using a
with statement. It automatically cleans up after the
simulation is complete. 

.. code-block:: python

>>> from pyswmm import Simulation
>>>
>>> with Simulation('./testmodel.inp') as sim:
... for step in sim:
... pass
... sim.report()

One feature that pyswmm adds to the modeling world is the simulation
stride function ``step_advance``.  Assuming a user has developed all
of their control rules in in a Python script, to reduce simulation
time a user can specify how often Python controls should be evaluated. 

For example, let's assume ``testmodel.inp`` has a 30 second routing step
(using the dynamic wave solver, this step could vary significantly).  If
complex control scenarios are developed, evaluating rules could add
significant time to the simulation.

.. code-block:: python

>>> from pyswmm import Simulation
>>>
>>> with Simulation('testmodel.inp') as sim:
... sim.step_advance(300)
... for step in sim:
... print(sim.current_time)
... # or here! sim.step_advance(newvalue)
... sim.report()

2015-11-01 14:05:00
2015-11-01 14:10:00
2015-11-01 14:15:00
2015-11-01 14:20:00

Nodes
-----

For interacting with nodes a :py:class:`pyswmm.nodes.Nodes` object must be initialized.
See the following example. Once the ``Nodes`` object is initialized,
you can then initialize a :py:class:`pyswmm.nodes.Node`

.. code-block:: python

>>> from pyswmm import Simulation, Nodes
>>>
>>> with Simulation('./testmodel.inp') as sim:
... node_object = Nodes(sim)
...
... #J1 node instantiation
... J1 = node_object["J1"]
... print(J1.invert_elevation)
... print(J1.is_junction())
...
... #Step through a simulation
... for step in sim:
... print(J1.total_inflow)
...
... sim.report()


Links
-----

For interacting with nodes a :py:class:`pyswmm.links.Links` object must be initialized.
See the following example. Once the ``Links`` object is initialized,
you can then initialize a :py:class:`pyswmm.links.Link`

.. code-block:: python


>>> from pyswmm import Simulation, Links
>>>
>>> with Simulation('./testmodel.inp') as sim:
... link_object = Links(sim)
...
... #C1:C2 link instantiation
... c1c2 = link_object["C1:C2"]
... print(c1c2.flow_limit)
... print(c1c2.is_conduit())
...
... #Step through a simulation
... for step in sim:
... print(c1c2.flow)
... if c1c2.flow > 10.0:
... c1c2.target_setting = 0.5
...
... sim.report()


Subcatchments
-------------

For interacting with subcatchments a :py:class:`pyswmm.subcatchments.Subcatchments`
object must be initialized. See the following example. Once the ``Subcatchments`` object is initialized,
you can then initialize a :py:class:`pyswmm.subcatchments.Subcatchment`

.. code-block:: python


>>> from pyswmm import Simulation, Subcatchments
>>>
>>> with Simulation('./testmodel.inp') as sim:
... subcatch_object = Subcatchments(sim)
...
... #SC1 subcatchment instantiation
... SC1 = subcatch_object["S1"]
... print(SC1.area)
...
... #Step through a simulation
... for step in sim:
... print(SC1.runoff)
...
... sim.report()


In the example above we introduce the option to change a link's settings.

PySWMM Controls
---------------

The pyswmm package exposes new possibility in interfacing with models.  All control
rules can now be removed from USEPA SWMM5 and brought into Python.  Now that this
functionality exists, open-source Python packages can now be used in conjunction
with pyswmm to bring even more complex control routines. 

The following example illustrates the use of functions for
comparing two depths.

.. code-block:: python

>>> from pyswmm import Simulation, Links, Nodes
>>>
>>> def TestDepth(node, node2):
>>> if node > node2:
>>> return True
>>> else:
>>> return False
>>>
>>> with Simulation('./testmodel.inp') as sim:
... link_object = Links(sim)
...
... #C1:C2 link instantiation
... c1c2 = link_object["C1:C2"]
...
... node_object = Nodes(sim)
... #J1 node instantiation
... J1 = node_object["J1"]
... #J2 node instantiation
... J2 = node_object["J2"]
...
... #Step through a simulation
... for step in sim:
... if TestDepth(J1.depth, J2.depth):
... c1c2.target_setting = 0.5
...
... sim.report()

If an EPA-SWMM5 Model has existing control actions within, any control
rules developed using pyswmm will have the highest priority.  All pyswmm
control actions are evaluated at the end of each simulation step, after
EPA-SWMM native controls have been evaluated.  If control actions are reported,
any control action updated by pyswmm will be output to the *.rpt file.


Generate Node Inflows
---------------------

Among the newest features pyswmm brings to SWMM5 modeling is the ability to
set a nodes inflow.  This can enable the user to model different behavior such as
runoff or seasonality. 

.. code-block:: python

>>> from pyswmm import Simulation, Nodes
>>>
>>> with Simulation('/testmodel.inp') as sim:
... j1 = Nodes(sim)["J1"]
... for step in sim:
... j1.generated_inflow(9)

How to Set your Model Preferences in InfoSWMM_SA

InfoSWMM SA Preferences

The Operation tab allows the user to control project settings relating to the management of various InfoSWMM SA operations. It should be noted that all the changes made on the Preferences dialog box must be saved prior to closing the dialog box. Many preferences set with this command will be reflected as the default choices on other dialog boxes.  You may change those settings as desired on those dialog boxes.

Contents of the Operation Settings tab of the preferences dialog is described below:

NAME	DESCRIPTION
Auto Length Calculation
Auto Length/Area Calculation
Auto Record Saving
Auto Database Packing
Delete Confirmation
Auto Link Node Inclusion
Check Update / Upgrade
Batch Import on Open
Batch Export on Save
Allow Duplicated Report
Store Absolute Conduit Invert
Allow Project Database Editing Buffer
Edit Diameter of (Filled) Circular Pipe in inch/mm in Attribute Browser
Include Extra Summary Data for Map Display
Use the Simulation Task Manager
Display Each Help Item in its Own Window	On some computers you may need to select this option for the online help to display properly
Text Editor
Length Scaling Factor
Area Scaling Factor
Auto Output Relate Update
Auto Output Retrieval
Auto Output Remembering
Enable Output Save As
Auto Link Delete
Single Output Report Loading
Display Calibration Data in Graph
Sewer Interface
Store Absolute Junction Rim

The Display tab controls visual effects and features within InfoSWMM SA. All changes made on the Preferences dialog box must be saved prior to closing the dialog box.  Many preferences set with this command will be reflected as the default choices on other dialog boxes.  You may change those settings as desired on those dialog boxes.

Contents of the Display Settings tab of the preferences dialog is described below:

NAME	DESCRIPTION
Node Locate Map Extent %
Pipe Locate Zoom Factor
Decimal Placement
Roughness Decimal Placement
Default Domain Highlighting Color
Selection Highlighting Color
Default Inactive Element Color
Default Google Maps Link	Google Maps can be launched from the Attribute Browser.  Choose the default view that the map will be opened in.
Height of Note Editing Box
Show Inactive Elements at a Default Display
Show Link Direction
Arrow Symbol Size
Show Domain With Color Coding
Show Inactive Elements with Color Coding
Show Subcatchment Linkage

The Preferences - Selection Settings command is used to set the user's preferences during selection of elements in InfoSWMM SA project using the available selection tools. All changes made on the Preferences dialog box must be saved (by clicking the OK button) prior to closing the dialog box.  Many preferences set with this command will be reflected as the default choices on other dialog boxes.  You may change those settings as desired on those dialog boxes.

A brief explanation of the elements of the contents:

NAME	DESCRIPTION
Selection Shape
Selection Mode
Selection Tolerance	This setting describes how many pixels away from the mouse pointer will be evaluated for the existence of map elements.  This setting has more significance when you are trying to single-click on a map element instead of trying to select it with the Window or Crossing methods.

The Default Symbols Size tab is used to adjust the size (graphical appearance) of the InfoSWMM SA data elements listed in the dialog editor shown below. All changes made on the Preferences dialog box must be saved prior to closing the dialog box. Many preferences set with this command will be reflected as the default choices on other dialog boxes. You may change those settings as desired on those dialog boxes.

The ID and Description tab allows the user to change InfoSWMM SA defaults for data element naming, data sets, curves, patterns, etc. All changes made on the Preferences dialog box must be saved prior to closing the dialog box. Many preferences set with this command will be reflected as the default choices on other dialog boxes.  You may change those settings as desired on those dialog boxes.

Contents of the ID and Description tab of the preferences dialog is described below:

NAME	DESCRIPTION
ID Suggestion
Data Type
ID Prefix
Next
Increment
Description

This section of the Preferences dialog box provides you with the tools to customize InfoSWMM's SA label settings. All changes made on the Preferences dialog box must be saved prior to closing the dialog box. Many preferences set with this command will be reflected as the default choices on other dialog boxes.  You may change those settings as desired on those dialog boxes.

Contents of the Labeling tab of the preferences dialog is described below:

NAME	DESCRIPTION
Node/Link/Subcatchment Label Settings	Enables specifying different node/link/Subcatchment label settings. Symbol - This launches the Symbol Selector dialog box. Specify the font, color, font size, font color, font fill style and other advanced features here. All the options that you specify here will be used while labelling the nodes/links/Subcatchments. Placement - Use this to specify the label placement. This option is not available for Subcatchments. Text Style Previewer - Once you specify the label text, a preview is displayed in this section. If you do not like the text properties you may change it by clicking on the Symbol button.
Label Node/Link/Subcatchment With IDs	Use this button to label all your nodes/links/Subcatchments.
Clear Node/Link/Subcatchment With IDs	Use this to clear your InfoSWMM SA map node/link/Subcatchment labels by clicking on  this button.

Robert Dickinson on Earn.Com

My http://Earn.COM  referral link. You can sign up to Earn Bitcoin from emails sent to your email address http://earn.com/dickinsonre/referral/?a=eq153i1b072smu3t …

The beautiful code for the Water Quality Treatment Functions in #SWMM5

//-----------------------------------------------------------------------------

//   treatmnt.c

//

//   Project:  EPA SWMM5

//   Version:  5.1

//   Date:     03/20/14   (Build 5.1.001)

//             03/19/15   (Build 5.1.008)

//   Author:   L. Rossman

//

//   Pollutant treatment functions.

//

//   Build 5.1.008:

//   - A bug in evaluating recursive calls to treatment functions was fixed.

//

//-----------------------------------------------------------------------------

#define _CRT_SECURE_NO_DEPRECATE

#include <stdlib.h>

#include <string.h>

#include "headers.h"

//-----------------------------------------------------------------------------

//  Constants

//-----------------------------------------------------------------------------

static const int PVMAX = 5;            // number of process variables

enum   ProcessVarType {pvHRT,          // hydraulic residence time

                       pvDT,           // time step duration

                       pvFLOW,         // flow rate

                       pvDEPTH,        // water height above invert

                       pvAREA};        // storage surface area

//-----------------------------------------------------------------------------

//  Shared variables

//-----------------------------------------------------------------------------

static int     ErrCode;                // treatment error code

static int     J;                      // index of node being analyzed

static double  Dt;                     // curent time step (sec)

static double  Q;                      // node inflow (cfs)

static double  V;                      // node volume (ft3)

static double* R;                      // array of pollut. removals

static double* Cin;                    // node inflow concentrations

//static TTreatment* Treatment; // defined locally in treatmnt_treat()         //(5.1.008)

//-----------------------------------------------------------------------------

//  External functions (declared in funcs.h)

//-----------------------------------------------------------------------------

//  treatmnt_open           (called from routing_open)

//  treatment_close         (called from routing_close)

//  treatmnt_readExpression (called from parseLine in input.c)

//  treatmnt_delete         (called from deleteObjects in project.c)

//  treatmnt_setInflow      (called from qualrout_execute)

//  treatmnt_treat          (called from findNodeQual in qualrout.c)

//-----------------------------------------------------------------------------

//  Local functions

//-----------------------------------------------------------------------------

static int    createTreatment(int node);

static double getRemoval(int pollut);

static int    getVariableIndex(char* s);

static double getVariableValue(int varCode);

//=============================================================================

int  treatmnt_open(void)

//

//  Input:   none

//  Output:  returns TRUE if successful, FALSE if not

//  Purpose: allocates memory for computing pollutant removals by treatment.

//

{

    R = NULL;

    Cin = NULL;

    if ( Nobjects[POLLUT] > 0 )

    {

        R = (double *) calloc(Nobjects[POLLUT], sizeof(double));

        Cin = (double *) calloc(Nobjects[POLLUT], sizeof(double));

        if ( R == NULL || Cin == NULL)

        {

            report_writeErrorMsg(ERR_MEMORY, "");

            return FALSE;

        }

    }

    return TRUE;

}

//=============================================================================

void treatmnt_close(void)

//

//  Input:   none

//  Output:  returns an error code

//  Purpose: frees memory used for computing pollutant removals by treatment.

//

{

    FREE(R);

    FREE(Cin);

}

//=============================================================================

int  treatmnt_readExpression(char* tok[], int ntoks)

//

//  Input:   tok[] = array of string tokens

//           ntoks = number of tokens

//  Output:  returns an error code

//  Purpose: reads a treatment expression from a tokenized line of input.

//

{

    char  s[MAXLINE+1];

    char* expr;

    int   i, j, k, p;

    MathExpr* equation;                // ptr. to a math. expression

    // --- retrieve node & pollutant

    if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");

    j = project_findObject(NODE, tok[0]);

    if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);

    p = project_findObject(POLLUT, tok[1]);

    if ( p < 0 ) return error_setInpError(ERR_NAME, tok[1]);

    // --- concatenate remaining tokens into a single string

    strcpy(s, tok[2]);

    for ( i=3; i<ntoks; i++)

    {

        strcat(s, " ");

        strcat(s, tok[i]);

    }

    // --- check treatment type

    if      ( UCHAR(s[0]) == 'R' ) k = 0;

    else if ( UCHAR(s[0]) == 'C' ) k = 1;

    else return error_setInpError(ERR_KEYWORD, tok[2]);

    // --- start treatment expression after equals sign

    expr = strchr(s, '=');

    if ( expr == NULL ) return error_setInpError(ERR_KEYWORD, "");

    else expr++;

    // --- create treatment objects at node j if they don't already exist

    if ( Node[j].treatment == NULL )

    {

        if ( !createTreatment(j) ) return error_setInpError(ERR_MEMORY, "");

    }

    // --- create a parsed expression tree from the string expr

    //     (getVariableIndex is the function that converts a treatment

    //      variable's name into an index number)

    equation = mathexpr_create(expr, getVariableIndex);

    if ( equation == NULL )

        return error_setInpError(ERR_TREATMENT_EXPR, "");

    // --- save the treatment parameters in the node's treatment object

    Node[j].treatment[p].treatType = k;

    Node[j].treatment[p].equation = equation;

    return 0;

}

//=============================================================================

void treatmnt_delete(int j)

//

//  Input:   j = node index

//  Output:  none

//  Purpose: deletes the treatment objects for each pollutant at a node.

//

{

    int p;

    if ( Node[j].treatment )

    {

        for (p=0; p<Nobjects[POLLUT]; p++)

            mathexpr_delete(Node[j].treatment[p].equation);

        free(Node[j].treatment);

    }

    Node[j].treatment = NULL;

}

//=============================================================================

void  treatmnt_setInflow(double qIn, double wIn[])

//

//  Input:   j = node index

//           qIn = flow inflow rate (cfs)

//           wIn = pollutant mass inflow rate (mass/sec)

//  Output:  none

//  Purpose: computes and saves array of inflow concentrations to a node.

//

{

    int    p;

    if ( qIn > 0.0 )

        for (p = 0; p < Nobjects[POLLUT]; p++) Cin[p] = wIn[p]/qIn;

    else

        for (p = 0; p < Nobjects[POLLUT]; p++) Cin[p] = 0.0;

}

//=============================================================================

void  treatmnt_treat(int j, double q, double v, double tStep)

//

//  Input:   j     = node index

//           q     = inflow to node (cfs)

//           v     = volume of node (ft3)

//           tStep = routing time step (sec)

//  Output:  none

//  Purpose: updates pollutant concentrations at a node after treatment.

//

{

    int    p;                          // pollutant index

    double cOut;                       // concentration after treatment

    double massLost;                   // mass lost by treatment per time step

    TTreatment* treatment;             // pointer to treatment object          //(5.1.008)

    // --- set locally shared variables for node j

    if ( Node[j].treatment == NULL ) return;

    ErrCode = 0;

    J  = j;                            // current node

    Dt = tStep;                        // current time step

    Q  = q;                            // current inflow rate

    V  = v;                            // current node volume

    // --- initialze each removal to indicate no value

    for ( p = 0; p < Nobjects[POLLUT]; p++) R[p] = -1.0;

    // --- determine removal of each pollutant

    for ( p = 0; p < Nobjects[POLLUT]; p++)

    {

        // --- removal is zero if there is no treatment equation

        treatment = &Node[j].treatment[p];                                     //(5.1.008)

        if ( treatment->equation == NULL ) R[p] = 0.0;                         //(5.1.008)

        // --- no removal for removal-type expression when there is no inflow

              else if ( treatment->treatType == REMOVAL && q <= ZERO ) R[p] = 0.0;   //(5.1.008)

        // --- otherwise evaluate the treatment expression to find R[p]

        else getRemoval(p);

    }

    // --- check for error condition

    if ( ErrCode == ERR_CYCLIC_TREATMENT )

    {

         report_writeErrorMsg(ERR_CYCLIC_TREATMENT, Node[J].ID);

    }

    // --- update nodal concentrations and mass balances

    else for ( p = 0; p < Nobjects[POLLUT]; p++ )

    {

        if ( R[p] == 0.0 ) continue;

        treatment = &Node[j].treatment[p];                                     //(5.1.008)

       // --- removal-type treatment equations get applied to inflow stream

        if ( treatment->treatType == REMOVAL )                                 //(5.1.008)

        {

            // --- if no pollutant in inflow then cOut is current nodal concen.

            if ( Cin[p] == 0.0 ) cOut = Node[j].newQual[p];

            // ---  otherwise apply removal to influent concen.

            else cOut = (1.0 - R[p]) * Cin[p];

            // --- cOut can't be greater than mixture concen. at node

            //     (i.e., in case node is a storage unit)

            cOut = MIN(cOut, Node[j].newQual[p]);

        }

        // --- concentration-type equations get applied to nodal concentration

        else

        {

            cOut = (1.0 - R[p]) * Node[j].newQual[p];

        }

        // --- mass lost must account for any initial mass in storage

        massLost = (Cin[p]*q*tStep + Node[j].oldQual[p]*Node[j].oldVolume -

                   cOut*(q*tStep + Node[j].oldVolume)) / tStep;

        massLost = MAX(0.0, massLost);

        // --- add mass loss to mass balance totals and revise nodal concentration

        massbal_addReactedMass(p, massLost);

        Node[j].newQual[p] = cOut;

    }

}

//=============================================================================

int  createTreatment(int j)

//

//  Input:   j = node index

//  Output:  returns TRUE if successful, FALSE if not

//  Purpose: creates a treatment object for each pollutant at a node.

//

{

    int p;

    Node[j].treatment = (TTreatment *) calloc(Nobjects[POLLUT],

                                              sizeof(TTreatment));

    if ( Node[j].treatment == NULL )

    {

        return FALSE;

    }

    for (p = 0; p < Nobjects[POLLUT]; p++)

    {

        Node[j].treatment[p].equation = NULL;

    }

    return TRUE;

}

//=============================================================================

int  getVariableIndex(char* s)

//

//  Input:   s = name of a process variable or pollutant

//  Output:  returns index of process variable or pollutant

//  Purpose: finds position of process variable/pollutant in list of names.

//

{

    // --- check for a process variable first

    int k;

    int m = PVMAX;                     // PVMAX is number of process variables

    k = findmatch(s, ProcessVarWords);

    if ( k >= 0 ) return k;

    // --- then check for a pollutant concentration

    k = project_findObject(POLLUT, s);

    if ( k >= 0 ) return (k + m);

    // --- finally check for a pollutant removal

    if ( UCHAR(s[0]) == 'R' && s[1] == '_')

    {

        k = project_findObject(POLLUT, s+2);

        if ( k >= 0 ) return (Nobjects[POLLUT] + k + m);

    }

    return -1;

}

//=============================================================================

double getVariableValue(int varCode)

//

//  Input:   varCode = code number of process variable or pollutant

//  Output:  returns current value of variable

//  Purpose: finds current value of a process variable or pollutant concen.,

//           making reference to the node being evaluated which is stored in

//           shared variable J.

//

{

    int    p;

    double a1, a2, y;

    TTreatment* treatment;                                                     //(5.1.008)

    // --- variable is a process variable

    if ( varCode < PVMAX )

    {

        switch ( varCode )

        {

          case pvHRT:                                 // HRT in hours

            if ( Node[J].type == STORAGE )

            {

                return Storage[Node[J].subIndex].hrt / 3600.0;

            }

            else return 0.0;

          case pvDT:

            return Dt;                                // time step in seconds

          case pvFLOW:

            return Q * UCF(FLOW);                     // flow in user's units

          case pvDEPTH:

            y = (Node[J].oldDepth + Node[J].newDepth) / 2.0;

            return y * UCF(LENGTH);                   // depth in ft or m

          case pvAREA:

            a1 = node_getSurfArea(J, Node[J].oldDepth);

            a2 = node_getSurfArea(J, Node[J].newDepth);

            return (a1 + a2) / 2.0 * UCF(LENGTH) * UCF(LENGTH);

           

          default: return 0.0;

        }

    }

    // --- variable is a pollutant concentration

    else if ( varCode < PVMAX + Nobjects[POLLUT] )

    {

        p = varCode - PVMAX;

        treatment = &Node[J].treatment[p];                                     //(5.1.008)

        if ( treatment->treatType == REMOVAL ) return Cin[p];                  //(5.1.008)

        return Node[J].newQual[p];

    }

    // --- variable is a pollutant removal

    else

    {

        p = varCode - PVMAX - Nobjects[POLLUT];

        if ( p >= Nobjects[POLLUT] ) return 0.0;

        return getRemoval(p);

    }

}

//=============================================================================

double  getRemoval(int p)

//

//  Input:   p = pollutant index

//  Output:  returns fractional removal of pollutant

//  Purpose: computes removal of a specific pollutant

//

{

    double c0 = Node[J].newQual[p];    // initial node concentration

    double r;                          // removal value

    TTreatment* treatment;                                                     //(5.1.008)

    // --- case where removal already being computed for another pollutant

    if ( R[p] > 1.0 || ErrCode )

    {

        ErrCode = 1;

        return 0.0;

    }

    // --- case where removal already computed

    if ( R[p] >= 0.0 && R[p] <= 1.0 ) return R[p];

    // --- set R[p] to value > 1 to show that value is being sought

    //     (prevents infinite recursive calls in case two removals

    //     depend on each other)

    R[p] = 10.0;

    // --- case where current concen. is zero

    if ( c0 == 0.0 )

    {

        R[p] = 0.0;

        return 0.0;

    }

    // --- apply treatment eqn.

    treatment = &Node[J].treatment[p];                                         //(5.1.008)

    r = mathexpr_eval(treatment->equation, getVariableValue);                  //(5.1.008)

    r = MAX(0.0, r);

    // --- case where treatment eqn. is for removal

    if ( treatment->treatType == REMOVAL )                                     //(5.1.008)

    {

        r = MIN(1.0, r);

        R[p] = r;

    }

    // --- case where treatment eqn. is for effluent concen.

    else

    {

        r = MIN(c0, r);

        R[p] = 1.0 - r/c0;

    }

    return R[p];

}

//=============================================================================