EPANET是美国环境保护局（环保署）公布的水力分析引擎，其源代码是开放的，但由于是用C写的，并且代码清晰度并不够，因此有必要解读一下。

EPANET当前最新版本是2.0，可以编译为动态链接库（DLL），也可以编译为独立的可执行程序（EXE）。编译DLL的最好方法是在VC++新建一个空的DLL项目，然后把*.c、*.h以及epanet.def文件加入到项目中，即可顺利完成编译。要编译EXE，可以新建一个空的EXE项目，同样加入文件后，定义CLE并进行编译。

DLL包含了一组导出函数，详细的用法在用户手册中有介绍，此处不作解释。EXE需要3个文件参数，分别是输入文件、报告文件和输出文件，下面从main函数开始根据正常执行流程进行说明。

main函数定义在epanet.c文件的起始位置。除了参数检查之外，最重要的是调用ENepanet函数（也是DLL中的导出函数，定义在epanet.c文件中）进行完整的模拟。ENepanet函数只包含几行代码，依次是调用ENopen、ENsolveH、ENsolveQ、ENreport和ENclose，这几个函数也都是从外部可调用的导出函数，也都定义在epanet.c文件中，分别解读如下。

ENopen需要的也是和main同样的3个文件参数。首先，函数初始化一些系统标志，包括Openflag OpenHflag OpenQflag SaveHflag SaveQflag Warnflag Messageflag Rptflag等，调用epanet.c中定义的initpointers函数初始化一系列指针（设置为NULL）。这些指针保存了管网的各种信息，其定义都在一个单独的文件var.h中。虽然指针的名字看起来有点乱，但因为非常重要，根据初始化的顺序详细列举如下：

D——各节点实际用水需求，为double*类型；

C——各节点实际水质，为double*类型；

H——各节点的水头（head ），为double*类型；

Q——各管段的流量，为double*类型；

R——各管段的反应率，为double*类型；

S——各管段的状态，为char*类型，值为CLOSED OPEN ACTIVE 等枚举状态（定义在types.h中）；

K——各管段的设置，如水泵速度、阀的设置等，为double*类型；

OldStat——各管段（水池）的原来状态；

Node——节点数据，为Snode*类型，Snode为一个结构，其定义在types.h中，包括一个字符串类型的ID、一个Sdemand*类型的需求指针、一个Ssource*类型的源指针、double类型的标高、初始水质、扩散器系数和一个char类型的报告标志；

Link——管段数据，为Slink*类型，Slink为一个结构，其定义在types.h中，包括一个字符串类型的ID、两个int类型的起止节点索引、double类型的直径、长度、粗糙度、最小损失系数、主流区反应速率系数Kb、管壁反应速率系数Kw、水流阻力以及char类型的管段类型、初始状态、报告标志；

Tank——水池数据，为Stank*类型，水池是一种特殊类型的节点，Stank除包含一个节点索引外，还包括double类型的面积、最小标高/水位、最大标高/水位、初始标高/水位、最小容量、最大容量、初始容量、反应系数、水池容量、浓度、混合室体积以及int类型的Fixed grade time pattern、容量-水位曲线索引和char类型的混合模型；

Pump——水泵数据，为Spump*类型，水泵是一种特殊类型的管段，Spump除包含一个管段索引外，还包括int类型的水泵曲线类型、水头-流量曲线索引、效率-流量曲线索引等，另外还有能量消耗、成本等属性（暂略）；

Valve——阀门数据，为Svalve*类型，阀门也是一种特殊类型的管段，没有单独的属性；

Pattern——时间模式，Curve——曲线，Control——控制，这几个是“非物理组件”，可以参见用户手册。

X——通用类型的double数组；Patlist——临时的时间模式列表；Curvelist——临时的曲线模型列表；

Adjlist——节点邻接列表指针，指向Sadjlist类型的节点邻接列表，Sadjlist类型包括节点索引、管段索引以及指向下一个节点邻接列表的指针。

Aii、Aij——矩阵方程A*H=F中矩阵A的对角线和非对角线元素，为double*类型。

F——矩阵方程A*H=F中的向量F，为double*类型。

P——管段水头损失关于流量求导的倒数，为double*类型，详细参见用户手册。

Y——流量校正因子，详细参见用户手册。

Order、Row、Ndx——与矩阵A相关的一些系数。

Nht、Lht——ID哈稀表。

EPANET代码解读2 

ENopen调用的initpointers除了初始化前述指针外，还调用定义在rules.c文件中的initrules函数将RuleState初始化为r_PRIORITY（枚举类型），将Rule变量初始化为NULL。Rule是一个指向aRule类型变量的指针，aRule类型包括字符串类型的标签、double类型的优先度、Premise类型的前提、Action类型的真假条件下的两个动作以及指向下一个规则的指针。

ENopen然后调用openfiles实际打开文件，然后调用netsize函数计算管网规模。netsize定义在input2.c文件中，其方法是逐行读取输入文件，并根据各段的有效行数计算节点、水管、水池、水泵、阀门、控制、规则、时间模式、曲线的最大数量。ENopen紧接着调用allocdata函数为管网数据分配内存。allocdata(位于rules.c)函数的功能包括创建Nht和Lht两个哈稀表(hash.c内)，分配节点、管段、水池、水泵、阀门、控制、时间模式、曲线的数据（设置前述指针），并初始化时间模式、曲线、节点实际用水需求指针，最后调用rules.c中的allocrules函数分配规则的内存。

ENopen在allocdata之后继续调用getdata函数从文件中读取输入数据。getdata定义在input1.c文件中，先调用setdefaults设置一些全局变量的初始值，再调用initreport初始化报告选项，然后调用rewind回转到输入文件起始位置，接着调用input2.c文件中定义的readdata函数读取输入文件的内容。readdata和netsize有点类似，也是逐行读入文件，并根据所在的节调用input2.c中的newline(int sect, char *line)函数处理一行。newline有一个长的switch-case分发列表，根据不同的节调用定义在input3.c中的不同的处理函数，包括juncdata() tankdata() pipedata() pumpdata() valvedata() patterndata() curvedata() demanddata() controldata() ruledata() sourcedata() emitterdata() qualdata() statusdata() energydata() reactdata() mixingdata() reportdata() timedata() optiondata()等。目前公开的epanet2.0代码并没有读入坐标等数据，对绘图有较大影响，需要添加相应的函数和变量。getdata在调用readdata之后，还调用了adjustdata initunits inittanks convertunits等函数进行一些处理，此处暂略。

ENopen在调用getdata之后，必要时会调用openhydfile打开水力模拟结果文件（不存在时创建、存在时检查格式但不读取结果内容），最后返回。

ENepanet函数中接着会根据是否有水力模拟结果文件决定是否调用ENsolveH。ENsolveH进行各时段的水力计算，包括调用ENopenH、ENinitH进行初始化，以及循环调用ENrunH、ENnextH进行各时段计算，最后调用ENcloseH结束水力计算。水力计算后续将深入分析，此处跳过。

ENepanet函数然后调用ENsolveQ进行各时段的水质计算。也包括了类似的初始化、循环计算等函数调用，后续再作分析。

最后ENepanet调用ENreport创建报告，调用ENclose释放内存和关闭文件。

EPANET代码解读3 

EPANET进行水力分析最重要的是ENrunH中调用的runhyd函数和ENnextH中调用的nexthyd函数。

（一）runhyd定义在hydraul.c文件中，依次调用了以下函数：

1、demands——计算节点当前时间各节点的用水需求。

2、controls——基于时间或水池水位控制各管段的设置，返回当前时间需要设置的管段数。

3、netsolve——水力分析的核心，求解水力方程A*H=F，采用的方法是Todini梯度方法，通过多次迭代调用linsolve(Njuncs,Aii,Aij,F)（定义在SMATRIX.C中）进行求解。netsolve在间隔一定次数的linsolve调用后，检查收敛条件。

int   runhyd(long *t)
/*
**--------------------------------------------------------------
**  Input:   none    
**  Output:  t = pointer to current time (in seconds)
**  Returns: error code                                         
**  Purpose: solves network hydraulics in a single time period
**--------------------------------------------------------------
*/
{
   int   iter;                          /* Iteration count   */
   int   errcode;                       /* Error code        */
   double relerr;                        /* Solution accuracy */

/* Find new demands & control actions */
   *t = Htime;
   demands();
   controls();

/* Solve network hydraulic equations */
   errcode = netsolve(&iter,&relerr);
   if (!errcode)
   {
      /* Report new status & save results */
      if (Statflag) writehydstat(iter,relerr);

/*** Updated 3/1/01 ***/
      /* If system unbalanced and no extra trials */
      /* allowed, then activate the Haltflag.     */
      if (relerr > Hacc && ExtraIter == -1) Haltflag = 1;

/* Report any warning conditions */
      if (!errcode) errcode = writehydwarn(iter,relerr);
   }
   return(errcode);
}                               /* end of runhyd */

（二）nexthyd定义在hydraul.c文件中，主要调用了以下函数：

1、timestep——计算到下一个水力模拟的时间（水池放空或注满、控制激活等）。

2、addenergy——增加水泵的耗能，以及计算水泵的一些统计信息。

int  nexthyd(long *tstep)
/*
**--------------------------------------------------------------
**  Input:   none    
**  Output:  tstep = pointer to time step (in seconds)
**  Returns: error code                                         
**  Purpose: finds length of next time step & updates tank
**           levels and rule-based contol actions
**--------------------------------------------------------------
*/
{
   long  hydstep;         /* Actual time step  */
   int   errcode = 0;     /* Error code        */

/*** Updated 3/1/01 ***/
   /* Save current results to hydraulics file and */
   /* force end of simulation if Haltflag is active */
   if (Saveflag) errcode = savehyd(&Htime);
   if (Haltflag) Htime = Dur;

/* Compute next time step & update tank levels */
   *tstep = 0;
   hydstep = 0;
   if (Htime < Dur) hydstep = timestep();
   if (Saveflag) errcode = savehydstep(&hydstep);

/* Compute pumping energy */
   if (Dur == 0) addenergy(0);
   else if (Htime < Dur) addenergy(hydstep);

/* Update current time. */
   if (Htime < Dur)  /* More time remains */
   {
      Htime += hydstep;
      if (Htime >= Rtime) Rtime += Rstep;
   }
   else
   {
      Htime++;          /* Force completion of analysis */
   }
   *tstep = hydstep;
   return(errcode);
}

本帖转自：http://blog.163.com/qisiliang274@126/blog/static/9394478620110391013980/