Tuesday, January 16, 2018

Robert Dickinson on Earn.Com

My referral link. You can sign up to Earn Bitcoin from emails sent to your email address

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];
}

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

Future and Backward Year Planning for an Interaction between SWMM and RED (me)

A forward and backward look at what i was doing in the future and past. The future is grossly overestimated for the number of years, but I like the view from the pivot year of 2018. It shows a range of SWMM and versions of SWMM involvement.
Forward Year
Forward Age
Backward Year
Backward Age
What was I doing?
2018
62
2018
62
Innovyze InfoSWMM/InfoSewer/ICM/XPSWMM and SWMM5/6
2019
63
2017
61
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2020
64
2016
60
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2021
65
2015
59
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2022
66
2014
58
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2023
67
2013
57
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2024
68
2012
56
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2025
69
2011
55
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2026
70
2010
54
Innovyze InfoSWMM/InfoSewer/ICM/SWMM5
2027
71
2009
53
MWHGlobal InfoSWMM/InfoSewer/SWMM5
2028
72
2008
52
MWHGlobal InfoSWMM/InfoSewer/SWMM5
2029
73
2007
51
CDM and SWMM5
2030
74
2006
50
CDM and SWMM5
2031
75
2005
49
CDM and SWMM5
2032
76
2004
48
CDM and SWMM5
2033
77
2003
47
CDM and SWMM4/5
2034
78
2002
46
CDM and SWMM4/5
2035
79
2001
45
CDM and SWMM4/5
2036
80
2000
44
Visual SWMM/SWMM4
2037
81
1999
43
Visual SWMM/SWMM4
2038
82
1998
42
XPSWMM/SWMM4
2039
83
1997
41
XPSWMM/SWMM4
2040
84
1996
40
XPSWMM/SWMM4
2041
85
1995
39
XPSWMM/SWMM4
2042
86
1994
38
XPSWMM/SWMM4
2043
87
1993
37
XPSWMM/SWMM4
2044
88
1992
36
XPSWMM/SWMM4
2045
89
1991
35
SWMM4
2046
90
1990
34
SWMM4
2047
91
1989
33
SWMM4
2048
92
1988
32
SWMM4
2049
93
1987
31
SWMM4
2050
94
1986
30
SWMM4
2051
95
1985
29
SWMM3
2052
96
1984
28
SWMM3
2053
97
1983
27
SWMM3
2054
98
1982
26
SWMM3
2055
99
1981
25
SWMM3
2056
100
1980
24
SWMM3
2057
101
1979
23
SWMM3
2058
102
1978
22
SWMM3
2059
103
1977
21
SWMM3

AI Rivers of Wisdom about ICM SWMM

Here's the text "Rivers of Wisdom" formatted with one sentence per line: [Verse 1] 🌊 Beneath the ancient oak, where shadows p...