Water Supply analysis program

[JohnNichols]

Dear mecej4:

I am teaching a class about water supply analysis, using the program Magni.  I have not played with the program since your very kind help in making the program.  I went back to the 3000 node model from an Australian town and have created the data files for Magni.  I am slowly debugging the input file, when it is working I will put up a copy.  I have about 350 of the nodes and pipes in and it appears to loop and work nicely.  Even when it fails to converge, it runs 1000 times in 24 seconds. 

I was racking my brain trying to remember where the heck the program came from -- I worked through this site reading the old posts when you got PARDISO working, but I was none the wiser. It has to come from somewhere - there is a C# version I think?  My students are doing a town of about 10000 people. I think I will need to look to add pump input.  

I wish I could remember.  I know you did a huge part of the program, I can see your fingerprints all over it. (Getting to old)

JMN

[JohnNichols]

https://ascelibrary.org/doi/10.1061/%28ASCE%29WR.1943-5452.0000596

People are trying to speed up the epanet version -- I wonder which is faster., c or fortran

[mecej4]

Here is a thread that you initiated in 2014 in the old MKL forum.

You had a C# program named WSA, which according to your description was derived from EPANET and used a linear equations solver from a C# library called MathNet. Based on recommendations from me and others in the MKL forum, you modified your C# program to use the MKL-Pardiso sparse solver, and then you were persuaded to redo the program in Fortran for speed. It was at this point that I intervened to get the Fortran version cleaned up and to improve the organization of the code, in particular the collecting of the contributions to the system matrix and passing correct arguments to Pardiso.

Your other question was regarding C versus Fortran. The two languages are equally effective, when used by experienced programmers. Fortran may yield 10 to 20 percent faster speed. C programs are usually much more difficult to debug because the language makes it easy to make a mess. Fortran has far better support for array operations. MKL supports Fortran, C and other languages such as Python. Considering that we are in 2020, a more appropriate comparison would be C++ vs Fortran.

I skimmed through the ASCE paper that you linked to (authors in an Austrian group). It is interesting, but the main conclusion struck me to be counter to what I remember from the 2014 exercise with introducing Pardiso as the solver -- that using Pardiso made the calculations significantly faster as compared to the solvers in MathNet or Lapack.

I see that a new version of EPANET, version 2.2, is now on Github , released in 2020. You probably built Magni based on the 2.0 (2008) version. Apart from the older version being in Fortran and the new version in C, I do not know what the modifications are.

Do you have EPANET input files for a problem that takes, say, between 10 s and 60 s to run? You wrote the file I/O parts of Magni. Is the input file format identical to that used in EPANET?

Let us also note that the Pardiso in MKL is an offshoot from Pardiso 4, and Schenk's Pardiso 5 and 6 (www.pardisoproject.org) are not quite the same as the Pardiso in MKL.

[JohnNichols]

Thank you for finding that stuff. I did not think to look on MKL, I do not really watch that forum.  

I remember the great problem getting the matrix into the correct form for PARDISO and I know you developed the final code for assembling the matrix in the correct form.  I use that code in a lot of programs now. Pardiso is just so fast and easy to use.   I cannot find the program I wrote in C# from the EPANET code, but I remember that I thought mathnet solver was a bit slow and that is when I tried to use PARDISO with C# -- but your suggestion was that FORTRAN would run twice as fast and I remember now a weekend coding from C# to Fortran and you were more than correct.  The structs however are straight from your ideas,  you showed me how to do them in some structures code from Harrison in 1973. 

I got onto accelerometers and put this to one side, and there was no real use for it.  Then this class came along and they are supposed to learn about management of facilities. I thought they might be interested in a water supply problem.  There are 100 of them so 10 nodes each and we have a 1000 node problem and I gave them a contour plan so they can get real heights. 

But I course I had to get MAGNI running on the latest IF compiler and VS. But that was no problem.  I do not like the EPANET text file, it is to complex and reading it in FORTRAN is a pain. So I gave them the simple text file format and they are going.  Yesterday, I thought what the heck, I will have a look at the old TAMWORTH file, it is in Autocad and I have LISP routines to create the data files.  2000+ nodes.

In the old days, AUTOLISP ran really slowly on Autocad, so it would take several hours to create a data file.  But now in runs on a core i5 in a minute.  So I have the 2000+ node data file.  The only issue is slowly adding nodes so you can debug the errors - adding reservoirs etc..  I like this program more than EPANET because I added a dxf file drawer and I can recreate the drawing from the data file so you have a visual check.  Only problem is you need Autocad.  

So I stopped typing this and went for a real hunt for the C# code, and finally found it. It runs really slowly compared to MAGNI, but it includes pumps so I can add pumps to MAGNI, which I need for TAMWORTH and the students.  

I include the complete node and pipe file for TAMWORTH, it just takes a while to add so many nodes and not get overwhelmed with errors.  

Here is the complete MAGNI at the moment, the T.inp is the start of the Tamworth file, I have added a bunch of fake reservoirs to make sure all pipes are connected to reservoirs, as I add pipes this will slowly get fixed.  

It seems to work a treat with the loops, it is fast.   The students will add the feature all good programs need to develop properly, a user to find the bugs and make suggestions.  

My only problem is the students do not have access to IF, so I have to give them the exe file, but never done that with PARDISO included.  HELP.

Anyway as always this was possible because of your help - thanks. 

My students are enjoying the experience. It helps to have real problems. 

Finally - I look at sone of my LISP programs and they are clean neat and easy to read, I look at PYTHON and I just see a mess.  If Fortran just had a few lisp features.  The interesting thing is the LISP commands you find in C.  Deep in K+R's minds must have been some dormant desire to take LISP and FORTRAN and create a beautiful language.  They failed. 

JMN

[JohnNichols]

Anyone can download a free Autocad viewer, this allows you to view the tam drawing.  If you tell it to show the overall view it works.   VIEW NAME = OVERALL

if you do an extents, you pick up the last node number at 0,0.  I have to store it somewhere, or else search all of the nodes each time I start.  Impossibly slow in the old days- but possible now. 

A system of using the XY coordinates to generate node numbers would be much better, tracking node numbers is a pain and causes most of the errors. The node numbers are in red and stored 0.25 m up and across in X for each node point. 

That was the interesting bit of LISP code to always recover the node numbers correctly.  I remember it took a while to get code that did not miss.  The thing to remember TAMwat dwg was created in 1988 and it still can be easily used today. 

Anyway MAGNI in Fortran is much faster than any other program for this work - I think -- the ASCE paper has got to be wrong MATHNET is not as fast as Pardiso, not by a long shot. 

[JohnNichols]

Three interesting errors. 

1. If you have two pipes going from the same nodes - the program stops execution and give me an access violation on the pardiso - need to check for this problem and signal error - do you think this is due to a singularity?

2. If I duplicate a pipe I get the same error, easy to do when you are slowly copying pipes over - need to check for this problem and signal error

3.  The pump curve works, but the checker to make sure you are connected to a reservoir gets fooled by pumps and if I have consecutive pipes a-b-c        g-h-i and later in the list c-d-e-f-g and then the access checker gets fooled. 

I had never really paid attention before, EPANET which was used to provide a simple C# algorithm did not use the node or the pipe numbers internally, the design from L Rossman generated an internal consecutive set, I copied this into the FORTRAN code and so the nodes and the pipe numbers in the file are just names.  Interesting, but the drawing package throws a fit if they are not ordered. 

This file has about 500 nodes and pipes.  it is just slow work going thru the system. 

I ran mathnet in C# as per that ASCE paper - somewhat - and it is slow as wet paper towel rolling down a hill in no wind. 

it is very fast - allowing that it is writing a dxf file as well that is slow. 

[mecej4]

JMN:

I could not build an EXE with the source files in Magni2019.zip because manicsr.f90 contains (between lines 110 and 120) a call to FILES(), and the sources do not include the code for that subroutine.

Please advise how to proceed.

Another questionable line of code: In misccsr.f90, we see the following statement

Qm = ((sqrt(temp3))*...

Since temp3 = (dh * dia * 2.0 * 9.806), and the head difference dh can be positive or negative, the calculation of Qm may cause the program to abort or take improper actions when dh is negative. 

[JohnNichols]

My apologies, I have several versions of the program, I have a small set of standard libraries I keep separate.  Here they are.  You need both files. 

If we are going to make the program fit for use, I need to go back and make sure it provides the same answers as the EPANET program and find the errors, such as you point out.  The problems come out as the size of the data files grows.   One slowly eliminates or traps the errors. 

I also need to bite the bullet and ensure the program reads the EPANET standard input files,  so we can run a number of standard problems.  

It is an interesting problem and working with the students, the Washington State Water Supply guy, the fire engineer and you make it fun. 

The Fire Engineer and the WatSupGuy are saying that the network has velocity limits, that I just cannot believe they achieve in the real network. 

Thanks for the help.  The three groups that sell this type of software must make a lot of money - it is an expensive program.  (I am not selling this - It is merely for people like students to learn with and to see how fast we can make it and try and set some real limits in the codes.  )

[JohnNichols]

dh does go negative , I put a stop code in - need to work out the work around.  

The interesting challenge is the flow equation has two main forms.  The Hazen Williams equation is an approximate equation, for reasons I have never fathomed the average US Civil engineer is really keen on the HW equation, the closer solution is the Colebrook White form of the Darcy Weisbach equation, they are harder to code, but much better, but in going from EPANET to MAGNI there is this additional change. 

I have a moral objection to using the HW - it is just a long way from reality. 

[jimdempseyatthecove]

John,

Does the head variable not only incorporate a static head value, but perceived head due to flow rate and/or momentum effects?

Jim Dempsey

[JohnNichols]

The static head are what we are after,  we assume for each iteration they are fixed at a node and we determine the loss between nodes as a pipe loss, we do not add in any of the other issues,  so a valve has a loss and a 90 degree band etc, but we ignore them. 

The negative head just means the water is flowing in the opposite direction 

[JohnNichols]

  
                        Qm = ((sqrt(temp3))*(dia*dia)*(3.1416/4))                 ! Check flow calcs.
                        if(dh .lt. 0.0) then 
                        
                         qm = 0.0
                         write(*,101)k2, qm
101                      Format( "Dh is negative :: ",i4,"   : ", f12.6)
                        endif 

QM is never actually used, I put in a trap to look at it, but the following equations catch negative heads  and calculate negative flows in the pipes.  

 spipes(k2)%fr = fr
                        spipes(k2)%q = sign(pi*dia*kvis*re*0.25D0, dh)
                        write(*,99)dh,spipes(k2)%q, Qm
                        write(*,'(F20.8, 5x, F10.8)')dh,spipes(k2)%q

The off-diagonal elements receive only two contribution per pipe,
! since it is assumed that the number of pipes connecting two nodes is
! either 0 or 1. The irn, jcn, jac arrays can be accumulated without
! regard to order.

This needs to be checked - as it easy to connect two nodes twice by mistake - 

[jimdempseyatthecove]

I haven't looked at your code, the code above is somewhat indicative that should the head be negative, that you intend to set the flow rate to 0. While this may be appropriate to do at an outlet, it may not be advisable when your pipes interconnect multiple water towers.

Jim Dempsey

[JohnNichols]

Qm was not used and has now been removed.  

The equation gives you a flow in the reverse direction as a negative flow and a negative head, but the equation is correct.  

I have started to add in a lot of output to a log file so that it can be checked and I am slowly working through thinking of all the errors to trap and trapping them.  

I went and had a look at the original EPANET code that the C# code was derived from - there is no error checking I could see, if your code gives you a problem that will not solve you get a ill conditioned problem.  

I also have a number of standard Hardy Cross problems that I have solved and they give the correct answers. 

The problem is at about 400 pipes, you start to get weird errors - usually a data entry error and trapping them is painful and tedious, better to write that code now on a small problem and then test my big problem. It is much more efficient.  

Just thinking of the potential errors in data is interesting exercise.  

My problem is I have no idea how to use the ifort CMD compiler to compile the code with MKL Pardiso.    I will need to find the manual and read it.  I am a VS person,  1988 last time I used MS compiler and even then I used VEDIT to run the compiler. 

But if I want the students to do a realistic problem we can get the student version of KY pipe which is node number limited or I can get Bentley that is a pain to load with a floating licence or I can use this one and enjoy talking to you all.  Option 3 much better and it is much faster.  

I have tried EPANET with some small problems like 10 nodes and it runs slower than the 400 node problems on PARDISO - observable slower.  

We are speed freaks otherwise we would code in Python.  

Plus telling the Fire Engineer and the Water supply guy they are wrong about their limits is fun.  

JMN

[JohnNichols]

One of the interesting challenges. Rossman who wrote the original EPA code in C - provided an ID number for each unit - pipe, pump etc, but in his code it is really a character name.  The problem in Fortran is a character name adds unneeded complexity and is fraught with errors.  So I used integers for the ID.  

If I count the pipes and store them in an array that is strictly monotonic ascending then one can allow the ID numbers to not sync with the pipe number, but that introduces a potential for mistaking the pipe number in the printout. 

I think I should enforce a strict one to one on the ID and the internal pipe number, but we lose some flexibility -- anyone got any thoughts - the real problem is the challenged character use in Fortran. 

[jimdempseyatthecove]

The ID/pipe number is not a property of the pipe (or other object)

Your diagram should be a collection of "things". A thing (you provide name) has multiple properties, one of which is the ID, another is a type (example a pipe), and a third is a portal connector array (array of size 2 for pipe, and n in the case of a manifold). A portal connector can be: not connected, capped, connected to thing+port.

Jim Dempsey

[JohnNichols]

Yes, I agree, but the problem is we do not live in a perfect world - how tightly do you suggest we limit the name of the things - EPANET calls it ID -- and it can be ABCD  or 1 or 34.  

[JohnNichols]

Fortran structs are interesting with this water supply model 

TYPE spipe
        INTEGER id  !/* Pipe ID number */
        INTEGER n1  !/* Start node index of pipe */
        INTEGER n2  !/* End node index of pipe */
        INTEGER :: zone = 0
        INTEGER pumpstatus
        INTEGER nodei, nodej
        
        INTEGER IDnode1, IDnode2
        REAL (KIND=dp) :: a0 = 0.0
        REAL (KIND=dp) a1
        REAL (KIND=dp) a2
        REAL (KIND=dp) a3
        REAL (KIND=dp) length
        REAL (KIND=dp) diameter !/* Diam. of pipe (in) */
        REAL (KIND=dp) c
        REAL (KIND=dp) wlrf !pipe wall roughness, m
        REAL (KIND=dp) e
        REAL (KIND=dp) pumpq
        REAL (KIND=dp) deltaq
        REAL (KIND=dp) q
        REAL (KIND=dp) h0
        REAL (KIND=dp) h1
        REAL (KIND=dp) h2
        REAL (KIND=dp) h3
        REAL (KIND=dp) h4
        REAL (KIND=dp) qt
        REAL (KIND=dp) v
        REAL (KIND=dp) :: k = 0.0
        REAL (KIND=dp) :: k1
        REAL (KIND=dp) :: RE
        REAL (KIND=dp) :: fr
        integer :: cv = 0
        integer :: prv = 0
        REAL (KIND=dp) :: prvcoeffic = 1.0
        REAL (KIND=dp) :: prvzeropressure
        REAL (KIND=dp) kvis !kinematic viscosity, m^2/s
        REAL (KIND=dp) tc !water temperature, celsius
    END TYPE

 - it turns out as soon as the model gets to beyond a size you can see the node numbers on a single plot - you start to lose node numbers, then you have to take the ID route and set up a computer model with its own numbering linked to the id's.  Son of a Gunn 

[jimdempseyatthecove]

Am I to assume that a node is what the pipe is attached to?

Where the node may be a coupler, a T, a hydrant, ...

Either the pipe needs to hold the "port" on the node as well as the node number, or the node needs to contain the pipeID (or ID of whatever is attache to the node, and in which case an identifier as to what that is).

The input (initialization) data should contain the interconnects, and the input procedure should detect errors. Example, multiple pipes attached to same port on node, and/or pipe attaching to non-existant node+port, and for post input, a test for open ports.

The back bone is connected to the hip bone, the hip bone is connected to the thigh bone,...

Jim Dempsey

[JohnNichols]

LINK : fatal error LNK1561: entry point must be defined

ifort -I"%MKLROOT%"\include -c misccsr.f90 draw.for magnicsr.f90 scotia.for karilib.for
link magnicsr.obj draw.obj scotia.obj karilib.obj misccsr.obj mkl_intel_c.lib mkl_intel_thread.lib mkl_core.lib libiomp5md.lib

--------------------------------------------------------------------------------------------------------------------

I am using a compiler 2020 - update 2 - IA32 cmd prompt. 

The batch file is above and I get the link entry point error, but magnicsr.f90 has a "program magni"

line and it gives the link error - if I drop off magni I get an error identifier missing - everything with Intel Fortran is stock standard.

---------------------------------------------------------------------------------------------------------------------

The only things we have are reservoirs and pipes at the moment.  A set of nodes are created -- each has a location and an ID, I use numbers for the ID.  Node 1 is a reservoir.  pipes span from node to node and some nodes have flow demands. 

As finally evolved by mecej4 a few years ago, the program works fine. The problems start when you get to large models, it is very easy to make mistakes as the model complexity grows and setting error traps to pick them all up is tedious.  The way forward is slow development of a large model and find all the errors.  There are some interesting ones.  

I put some of the standard errors in the epanet program and it simply tells you illconditioned matrix, which is ok with 100 nodes, but 1000 nodes it is tedious to find and eliminate.  

But is is getting better. 

[mecej4]

Do not use LINK directly as you reported, but use IFORT to do compiling and linking:

ifort /Qmkl misccsr.f90 draw.for magnicsr.f90 scotia.for karilib.for /Femagnicsr

If you have compiled the source files individually and all the OBJ files have been generated earlier, leaving only linking to be done, the command would be

ifort /Qmkl misccsr.obj draw.obj magnicsr.obj scotia.obj karilib.obj /Femagnicsr