JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases Wiki Activity

Merging changes from upstream dev

Browse Source
This commit is contained in:
Michael Tryby
2018-08-22 15:09:33 -04:00
parent fafad2a070 761199b08e
commit 75e400cf44
21 changed files with 1978 additions and 470 deletions
Download Patch File Download Diff File Expand all files Collapse all files

183
.github/BuildAndTest.md vendored Normal file
View File

@@ -0,0 +1,183 @@
## Building OWA EPANET From Source on Windows
by Michael E. Tryby
Created on: March 13, 2018
### Introduction
Building OWA's fork of EPANET from source is a basic skill that all developers
interested in contributing to the project should know how to perform. This
document describes the build process step-by-step. You will learn 1) how to
configure your machine to build the project locally; 2) how to obtain the
project files using git; 3) how to use cmake to generate build files and build
the project; and 4) how to use ctest and nrtest to perform unit and regression
testing on the build artifacts produced. Be advised, you will need local admin
privileges on your machine to follow this tutorial. Let’s begin!
### Dependencies
Before the project can be built the required tools must be installed. The OWA
EPANET project adheres to a platform compiler policy - for each platform there
is a designated compiler. The platform compiler for Windows is Visual
Studio cl, for Linux gcc, and for Mac OS clang. These instructions describe how
to build EPANET on Windows. CMake is a cross platform build, testing, and packaging
tool that is used to automate the EPANET build workflow. Boost is a free portable
peer-reviewed C++ library. Unit tests are linked with Boost unit test libraries.
Lastly, git is a free and open source distributed version control system. Git must
be installed to obtain the project source code from the OWA EPANET repository
found on GitHub.
### Summary of Build Dependencies
- Platform Compiler
- Windows: Visual Studio 10.0 32-bit cl (version 16.00.40219.01 for 80x86)
- CMake (version 3.0.0 or greater)
- Boost Libraries (version 1.58 or greater)
- git (version 2.6.0 or greater)
### Build Procedure
1. Install Dependencies
* A. Install Visual Studio 2010 Express and SP1
Our current benchmark platform and compiler is Windows 32-bit Visual Studio 10
2010. Older versions of Visual Studio are available for download here:
https://www.visualstudio.com/vs/older-downloads/
A service pack for Visual Studio 10 2010 is available here:
https://www.microsoft.com/en-us/download/details.aspx?id=34677
* B. Install Boost
Boost binaries for Windows offer a convenient installation solution. Be sure to
select for download the boost installer exe that corresponds to the version of Visual Studio you have installed.
https://sourceforge.net/projects/boost/files/boost-binaries/1.58.0/
Although newer version of Boost are available, a link to Boost 1.58 is provided. This is the library version that the unit tests have been written against. Older versions of Boost may not work. The default install location for the Boost
Libraries is C:\local\boost_1_58_0
* C. Install Chocolatey, CMake, and git
Chocolatey is a Windows package manager that makes installing some of these
dependencies a little easier. When working with Chocolatey it is useful to have
local admin privileges. Chocolatey is available here:
https://chocolatey.org/install
Once Chocolately is installed, from a command prompt running with admin privileges
issue the following commands
```
\>choco install -y cmake --installargs 'ADD_CMAKE_TO_PATH=User'
\>choco install -y git --installargs /GitOnlyOnPath
\>refreshenv
```
* D. Common Problems
Using chocolatey requires a command prompt with admin privileges.
Check to make sure installed applications are on the command path.
Make note of the Boost Library install location.
2. Build The Project
As administrator open a Visual Studio 2010 Command Prompt. Change directories to
the location where you wish to build the EPANET project. Now we will issue a series
of commands to create a parent directory for the project root and clone the project
from OWA's GitHub repository.
* A. Clone the EPANET Repository
```
\>mkdir OWA
\>cd OWA
\>git clone --branch=dev https://github.com/OpenWaterAnalytics/EPANET.git
\>cd EPANET
```
The present working directory is now the project root EPANET. The directory contains
the same files that are visibly present in the GitHub Repo by browsing to the URL
https://github.com/OpenWaterAnalytics/EPANET/tree/dev.
Now we will create a build products directory and generate the platform build
file using cmake.
* B. Generate the build files
```
\>mkdir buildprod
\>cd buildprod
\>set BOOST_ROOT=C:\local\boost_1_58_0
\>cmake -G "Visual Studio 10 2010" -DBOOST_ROOT="%BOOST_ROOT%" -DBoost_USE_STATIC_LIBS="ON" ..
```
Now that the dependencies have been installed and the build system has been
generated, building EPANET is a simple CMake command.
* C. Build EPANET
\>cmake --build . --config Debug
* D. Common Problems
CMake may not be able to find the project CMakeLists.txt file or the Boost
library install location.
3. Testing
Unit Testing uses Boost Unit Test library and CMake ctest as the test runner.
Cmake has been configured to register tests with ctest as part of the build process.
* A. Unit Testing
```
\>cd tests
\>ctest -C Debug
```
The unit tests run quietly. Ctest redirects stdout to a log file which can be
found in the "tests\Testing\Temporary" folder. This is useful when a test fails.
Regression testing is somewhat more complicated because it relies on Python
to execute EPANET for each test and compare the binary files and report files.
To run regression tests first python and any required packages must be installed.
If Python is already installed on your local machine the installation of
miniconda can be skipped.
* B. Installing Regression Testing Dependencies
```
cd ..\..
\>choco install -y miniconda --installargs '/AddToPath=1'
\>choco install -y curl
\>choco install -y 7zip
\>refreshenv
\>pip install -r tools/requirements-appveyor.txt
```
With Python and the necessary dependencies installed, regression testing can be run
using the before-test and run-nrtest helper scripts found in the tools folder. The script
before-test stages the test and benchmark files for nrtest. The script run-nrtest calls
nrtest execute and nrtest compare to perform the regression test.
To run the executable under test, nrtest needs the absolute path to it and a
unique identifier for it such as the version number. The project cmake build places build
artifacts in the buildprod\bin\ folder. On Windows the build configuration "Debug" or
"Release" must also be indicated. On Windows it is also necessary to specify the path to
the Python Scripts folder so the nrtest execute and compare commands can be found. You
need to substitute bracketed fields below like "<build identifier>" with the values for
your setup.
* C. Regression Testing
```
\>tools\before-test.cmd <relative path to regression test location> <absolute path to exe under test> <build identifier>
\>tools\run-nrtest.cmd <absolute path to python scripts> <relative path to regression test location> <build identifier>
```
* D. Common Problems
The batch file before-test.cmd needs to run with admin privileges. The nrtest script complains when it can't find manifest files.
That concludes this tutorial on building OWA EPANET from source on Windows.
You have learned how to configure your machine satisfying project dependencies
and how to acquire, build, and test EPANET on your local machine. To be sure,
there is a lot more to learn, but this is a good start! Learn more about project
build and testing dependencies by following the links provided below.
### Further Reading
* Visual Studio - https://msdn.microsoft.com/en-us/library/dd831853(v=vs.100).aspx
* CMake - https://cmake.org/documentation/
* Boost - http://www.boost.org/doc/
* git - https://git-scm.com/doc
* Miniconda - https://conda.io/docs/user-guide/index.html
* curl - https://curl.haxx.se/
* 7zip - https://www.7-zip.org/
* nrtest - https://nrtest.readthedocs.io/en/latest/

10
appveyor.yml
View File

@@ -17,13 +17,13 @@ matrix:
environment:
matrix:
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2013
GENERATOR: "Visual Studio 10 2010"
GROUP: "SUPPORTED"
BOOST_ROOT: "C:/Libraries/boost"
# New build on Visual Studio 15 2017
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
GENERATOR: "Visual Studio 15 2017"
GROUP: "SUPPORTED"
BOOST_ROOT: "C:/Libraries/boost_1_67_0"
# New build on Visual Studio 15 2017
- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
GENERATOR: "Visual Studio 15 2017 Win64"
GROUP: "EXPERIMENTAL"
BOOST_ROOT: "C:/Libraries/boost_1_67_0"

11
include/epanet2.bas
View File

@@ -56,6 +56,9 @@ Public Const EN_EFFICIENCY = 16
Public Const EN_HEADCURVE = 17
Public Const EN_EFFICIENCYCURVE = 18
Public Const EN_PRICEPATTERN = 19
Public Const EN_STATE = 20
Public Const EN_CONST_POWER = 21
Public Const EN_SPEED = 22
Public Const EN_DURATION = 0 ' Time parameters
Public Const EN_HYDSTEP = 1
@@ -154,6 +157,13 @@ Public Const EN_INITFLOW = 10 ' Re-initialize flow flag
Public Const EN_CONST_HP = 0 ' constant horsepower
Public Const EN_POWER_FUNC = 1 ' power function
Public Const EN_CUSTOM = 2 ' user-defined custom curve
Public Const EN_NOCURVE = 3 ' no curve
Public Const EN_V_CURVE = 0 ' volume curve
Public Const EN_P_CURVE = 1 ' pump curve
Public Const EN_E_CURVE = 2 ' efficiency curve
Public Const EN_H_CURVE = 3 ' head loss curve
Public Const EN_G_CURVE = 4 ' General\default curve
'These are the external functions that comprise the DLL
@@ -219,6 +229,7 @@ Public Const EN_CUSTOM = 2 ' user-defined custom curve
Declare Function ENgetcurve Lib "epanet2.dll" (ByVal curveIndex As Long, ByVal CurveID As String, nValues As Long, xValues As Any, yValues As Any) As Long
Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal pumpIndex As Long, curveIndex As Long) As Long
Declare Function ENgetpumptype Lib "epanet2.dll" (ByVal index As Long, PumpType As Long) As Long
Declare Function ENgetcurvetype Lib "epanet2.dll" (ByVal curveindex As Long, CurveType As Long) As Long
Declare Function ENgetversion Lib "epanet2.dll" (value As Long) As Long

79
include/epanet2.h
View File

@@ -115,7 +115,10 @@ typedef enum {
EN_EFFICIENCY = 16,
EN_HEADCURVE = 17,
EN_EFFICIENCYCURVE = 18,
EN_PRICEPATTERN = 19
EN_PRICEPATTERN = 19,
EN_STATE = 20,
EN_CONST_POWER = 21,
EN_SPEED = 22
} EN_LinkProperty;
/// Time parameter codes
@@ -224,8 +227,6 @@ typedef enum {
EN_TIMEOFDAY = 3
} EN_ControlType;
typedef enum {
EN_AVERAGE = 1, /* Time statistic types. */
EN_MINIMUM = 2, /* See TstatType in TYPES.H */
@@ -233,8 +234,6 @@ typedef enum {
EN_RANGE = 4
} EN_StatisticType;
typedef enum {
EN_MIX1 = 0, /* Tank mixing models */
EN_MIX2 = 1,
@@ -242,8 +241,6 @@ typedef enum {
EN_LIFO = 3
} EN_MixingModel;
typedef enum {
EN_NOSAVE = 0,
EN_SAVE = 1,
@@ -251,16 +248,21 @@ typedef enum {
EN_SAVE_AND_INIT = 11
} EN_SaveOption;
typedef enum {
EN_CONST_HP = 0, /* constant horsepower */
EN_POWER_FUNC = 1, /* power function */
EN_CUSTOM = 2 /* user-defined custom curve */
EN_CUSTOM = 2, /* user-defined custom curve */
EN_NOCURVE = 3 /* no curve */
} EN_PumpType;
typedef enum {
EN_V_CURVE = 0, /* volume curve */
EN_P_CURVE = 1, /* pump curve */
EN_E_CURVE = 2, /* efficiency curve */
EN_H_CURVE = 3, /* head loss curve */
EN_G_CURVE = 4 /* General\default curve */
} EN_CurveType;
// --- Declare the EPANET toolkit functions
#if defined(__cplusplus)
extern "C" {
@@ -763,10 +765,19 @@ extern "C" {
@param linkIndex The index of the pump element
@param[out] outType The integer-typed pump curve type signifier (output parameter)
@return Error code
@see EN_CurveType
@see EN_PumpType
*/
int DLLEXPORT ENgetpumptype(int linkIndex, int *outType);
/**
@brief Get the type of a curve
@param curveIndex The index of the curve element
@param[out] outType The integer-typed curve curve type signifier (output parameter)
@return Error code
@see EN_CurveType
*/
int DLLEXPORT ENgetcurvetype(int curveIndex, int *outType);
/**
@brief Get the version number. This number is to be interpreted with implied decimals, i.e., "20100" == "2(.)01(.)00"
@param[out] version The version of EPANET
@@ -1184,26 +1195,28 @@ extern "C" {
int DLLEXPORT EN_getqualtype(EN_ProjectHandle ph, int *qualcode, int *tracenode);
int DLLEXPORT EN_geterror(int errcode, char *errmsg, int maxLen);
int DLLEXPORT EN_getstatistic(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE* value);
int DLLEXPORT EN_getnodeindex(EN_ProjectHandle ph, char *id, int *index);
int DLLEXPORT EN_getnodeid(EN_ProjectHandle ph, int index, char *id);
int DLLEXPORT EN_getnodetype(EN_ProjectHandle ph, int index, int *code);
int DLLEXPORT EN_getnodevalue(EN_ProjectHandle ph, int index, int code, EN_API_FLOAT_TYPE *value);
int DLLEXPORT EN_getcoord(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y);
int DLLEXPORT EN_setcoord(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y);
int DLLEXPORT EN_getnumdemands(EN_ProjectHandle ph, int nodeIndex, int *numDemands);
int DLLEXPORT EN_getbasedemand(EN_ProjectHandle ph, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand);
int DLLEXPORT EN_getdemandpattern(EN_ProjectHandle ph, int nodeIndex, int demandIndex, int *pattIndex);
int DLLEXPORT EN_getlinkindex(EN_ProjectHandle ph, char *id, int *index);
int DLLEXPORT EN_getlinkid(EN_ProjectHandle ph, int index, char *id);
int DLLEXPORT EN_getlinktype(EN_ProjectHandle ph, int index, EN_LinkType *code);
int DLLEXPORT EN_setlinktype(EN_ProjectHandle ph, char *id, EN_LinkType type);
int DLLEXPORT EN_getlinknodes(EN_ProjectHandle ph, int index, int *node1, int *node2);
int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value);
int DLLEXPORT EN_getcurve(EN_ProjectHandle ph, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues);
int DLLEXPORT EN_getheadcurveindex(EN_ProjectHandle ph, int pumpIndex, int *curveIndex);
int DLLEXPORT EN_setheadcurveindex(EN_ProjectHandle ph, int pumpIndex, int curveIndex);
int DLLEXPORT EN_getpumptype(EN_ProjectHandle ph, int linkIndex, int *outType);
int DLLEXPORT EN_getstatistic(EN_Project *p, int code, EN_API_FLOAT_TYPE* value);
int DLLEXPORT EN_getnodeindex(EN_Project *p, char *id, int *index);
int DLLEXPORT EN_getnodeid(EN_Project *p, int index, char *id);
int DLLEXPORT EN_getnodetype(EN_Project *p, int index, int *code);
int DLLEXPORT EN_getnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE *value);
int DLLEXPORT EN_getcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y);
int DLLEXPORT EN_setcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y);
int DLLEXPORT EN_getnumdemands(EN_Project *p, int nodeIndex, int *numDemands);
int DLLEXPORT EN_getbasedemand(EN_Project *p, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand);
int DLLEXPORT EN_getdemandpattern(EN_Project *p, int nodeIndex, int demandIndex, int *pattIndex);
int DLLEXPORT EN_getlinkindex(EN_Project *p, char *id, int *index);
int DLLEXPORT EN_getlinkid(EN_Project *p, int index, char *id);
int DLLEXPORT EN_getlinktype(EN_Project *p, int index, EN_LinkType *code);
int DLLEXPORT EN_setlinktype(EN_Project *p, char *id, EN_LinkType type);
int DLLEXPORT EN_getlinknodes(EN_Project *p, int index, int *node1, int *node2);
int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value);
int DLLEXPORT EN_getcurve(EN_Project *p, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues);
int DLLEXPORT EN_getheadcurveindex(EN_Project *p, int pumpIndex, int *curveIndex);
int DLLEXPORT EN_setheadcurveindex(EN_Project *p, int pumpIndex, int curveIndex);
int DLLEXPORT EN_getpumptype(EN_Project *p, int linkIndex, int *outType);
int DLLEXPORT EN_getcurvetype(EN_Project *p, int curveIndex, int *outType);
int DLLEXPORT EN_getversion(int *version);
int DLLEXPORT EN_setcontrol(EN_ProjectHandle ph, int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level);
int DLLEXPORT EN_setnodevalue(EN_ProjectHandle ph, int index, int code, EN_API_FLOAT_TYPE v);

61
src/epanet.c
View File

@@ -436,6 +436,10 @@ int DLLEXPORT ENgetpumptype(int index, int *type) {
return EN_getpumptype(_defaultModel, index, type);
}
int DLLEXPORT ENgetcurvetype(int curveindex, int *type) {
return EN_getcurvetype(_defaultModel, curveindex, type);
}
int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands) {
return EN_getnumdemands(_defaultModel, nodeIndex, numDemands);
}
@@ -2231,6 +2235,7 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co
EN_API_FLOAT_TYPE *value) {
double a, h, q, v = 0.0;
int returnValue = 0;
int pmp;
EN_Project *p = (EN_Project*)ph;
@@ -2246,8 +2251,6 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co
double *LinkFlows = hyd->LinkFlows;
double *LinkSetting = hyd->LinkSetting;
/* Check for valid arguments */
*value = 0.0;
if (!p->Openflag)
@@ -2364,6 +2367,38 @@ int DLLEXPORT EN_getlinkvalue(EN_ProjectHandle ph, int index, EN_LinkProperty co
v = 1.0;
break;
case EN_STATE:
v = hyd->LinkStatus[index];
if (Link[index].Type == EN_PUMP) {
pmp = findpump(net, index);
if (hyd->LinkStatus[index] >= OPEN) {
if (hyd->LinkFlows[index] > hyd->LinkSetting[index] * Pump[pmp].Qmax)
v = XFLOW;
if (hyd->LinkFlows[index] < 0.0)
v = XHEAD;
}
}
break;
case EN_CONST_POWER:
v = 0;
if (Link[index].Type == EN_PUMP) {
pmp = findpump(net, index);
if (Pump[pmp].Ptype == CONST_HP) {
v = Link[index].Km; // Power in HP
}
}
break;
case EN_SPEED:
v = 0;
if (Link[index].Type == EN_PUMP) {
pmp = findpump(net, index);
v = Link[index].Kc;
}
break;
case EN_SETTING:
if (Link[index].Type == EN_PIPE || Link[index].Type == EN_CVPIPE) {
return (ENgetlinkvalue(index, EN_ROUGHNESS, value));
@@ -3187,7 +3222,7 @@ int DLLEXPORT EN_addcurve(EN_ProjectHandle ph, char *id) {
strcpy(tmpCur[n].ID, id);
tmpCur[n].Npts = 1;
tmpCur[n].Type = -1;
tmpCur[n].Type = G_CURVE;
tmpCur[n].X = (double *)calloc(tmpCur[n].Npts, sizeof(double));
tmpCur[n].Y = (double *)calloc(tmpCur[n].Npts, sizeof(double));
if (tmpCur[n].X == NULL)
@@ -3231,8 +3266,6 @@ int DLLEXPORT EN_setcurve(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE *x,
EN_Network *net = &p->network;
Scurve *Curve = net->Curve;
int j;
/* Check for valid arguments */
@@ -3265,12 +3298,9 @@ int DLLEXPORT EN_setcurvevalue(EN_ProjectHandle ph, int index, int pnt, EN_API_F
EN_Project *p = (EN_Project*)ph;
EN_Network *net = &p->network;
Scurve *Curve = net->Curve;
const int Ncurves = net->Ncurves;
if (!p->Openflag)
return (102);
if (index <= 0 || index > Ncurves)
@@ -3595,6 +3625,7 @@ int DLLEXPORT EN_setheadcurveindex(EN_ProjectHandle ph, int index, int curveinde
pump->Qmax /= Ucf[FLOW];
pump->Hmax /= Ucf[HEAD];
p->network.Curve[curveindex].Type = P_CURVE;
return (0);
}
@@ -3617,6 +3648,18 @@ int DLLEXPORT EN_getpumptype(EN_ProjectHandle ph, int index, int *type) {
return (0);
}
int DLLEXPORT EN_getcurvetype(EN_Project *p, int curveindex, int *type) {
EN_Network *net = &p->network;
if (!p->Openflag)
return (102);
if (curveindex < 1 || curveindex > net->Ncurves)
return (206);
*type = net->Curve[curveindex].Type;
return (0);
}
/*
----------------------------------------------------------------
Functions for opening files
@@ -3986,7 +4029,7 @@ int allocdata(EN_Project *p)
}
for (n = 0; n <= par->MaxCurves; n++) {
net->Curve[n].Npts = 0;
net->Curve[n].Type = -1;
net->Curve[n].Type = G_CURVE;
net->Curve[n].X = NULL;
net->Curve[n].Y = NULL;
}

18
src/funcs.h
View File

@@ -196,24 +196,8 @@ void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs.
/* ----------- SMATRIX.C ---------------*/
int createsparse(EN_Project *pr); /* Creates sparse matrix */
int allocsparse(EN_Project *pr); /* Allocates matrix memory */
void freesparse(EN_Project *pr); /* Frees matrix memory */
int buildlists(EN_Project *pr, int); /* Builds adjacency lists */
int paralink(EN_Project *pr, int, int, int); /* Checks for parallel links */
void xparalinks(EN_Project *pr); /* Removes parallel links */
void freelists(EN_Project *pr); /* Frees adjacency lists */
void countdegree(EN_Project *pr); /* Counts links at each node */
int reordernodes(EN_Project *pr); /* Finds a node re-ordering */
int mindegree(solver_t *s, int, int); /* Finds min. degree node */
int growlist(EN_Project *pr, int); /* Augments adjacency list */
int newlink(EN_Project *pr, Padjlist); /* Adds fill-ins for a node */
int linked(EN_Network *net, int, int); /* Checks if 2 nodes linked */
int addlink(EN_Network *net, int, int, int); /* Creates new fill-in */
int storesparse(EN_Project *pr, int); /* Stores sparse matrix */
int ordersparse(hydraulics_t *h, int); /* Orders matrix storage */
void transpose(int,int *,int *, /* Transposes sparse matrix */
int *,int *,int *,int *,int *);
int linsolve(solver_t *s, int); /* Solves set of inear eqns. */
int linsolve(EN_Project *pr, int); /* Solves set of linear eqns. */
/* ----------- QUALITY.C ---------------*/
int openqual(EN_Project *pr); /* Opens WQ solver system */

1001
src/genmmd.c Normal file
View File

File diff suppressed because it is too large Load Diff

2
src/hydsolver.c
View File

@@ -122,7 +122,7 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
hlosscoeff(pr, i);
}
matrixcoeffs(pr);
errcode = linsolve(&hyd->solver, net->Njuncs);
errcode = linsolve(pr, net->Njuncs);
/* Take action depending on error code */
if (errcode < 0) {

520
src/smatrix.c
View File

@@ -15,17 +15,21 @@ module are:
freesparse() -- called from closehyd() in HYDRAUL.C
linsolve() -- called from netsolve() in HYDRAUL.C
Createsparse() does the following:
createsparse() does the following:
1. for each node, builds an adjacency list that identifies
all links connected to the node (see buildlists())
2. re-orders the network's nodes to minimize the number
of non-zero entries in the hydraulic solution matrix
(see reorder())
3. converts the adjacency lists into a compact scheme
3. symbolically factorizes the solution matrix
(see factorize())
4. converts the adjacency lists into a compact scheme
for storing the non-zero coeffs. in the lower diagonal
portion of the solution matrix (see storesparse())
Freesparse() frees the memory used for the sparse matrix.
Linsolve() solves the linearized system of hydraulic equations.
freesparse() frees the memory used for the sparse matrix.
linsolve() solves the linearized system of hydraulic equations.
********************************************************************
*/
@@ -38,13 +42,54 @@ Linsolve() solves the linearized system of hydraulic equations.
#include <stdlib.h>
#endif
#include <math.h>
#include "hash.h"
#include <limits.h>
#include <time.h>
#include "text.h"
#include "types.h"
#include "epanet2.h"
#include "funcs.h"
#define EXTERN extern
#include "vars.h"
// The multiple minimum degree re-ordering routine (see genmmd.c)
extern int genmmd(int *neqns, int *xadj, int *adjncy, int *invp, int *perm,
int *delta, int *dhead, int *qsize, int *llist, int *marker,
int *maxint, int *nofsub);
// Local functions
static int allocsparse(EN_Project *pr);
static int buildlists(EN_Project *pr, int);
static int paralink(EN_Project *pr, int, int, int);
static void xparalinks(EN_Project *pr);
static void freelists(EN_Project *pr);
static void countdegree(EN_Project *pr);
static int reordernodes(EN_Project *pr);
static int factorize(EN_Project *pr);
static int growlist(EN_Project *pr, int);
static int newlink(EN_Project *pr, Padjlist);
static int linked(EN_Network *net, int, int);
static int addlink(EN_Network *net, int, int, int);
static int storesparse(EN_Project *pr, int);
static int sortsparse(EN_Project *pr, int);
static void transpose(int, int *, int *, int *, int *,
int *, int *, int *);
/*************************************************************************
* Timer macros
**************************************************************************/
//#define cleartimer(tmr) (tmr = 0.0)
//#define starttimer(tmr) (tmr -= ((double) clock()/CLOCKS_PER_SEC));
//#define stoptimer(tmr) (tmr += ((double) clock()/CLOCKS_PER_SEC));
//#define gettimer(tmr) (tmr)
/*************************************************************************
* The following data type implements a timer
**************************************************************************/
// typedef double timer;
// timer SmatrixTimer;
int createsparse(EN_Project *pr)
/*
@@ -56,11 +101,14 @@ int createsparse(EN_Project *pr)
*/
{
int errcode = 0;
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *solver = &pr->hydraulics.solver;
// cleartimer(SmatrixTimer);
// starttimer(SmatrixTimer);
/* Allocate data structures */
ERRCODE(allocsparse(pr));
@@ -70,37 +118,41 @@ int createsparse(EN_Project *pr)
}
/* Build node-link adjacency lists with parallel links removed. */
s->Degree = (int *) calloc(n->Nnodes+1, sizeof(int));
ERRCODE(MEMCHECK(s->Degree));
ERRCODE(buildlists(pr,TRUE));
if (!errcode){
xparalinks(pr); /* Remove parallel links */
countdegree(pr); /* Find degree of each junction */
} /* (= # of adjacent links) */
solver->Degree = (int *) calloc(net->Nnodes+1, sizeof(int));
ERRCODE(MEMCHECK(solver->Degree));
ERRCODE(buildlists(pr, TRUE));
if (!errcode)
{
xparalinks(pr); // Remove parallel links
countdegree(pr); // Find degree of each junction
} // (= # of adjacent links)
/* Re-order nodes to minimize number of non-zero coeffs. */
/* in factorized solution matrix. At same time, adjacency */
/* list is updated with links representing non-zero coeffs. */
hyd->Ncoeffs = n->Nlinks;
// Re-order nodes to minimize number of non-zero coeffs.
// in factorized solution matrix.
hyd->Ncoeffs = net->Nlinks;
ERRCODE(reordernodes(pr));
/* Allocate memory for sparse storage of positions of non-zero */
/* coeffs. and store these positions in vector NZSUB. */
ERRCODE(storesparse(pr,net->Njuncs));
// Factorize solution matrix by updating adjacency lists
// with non-zero connections due to fill-ins.
ERRCODE(factorize(pr));
/* Free memory used for adjacency lists and sort */
/* row indexes in NZSUB to optimize linsolve(). */
// Allocate memory for sparse storage of positions of non-zero
// coeffs. and store these positions in vector NZSUB.
ERRCODE(storesparse(pr, net->Njuncs));
// Free memory used for adjacency lists and sort
// row indexes in NZSUB to optimize linsolve().
if (!errcode) {
freelists(pr);
}
ERRCODE(ordersparse(hyd,net->Njuncs));
ERRCODE(sortsparse(pr, net->Njuncs));
/* Re-build adjacency lists without removing parallel */
/* links for use in future connectivity checking. */
// Re-build adjacency lists without removing parallel
// links for use in future connectivity checking.
ERRCODE(buildlists(pr,FALSE));
/* Free allocated memory */
free(s->Degree);
// Free allocated memory
free(solver->Degree);
return(errcode);
} /* End of createsparse */
@@ -114,18 +166,18 @@ int allocsparse(EN_Project *pr)
**--------------------------------------------------------------
*/
{
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
int errcode = 0;
n->Adjlist = (Padjlist *) calloc(n->Nnodes+1, sizeof(Padjlist));
s->Order = (int *) calloc(n->Nnodes+1, sizeof(int));
s->Row = (int *) calloc(n->Nnodes+1, sizeof(int));
s->Ndx = (int *) calloc(n->Nlinks+1, sizeof(int));
ERRCODE(MEMCHECK(n->Adjlist));
ERRCODE(MEMCHECK(s->Order));
ERRCODE(MEMCHECK(s->Row));
ERRCODE(MEMCHECK(s->Ndx));
net->Adjlist = (Padjlist *) calloc(net->Nnodes+1, sizeof(Padjlist));
solver->Order = (int *) calloc(net->Nnodes+1, sizeof(int));
solver->Row = (int *) calloc(net->Nnodes+1, sizeof(int));
solver->Ndx = (int *) calloc(net->Nlinks+1, sizeof(int));
ERRCODE(MEMCHECK(net->Adjlist));
ERRCODE(MEMCHECK(solver->Order));
ERRCODE(MEMCHECK(solver->Row));
ERRCODE(MEMCHECK(solver->Ndx));
return(errcode);
}
@@ -139,17 +191,22 @@ void freesparse(EN_Project *pr)
**----------------------------------------------------------------
*/
{
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
// stoptimer(SmatrixTimer);
// printf("\n");
// printf("\n Processing Time = %7.3f s", gettimer(SmatrixTimer));
// printf("\n");
freelists(pr);
free(n->Adjlist);
free(s->Order);
free(s->Row);
free(s->Ndx);
free(s->XLNZ);
free(s->NZSUB);
free(s->LNZ);
FREE(net->Adjlist);
FREE(solver->Order);
FREE(solver->Row);
FREE(solver->Ndx);
FREE(solver->XLNZ);
FREE(solver->NZSUB);
FREE(solver->LNZ);
} /* End of freesparse */
@@ -167,34 +224,34 @@ int buildlists(EN_Project *pr, int paraflag)
int errcode = 0;
Padjlist alink;
EN_Network *n = &pr->network;
EN_Network *net = &pr->network;
/* For each link, update adjacency lists of its end nodes */
for (k=1; k <= n->Nlinks; k++)
// For each link, update adjacency lists of its end nodes
for (k=1; k <= net->Nlinks; k++)
{
i = n->Link[k].N1;
j = n->Link[k].N2;
i = net->Link[k].N1;
j = net->Link[k].N2;
if (paraflag) {
pmark = paralink(pr,i,j,k); /* Parallel link check */
pmark = paralink(pr, i, j, k); // Parallel link check
}
/* Include link in start node i's list */
// Include link in start node i's list
alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist));
if (alink == NULL) return(101);
if (!pmark) alink->node = j;
else alink->node = 0; /* Parallel link marker */
else alink->node = 0; // Parallel link marker
alink->link = k;
alink->next = n->Adjlist[i];
n->Adjlist[i] = alink;
alink->next = net->Adjlist[i];
net->Adjlist[i] = alink;
/* Include link in end node j's list */
// Include link in end node j's list
alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist));
if (alink == NULL) return(101);
if (!pmark) alink->node = i;
else alink->node = 0; /* Parallel link marker */
else alink->node = 0; // Parallel link marker
alink->link = k;
alink->next = n->Adjlist[j];
n->Adjlist[j] = alink;
alink->next = net->Adjlist[j];
net->Adjlist[j] = alink;
}
return(errcode);
} /* End of buildlists */
@@ -215,13 +272,16 @@ int paralink(EN_Project *pr, int i, int j, int k)
Padjlist alink;
for (alink = pr->network.Adjlist[i]; alink != NULL; alink = alink->next)
{
if (alink->node == j) /* Link || to k (same end nodes) */
// Link || to k (same end nodes)
if (alink->node == j)
{
pr->hydraulics.solver.Ndx[k] = alink->link; /* Assign Ndx entry to this link */
// Assign Ndx entry to this link
pr->hydraulics.solver.Ndx[k] = alink->link;
return(1);
}
}
pr->hydraulics.solver.Ndx[k] = k; /* Ndx entry if link not parallel */
// Ndx entry if link not parallel
pr->hydraulics.solver.Ndx[k] = k;
return(0);
} /* End of paralink */
@@ -236,35 +296,35 @@ void xparalinks(EN_Project *pr)
*/
{
int i;
Padjlist alink, /* Current item in adjacency list */
blink; /* Previous item in adjacency list */
EN_Network *n = &pr->network;
Padjlist alink, // Current item in adjacency list
blink; // Previous item in adjacency list
EN_Network *net = &pr->network;
/* Scan adjacency list of each node */
for (i=1; i <= n->Nnodes; i++)
// Scan adjacency list of each node
for (i=1; i <= net->Nnodes; i++)
{
alink = n->Adjlist[i]; /* First item in list */
alink = net->Adjlist[i]; // First item in list
blink = NULL;
while (alink != NULL)
{
if (alink->node == 0) /* Parallel link marker found */
if (alink->node == 0) // Parallel link marker found
{
if (blink == NULL) /* This holds at start of list */
if (blink == NULL) // This holds at start of list
{
n->Adjlist[i] = alink->next;
free(alink); /* Remove item from list */
alink = n->Adjlist[i];
net->Adjlist[i] = alink->next;
free(alink); // Remove item from list
alink = net->Adjlist[i];
}
else /* This holds for interior of list */
else // This holds for interior of list
{
blink->next = alink->next;
free(alink); /* Remove item from list */
free(alink); // Remove item from list
alink = blink->next;
}
}
else
{
blink = alink; /* Move to next item in list */
blink = alink; // Move to next item in list
alink = alink->next;
}
}
@@ -283,14 +343,14 @@ void freelists(EN_Project *pr)
{
int i;
Padjlist alink;
EN_Network *n = &pr->network;
EN_Network *net = &pr->network;
for (i=0; i <= n->Nnodes; i++)
for (i=0; i <= net->Nnodes; i++)
{
for (alink = n->Adjlist[i]; alink != NULL; alink = n->Adjlist[i])
for (alink = net->Adjlist[i]; alink != NULL; alink = net->Adjlist[i])
{
n->Adjlist[i] = alink->next;
net->Adjlist[i] = alink->next;
free(alink);
}
}
@@ -308,14 +368,15 @@ void countdegree(EN_Project *pr)
{
int i;
Padjlist alink;
EN_Network *n = &pr->network;
memset(pr->hydraulics.solver.Degree,0,(n->Nnodes+1) * sizeof(int));
EN_Network *net = &pr->network;
/* NOTE: For purposes of node re-ordering, Tanks (nodes with */
/* indexes above Njuncs) have zero degree of adjacency. */
memset(pr->hydraulics.solver.Degree, 0, (net->Nnodes+1) * sizeof(int));
for (i=1; i <= n->Njuncs; i++) {
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) {
// NOTE: For purposes of node re-ordering, Tanks (nodes with
// indexes above Njuncs) have zero degree of adjacency.
for (i=1; i <= net->Njuncs; i++) {
for (alink = net->Adjlist[i]; alink != NULL; alink = alink->next) {
if (alink->node > 0) {
pr->hydraulics.solver.Degree[i]++;
}
@@ -334,59 +395,109 @@ int reordernodes(EN_Project *pr)
**--------------------------------------------------------------
*/
{
int k, knode, m, n;
EN_Network *net = &pr->network;
solver_t *s = &pr->hydraulics.solver;
int k, knode, m, njuncs, nlinks;
int delta = -1;
int nofsub = 0;
int maxint = INT_MAX; //defined in limits.h
int errcode;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
Padjlist alink;
// Local versions of node adjacency lists
int *adjncy = NULL;
int *xadj = NULL;
// Work arrays
int *dhead = NULL;
int *qsize = NULL;
int *llist = NULL;
int *marker = NULL;
// Default ordering
for (k=1; k <= net->Nnodes; k++)
{
s->Row[k] = k;
s->Order[k] = k;
solver->Row[k] = k;
solver->Order[k] = k;
}
n = net->Njuncs;
for (k=1; k<=n; k++) /* Examine each junction */
njuncs = net->Njuncs;
nlinks = net->Nlinks;
// Allocate memory
adjncy = (int *) calloc(2*nlinks+1, sizeof(int));
xadj = (int *) calloc(njuncs+2, sizeof(int));
dhead = (int *) calloc(njuncs+1, sizeof(int));
qsize = (int *) calloc(njuncs + 1, sizeof(int));
llist = (int *) calloc(njuncs + 1, sizeof(int));
marker = (int *) calloc(njuncs + 1, sizeof(int));
if (adjncy && xadj && dhead && qsize && llist && marker)
{
m = mindegree(s,k,n); /* Node with lowest degree */
knode = s->Order[m]; /* Node's index */
if (!growlist(pr,knode)) {
return(101); /* Augment adjacency list */
// Create local versions of node adjacency lists
xadj[1] = 1;
m = 1;
for (k = 1; k <= njuncs; k++)
{
for (alink = net->Adjlist[k]; alink != NULL; alink = alink->next)
{
knode = alink->node;
if (knode <= njuncs)
{
adjncy[m] = knode;
m++;
}
s->Order[m] = s->Order[k]; /* Switch order of nodes */
s->Order[k] = knode;
s->Degree[knode] = 0; /* In-activate node */
}
for (k=1; k<=n; k++) { /* Assign nodes to rows of */
s->Row[s->Order[k]] = k; /* coeff. matrix */
xadj[k+1] = m;
}
return(0);
// Generate a multiple minimum degree node re-ordering
genmmd(&njuncs, xadj, adjncy, solver->Row, solver->Order, &delta,
dhead, qsize, llist, marker, &maxint, &nofsub);
errcode = 0;
}
else errcode = 101; //insufficient memory
// Free memory
FREE(adjncy);
FREE(xadj);
FREE(dhead);
FREE(qsize);
FREE(llist);
FREE(marker);
return errcode;
} /* End of reordernodes */
int mindegree(solver_t *s, int k, int n)
int factorize(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: k = first node in list of active nodes
** n = total number of junction nodes
** Output: returns node index with fewest direct connections
** Purpose: finds active node with fewest direct connections
** Input: none
** Output: returns error code
** Purpose: symbolically factorizes the solution matrix in
** terms of its adjacency lists
**--------------------------------------------------------------
*/
{
int i, m;
int min = n,
imin = n;
int k, knode;
int errcode = 0;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
for (i=k; i<=n; i++)
// Augment each junction's adjacency list to account for
// new connections created when solution matrix is solved.
// NOTE: Only junctions (indexes <= Njuncs) appear in solution matrix.
for (k = 1; k <= net->Njuncs; k++) // Examine each junction
{
m = s->Degree[s->Order[i]];
if (m < min)
knode = solver->Order[k]; // Re-ordered index
if (!growlist(pr, knode)) // Augment adjacency list
{
min = m;
imin = i;
errcode = 101;
break;
}
solver->Degree[knode] = 0; // In-activate node
}
return(imin);
} /* End of mindegree */
return(errcode);
} /* End of factorize */
int growlist(EN_Project *pr, int knode)
@@ -402,17 +513,17 @@ int growlist(EN_Project *pr, int knode)
{
int node;
Padjlist alink;
EN_Network *n = &pr->network;
solver_t *s = &pr->hydraulics.solver;
EN_Network *net = &pr->network;
solver_t *solver = &pr->hydraulics.solver;
/* Iterate through all nodes connected to knode */
for (alink = n->Adjlist[knode]; alink != NULL; alink = alink -> next)
// Iterate through all nodes connected to knode
for (alink = net->Adjlist[knode]; alink != NULL; alink = alink -> next)
{
node = alink->node; /* End node of connecting link */
if (s->Degree[node] > 0) /* End node is active */
node = alink->node; // End node of connecting link
if (solver->Degree[node] > 0) // End node is active
{
s->Degree[node]--; /* Reduce degree of adjacency */
if (!newlink(pr,alink)) { /* Add to adjacency list */
solver->Degree[node]--; // Reduce degree of adjacency
if (!newlink(pr, alink)) { // Add to adjacency list
return(0);
}
}
@@ -433,35 +544,32 @@ int newlink(EN_Project *pr, Padjlist alink)
{
int inode, jnode;
Padjlist blink;
EN_Network *n = &pr->network;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *s = &pr->hydraulics.solver;
solver_t *solver = &pr->hydraulics.solver;
/* Scan all entries in adjacency list that follow anode. */
inode = alink->node; /* End node of connection to anode */
// Scan all entries in adjacency list that follow anode.
inode = alink->node; // End node of connection to anode
for (blink = alink->next; blink != NULL; blink = blink->next)
{
jnode = blink->node; /* End node of next connection */
jnode = blink->node; // End node of next connection
/* If jnode still active, and inode not connected to jnode, */
/* then add a new connection between inode and jnode. */
if (s->Degree[jnode] > 0) /* jnode still active */
// If jnode still active, and inode not connected to jnode,
// then add a new connection between inode and jnode.
if (solver->Degree[jnode] > 0) // jnode still active
{
if (!linked(n, inode,jnode)) { /* inode not linked to jnode */
/* Since new connection represents a non-zero coeff. */
/* in the solution matrix, update the coeff. count. */
if (!linked(net, inode, jnode)) // inode not linked to jnode
{
// Since new connection represents a non-zero coeff.
// in the solution matrix, update the coeff. count.
hyd->Ncoeffs++;
/* Update adjacency lists for inode & jnode to */
/* reflect the new connection. */
if (!addlink(n,inode,jnode,hyd->Ncoeffs)) {
return(0);
}
if (!addlink(n,jnode,inode,hyd->Ncoeffs)) {
return(0);
}
s->Degree[inode]++;
s->Degree[jnode]++;
// Update adjacency lists for inode & jnode to
// reflect the new connection.
if (!addlink(net, inode, jnode, hyd->Ncoeffs)) return(0);
if (!addlink(net, jnode, inode, hyd->Ncoeffs)) return(0);
solver->Degree[inode]++;
solver->Degree[jnode]++;
}
}
}
@@ -480,10 +588,9 @@ int linked(EN_Network *n, int i, int j)
*/
{
Padjlist alink;
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next) {
if (alink->node == j) {
return(1);
}
for (alink = n->Adjlist[i]; alink != NULL; alink = alink->next)
{
if (alink->node == j) return(1);
}
return(0);
} /* End of linked */
@@ -527,43 +634,43 @@ int storesparse(EN_Project *pr, int n)
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *s = &pr->hydraulics.solver;
solver_t *solver = &pr->hydraulics.solver;
/* Allocate sparse matrix storage */
s->XLNZ = (int *) calloc(n+2, sizeof(int));
s->NZSUB = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
s->LNZ = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
ERRCODE(MEMCHECK(s->XLNZ));
ERRCODE(MEMCHECK(s->NZSUB));
ERRCODE(MEMCHECK(s->LNZ));
if (errcode) {
return(errcode);
}
solver->XLNZ = (int *) calloc(n+2, sizeof(int));
solver->NZSUB = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
solver->LNZ = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
ERRCODE(MEMCHECK(solver->XLNZ));
ERRCODE(MEMCHECK(solver->NZSUB));
ERRCODE(MEMCHECK(solver->LNZ));
if (errcode) return(errcode);
/* Generate row index pointers for each column of matrix */
// Generate row index pointers for each column of matrix
k = 0;
s->XLNZ[1] = 1;
for (i=1; i<=n; i++) { /* column */
solver->XLNZ[1] = 1;
for (i=1; i<=n; i++) // column
{
m = 0;
ii = s->Order[i];
ii = solver->Order[i];
for (alink = net->Adjlist[ii]; alink != NULL; alink = alink->next)
{
j = s->Row[alink->node]; /* row */
j = solver->Row[alink->node]; // row
l = alink->link;
if (j > i && j <= n) {
if (j > i && j <= n)
{
m++;
k++;
s->NZSUB[k] = j;
s->LNZ[k] = l;
solver->NZSUB[k] = j;
solver->LNZ[k] = l;
}
}
s->XLNZ[i+1] = s->XLNZ[i] + m;
solver->XLNZ[i+1] = solver->XLNZ[i] + m;
}
return(errcode);
} /* End of storesparse */
int ordersparse(hydraulics_t *h, int n)
int sortsparse(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: n = number of rows in solution matrix
@@ -575,40 +682,45 @@ int ordersparse(hydraulics_t *h, int n)
int i, k;
int *xlnzt, *nzsubt, *lnzt, *nzt;
int errcode = 0;
solver_t *s = &h->solver;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *solver = &pr->hydraulics.solver;
int *LNZ = solver->LNZ;
int *XLNZ = solver->XLNZ;
int *NZSUB = solver->NZSUB;
xlnzt = (int *) calloc(n+2, sizeof(int));
nzsubt = (int *) calloc(h->Ncoeffs+2, sizeof(int));
lnzt = (int *) calloc(h->Ncoeffs+2, sizeof(int));
nzsubt = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
lnzt = (int *) calloc(hyd->Ncoeffs+2, sizeof(int));
nzt = (int *) calloc(n+2, sizeof(int));
ERRCODE(MEMCHECK(xlnzt));
ERRCODE(MEMCHECK(nzsubt));
ERRCODE(MEMCHECK(lnzt));
ERRCODE(MEMCHECK(nzt));
if (!errcode) {
/* Count # non-zeros in each row */
for (i=1; i<=n; i++) {
nzt[i] = 0;
}
for (i=1; i<=n; i++) {
for (k = s->XLNZ[i]; k < s->XLNZ[i+1]; k++) nzt[s->NZSUB[k]]++;
if (!errcode)
{
// Count # non-zeros in each row
for (i=1; i<=n; i++) nzt[i] = 0;
for (i=1; i<=n; i++)
{
for (k = XLNZ[i]; k < XLNZ[i+1]; k++) nzt[NZSUB[k]]++;
}
xlnzt[1] = 1;
for (i=1; i<=n; i++) xlnzt[i+1] = xlnzt[i] + nzt[i];
/* Transpose matrix twice to order column indexes */
transpose(n,s->XLNZ,s->NZSUB,s->LNZ,xlnzt,nzsubt,lnzt,nzt);
transpose(n,xlnzt,nzsubt,lnzt,s->XLNZ,s->NZSUB,s->LNZ,nzt);
// Transpose matrix twice to order column indexes
transpose(n, XLNZ, NZSUB, LNZ, xlnzt, nzsubt, lnzt, nzt);
transpose(n, xlnzt, nzsubt, lnzt, XLNZ, NZSUB, LNZ, nzt);
}
/* Reclaim memory */
free(xlnzt);
free(nzsubt);
free(lnzt);
free(nzt);
// Reclaim memory
FREE(xlnzt);
FREE(nzsubt);
FREE(lnzt);
FREE(nzt);
return(errcode);
} /* End of ordersparse */
} /* End of sortsparse */
void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
@@ -640,7 +752,7 @@ void transpose(int n, int *il, int *jl, int *xl, int *ilt, int *jlt,
} /* End of transpose */
int linsolve(solver_t *s, int n)
int linsolve(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: s = solver struct
@@ -666,13 +778,13 @@ int linsolve(solver_t *s, int n)
**--------------------------------------------------------------
*/
{
double *Aii = s->Aii;
double *Aij = s->Aij;
double *B = s->F;
int *LNZ = s->LNZ;
int *XLNZ = s->XLNZ;
int *NZSUB = s->NZSUB;
solver_t *solver = &pr->hydraulics.solver;
double *Aii = solver->Aii;
double *Aij = solver->Aij;
double *B = solver->F;
int *LNZ = solver->LNZ;
int *XLNZ = solver->XLNZ;
int *NZSUB = solver->NZSUB;
int *link, *first;
int i, istop, istrt, isub, j, k, kfirst, newk;

7
src/types.h
View File

@@ -102,7 +102,7 @@ typedef int INT4;
---------------------------------------------------------------------
*/
#define MEMCHECK(x) (((x) == NULL) ? 101 : 0 )
#define FREE(x) (free((x)))
#define FREE(x) if ((x)) free((x))
/*
---------------------------------------------------------------------
@@ -153,8 +153,9 @@ typedef enum {
V_CURVE, /* volume curve */
P_CURVE, /* pump curve */
E_CURVE, /* efficiency curve */
H_CURVE
} CurveType; /* head loss curve */
H_CURVE, /* head loss curve */
G_CURVE /* General\default curve */
} CurveType;
typedef enum {
CONST_HP, /* constant horsepower */

4
tools/before-test.cmd
View File

@@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0
set TEST_HOME=%~1
set EXAMPLES_VER=1.0.2-dev
set BENCHMARK_VER=220dev-vs17
set EXAMPLES_VER=1.0.2-dev.1
set BENCHMARK_VER=220dev1
set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip

4
tools/before-test.sh
View File

@@ -22,8 +22,8 @@
SCRIPT_HOME="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
TEST_HOME=$1
EXAMPLES_VER="1.0.1"
BENCHMARK_VER="2012vs10"
EXAMPLES_VER="1.0.2-dev.1"
BENCHMARK_VER="220dev1"
TEST_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v${EXAMPLES_VER}.tar.gz"
BENCH_URL="https://github.com/OpenWaterAnalytics/epanet-example-networks/releases/download/v${EXAMPLES_VER}/epanet-benchmark-${BENCHMARK_VER}.tar.gz"

43
tools/nrtest-epanet/main.py
View File

@@ -4,7 +4,9 @@ import time
import cStringIO
import itertools as it
import epanet_reader as er
# project import
import nrtest_epanet.output_reader as er
def result_compare(path_test, path_ref, comp_args):
@@ -15,28 +17,35 @@ def result_compare(path_test, path_ref, comp_args):
total = 0
output = cStringIO.StringIO()
eps = np.finfo(float).eps
min_cdd = 100.0
start = time.time()
test_reader = er.reader(path_test)
ref_reader = er.reader(path_ref)
test_reader = er.output_generator(path_test)
ref_reader = er.output_generator(path_ref)
for test, ref in it.izip(test_reader, ref_reader):
total += 1
if total%100000 == 0:
print(total)
if test.size != ref.size:
if len(test[0]) != len(ref[0]):
raise ValueError('Inconsistent lengths')
# Skip results if they are zero or equal
if np.array_equal(test, ref):
equal += 1
continue
#if np.array_equal(test, ref):
# equal += 1
# continue
else:
try:
np.testing.assert_allclose(test, ref, 1.0e-06, 2*eps)
close += 1
diff = np.fabs(np.subtract(test[0], ref[0]))
idx = np.unravel_index(np.argmax(diff), diff.shape)
if diff[idx] != 0.0:
tmp = - np.log10(diff[idx])
if tmp < min_cdd:
min_cdd = tmp;
except AssertionError as ae:
notclose += 1
@@ -49,8 +58,9 @@ def result_compare(path_test, path_ref, comp_args):
print(output.getvalue())
output.close()
print('equal: %d close: %d notclose: %d total: %d in %f (sec)\n' %
(equal, close, notclose, total, (stop - start)))
print('mincdd: %d in %f (sec)' % (np.floor(min_cdd), (stop - start)))
#print('equal: %d close: %d notclose: %d total: %d in %f (sec)\n' %
# (equal, close, notclose, total, (stop - start)))
if notclose > 0:
print('%d differences found\n' % notclose)
@@ -120,6 +130,8 @@ def nrtest_execute(app_path, test_path, output_path):
exit(not success)
import report_diff as rd
if __name__ == "__main__":
# app_path = "apps\\swmm-5.1.11.json"
# test_path = "tests\\examples\\example1.json"
@@ -130,7 +142,10 @@ if __name__ == "__main__":
# ref_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012"
# print(nrtest_compare(test_path, ref_path, (0.001, 0.0)))
benchmark_path = "C:\\Users\\mtryby\\Workspace\\GitRepo\\michaeltryby\\epanet-lr\\nrtestsuite\\benchmarks\\"
path_test = benchmark_path + "epanet-220dev\\example2\\example2.out"
path_ref = benchmark_path + "epanet-2012\\example2\\example2.out"
#result_compare(path_test, path_ref, (0.001, 0.0))
rd.report_diff(path_test, path_ref, 2)
path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2011a\\Example_3\\example3.out"
path_ref = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\epanet-testsuite\\benchmarks\\v2012\\Example_3\\example3.out"
result_compare(path_test, path_ref, (0.001, 0.0))

65
tools/nrtest-epanet/nrtest_epanet/__init__.py
View File

@@ -28,7 +28,7 @@ __copyright__ = "None"
__credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell"
__license__ = "CC0 1.0 Universal"
__version__ = "0.4.0"
__version__ = "0.5.0"
__date__ = "September 6, 2017"
__maintainer__ = "Michael Tryby"
@@ -38,7 +38,7 @@ __status = "Development"
def epanet_allclose_compare(path_test, path_ref, rtol, atol):
'''
Compares results in two EPANET binary files. Using the comparison criteria
Compares results in two EPANET binary files using the comparison criteria
described in the numpy assert_allclose documentation.
(test_value - ref_value) <= atol + rtol * abs(ref_value)
@@ -67,22 +67,67 @@ def epanet_allclose_compare(path_test, path_ref, rtol, atol):
for (test, ref) in it.izip(ordr.output_generator(path_test),
ordr.output_generator(path_ref)):
if len(test) != len(ref):
if len(test[0]) != len(ref[0]):
raise ValueError('Inconsistent lengths')
# Skip over arrays that are equal
if np.array_equal(test, ref):
if np.array_equal(test[0], ref[0]):
continue
else:
np.testing.assert_allclose(test, ref, rtol, atol)
np.testing.assert_allclose(test[0], ref[0], rtol, atol)
return True
# def epanet_better_compare(path_test, path_ref, rtol, atol):
# '''
# If you don't like assert_allclose you can add another function here.
# '''
# pass
def epanet_mincdd_compare(path_test, path_ref, rtol, atol):
'''
Compares the results of two EPANET binary files using a correct decimal
digits (cdd) comparison criteria:
min cdd(test, ref) >= atol
Returns true if min cdd in the file is greater than or equal to atol,
otherwise an AssertionError is thrown.
Arguments:
path_test - path to result file being testedgit
path_ref - path to reference result file
rtol - ignored
atol - minimum allowable cdd value (i.e. 3)
Returns:
True
Raises:
ValueError()
AssertionError()
'''
min_cdd = 100.0
for (test, ref) in it.izip(ordr.output_generator(path_test),
ordr.output_generator(path_ref)):
if len(test[0]) != len(ref[0]):
raise ValueError('Inconsistent lengths')
# Skip over arrays that are equal
if np.array_equal(test[0], ref[0]):
continue
else:
diff = np.fabs(np.subtract(test[0], ref[0]))
idx = np.unravel_index(np.argmax(diff), diff.shape)
if diff[idx] != 0.0:
tmp = - np.log10(diff[idx])
if tmp < min_cdd:
min_cdd = tmp;
if np.floor(min_cdd) >= atol:
return True
else:
raise AssertionError('min_cdd=%d less than atol=%g' % (min_cdd, atol))
def epanet_report_compare(path_test, path_ref, rtol, atol):
'''

6
tools/nrtest-epanet/nrtest_epanet/output_reader.py
View File

@@ -23,7 +23,8 @@ def output_generator(path_ref):
yield element attributes. It is useful for comparing contents of binary
files for numerical regression testing.
The generator yields a Python list containing element attributes.
The generator yields a Python tuple containing an array of element
attributes and a tuple containing the element type, period, and attribute.
Arguments:
path_ref - path to result file
@@ -38,7 +39,8 @@ def output_generator(path_ref):
for element_type in oapi.ElementType:
for attribute in br.elementAttributes[element_type]:
yield br.element_attribute(element_type, period_index, attribute)
yield (br.element_attribute(element_type, period_index, attribute),
(element_type, period_index, attribute))
class OutputReader():

96
tools/nrtest-epanet/report_diff.py Normal file
View File

@@ -0,0 +1,96 @@
# -*- coding: utf-8 -*-
#
# report_diff.py
#
# Date Created: July 11, 2018
#
# Author: Michael E. Tryby
# US EPA - ORD/NRMRL
#
# system imports
import itertools as it
# third party imports
import numpy as np
# project imports
import nrtest_epanet.output_reader as ordr
def _binary_diff(path_test, path_ref, min_cdd):
for (test, ref) in it.izip(ordr.output_generator(path_test),
ordr.output_generator(path_ref)):
if len(test[0]) != len(ref[0]):
raise ValueError('Inconsistent lengths')
# Skip over arrays that are equal
if np.array_equal(test[0], ref[0]):
continue
else:
lre = _log_relative_error(test[0], ref[0])
idx = np.unravel_index(np.argmin(lre), lre.shape)
if lre[idx] < min_cdd:
_print_diff(idx, lre, test, ref)
return
def _log_relative_error(q, c):
'''
Computes log relative error, a measure of numerical accuracy.
Single precision machine epsilon is between 2^-24 and 2^-23.
Reference:
McCullough, B. D. "Assessing the Reliability of Statistical Software: Part I."
The American Statistician, vol. 52, no. 4, 1998, pp. 358-366.
'''
diff = np.subtract(q, c)
tmp_c = np.copy(c)
# If ref value is small compute absolute error
tmp_c[np.fabs(tmp_c) < 1.0e-6] = 1.0
re = np.fabs(diff)/np.fabs(tmp_c)
# If re is tiny set lre to number of digits
re[re < 1.0e-7] = 1.0e-7
# If re is very large set lre to zero
re[re > 2.0] = 1.0
return np.negative(np.log10(re))
def _print_diff(idx, lre, test, ref):
idx_val = (idx[0], ref[1])
test_val = (test[0][idx[0]])
ref_val = (ref[0][idx[0]])
diff_val = (test_val - ref_val)
lre_val = (lre[idx[0]])
print("Idx: %s\nSut: %f Ref: %f Diff: %f LRE: %f\n"
% (idx_val, test_val, ref_val, diff_val, lre_val))
def report(args):
_binary_diff(args.test, args.ref, args.mincdd)
if __name__ == '__main__':
from argparse import ArgumentParser
parser = ArgumentParser(description='EPANET benchmark difference reporting')
parser.set_defaults(func=report)
parser.add_argument('-t', '--test', default=None,
help='Path to test benchmark')
parser.add_argument('-r', '--ref', default=None,
help='Path to reference benchmark')
parser.add_argument('-mc', '--mincdd', type=int, default=3,
help='Minimum correct decimal digits')
args = parser.parse_args()
args.func(args)

3
tools/nrtest-epanet/setup.py
View File

@@ -17,6 +17,7 @@ except ImportError:
entry_points = {
'nrtest.compare': [
'epanet allclose = nrtest_epanet:epanet_allclose_compare',
'epanet mincdd = nrtest_epanet:epanet_mincdd_compare',
'epanet report = nrtest_epanet:epanet_report_compare',
# Add entry point for new comparison functions here
]
@@ -24,7 +25,7 @@ entry_points = {
setup(
name='nrtest-epanet',
version='0.4.0',
version='0.5.0',
description="EPANET extension for nrtest",
author="Michael E. Tryby",

2
tools/run-nrtest.cmd
View File

@@ -18,7 +18,7 @@ setlocal
set NRTEST_SCRIPT_PATH=%~1
set TEST_SUITE_PATH=%~2
set BENCHMARK_VER=220dev-vs17
set BENCHMARK_VER=220dev1
set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute

2
tools/run-nrtest.sh
View File

@@ -19,7 +19,7 @@ run-nrtest()
return_value=0
test_suite_path=$1
benchmark_ver="2012vs10"
benchmark_ver="220dev1"
nrtest_execute_cmd="nrtest execute"

8
win_build/WinSDK/Makefile.bat
View File

@@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp
If %ERRORLEVEL% == 1 (
CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release
rem : create EPANET2.DLL
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
md "%Build_PATH%"\64bit
move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit
move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit
@@ -35,9 +35,9 @@ rem : 32 bit with DEF
CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release
echo "32 bit with epanet2.def mapping"
rem : create EPANET2.DLL
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
md "%Build_PATH%"\32bit
move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit
move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit

1
win_build/WinSDK/epanet2.def
View File

@@ -93,3 +93,4 @@ EXPORTS
ENdeletelink = _ENdeletelink@4
ENdeletenode = _ENdeletenode@4
ENsetlinktype = _ENsetlinktype@8
ENgetcurvetype = _ENgetcurvetype@8