Sparse Matrix Format in smatrix.c

[NicoRenaud]

Dear EPANET devs, thanks a lot for all the work you are putting in developing and maintaining this code !

I'm trying to understand the sparse matrix format used in smatrix.c. When I try to reconstruct the matrix I obtain a non-symmetric matrix. However I think that the Cholesky factorization used to solve the linear system assumes that the matrix is symmetric so I'm missing something.

I am using a very simple network (see below). I'm using EPANET through the WNTR interface with the following code

import wntr

# Create a water network model
inp_file = "networks/Net0.inp"
wn = wntr.network.WaterNetworkModel(inp_file)

# Graph the network
wntr.graphics.plot_network(wn, title=wn.name, node_labels=True)

# run the simulation
sim = wntr.sim.EpanetSimulator(wn)  #
results = sim.run_sim()

At the first call of linsolve I get this matrix for the linear system :

A = [ 0.116  -2.339]
    [-0.116   2.454]

with this b vector:

b = [ -1.614  230.277]

If I compute the Cholesky decomposition of this matrix and use the forward/backward substitution just as in linsolve I do obtain the solution given by linsolve, i.e.:

x =  [83.8   97.772]

However this solution vector is not a solution of the linear system defined above:

A @ x = [-218.967,  230.211] != b

While the second component of the vector is nearly correct, the first one is not. Is it an expected behavior or am I doing something wrong ?

Sorry for this rather basic question about the code but I couldn't find the answer anywhere. Thanks a lot for your help !

---

[TITLE]

[JUNCTIONS] ;ID Elev Demand Pattern
J1 0 0 ; D1 0 50 ;

[RESERVOIRS] ;ID Head Pattern
R1 30 ;

[TANKS] ;ID Elevation InitLevel MinLevel MaxLevel Diameter MinVol VolCurve Overflow

[PIPES] ;ID Node1 Node2 Length Diameter Roughness MinorLoss Status P1 R1 J1 100 250 0.05 0 Open ; P2 J1 D1 1000 200 0.04 0 Open ;

[PUMPS] ;ID Node1 Node2 Parameters

[VALVES] ;ID Node1 Node2 Diameter Type Setting MinorLoss

[TAGS]

[DEMANDS] ;Junction Demand Pattern Category

[STATUS] ;ID Status/Setting

[PATTERNS] ;ID Multipliers

[CURVES] ;ID X-Value Y-Value

[CONTROLS]

[RULES]

[ENERGY] Global Efficiency 75 Global Price 0 Demand Charge 0

[EMITTERS] ;Junction Coefficient

[QUALITY] ;Node InitQual

[SOURCES] ;Node Type Quality Pattern

[REACTIONS] ;Type Pipe/Tank Coefficient

[REACTIONS] Order Bulk 1 Order Tank 1 Order Wall 1 Global Bulk 0 Global Wall 0 Limiting Potential 0 Roughness Correlation 0

[MIXING] ;Tank Model

[TIMES] Duration 0 Hydraulic Timestep 1:00 Quality Timestep 0:05 Pattern Timestep 1:00 Pattern Start 0:00 Report Timestep 1:00 Report Start 0:00 Start ClockTime 12 am Statistic None

[REPORT] Status No Summary No Page 0

[OPTIONS] Units LPS Headloss D-W Specific Gravity 1 Viscosity 1 Trials 40 Accuracy 0.001 CHECKFREQ 2 MAXCHECK 10 DAMPLIMIT 0 Unbalanced Continue 10 Pattern 1 Demand Multiplier 1.0 Emitter Exponent 0.5 Quality None mg/L Diffusivity 1 Tolerance 0.01

[COORDINATES] ;Node X-Coord Y-Coord J1 1000.000 6000.000
D1 11000.000 6000.000
R1 -1172.214 7424.023

[VERTICES] ;Link X-Coord Y-Coord

[LABELS] ;X-Coord Y-Coord Label & Anchor Node

[BACKDROP] DIMENSIONS 0.000 0.000 10000.000 10000.000
UNITS None FILE
OFFSET 0.00 0.00

[END]

[NicoRenaud]

I've explored a bit more and I may have answered my question. I'm guessing that the elements on the upper triangular of the matrix are irrelevant. The Cholesky decomposition computes the L matrix as if the sparse matrix was symmetric with the elements of the lower triangular transposed on the upper triangular. Am I correct ? Thank for your help

[LRossman]

I suggest you read Chapters 1, 2, 3 and 5 of the Computer Solution of Sparse Linear Systems book on which the EPANET sparse linear solver is based.

[NicoRenaud]

Thanks !