zoukankan      html  css  js  c++  java
  • [水力建模]EPANET代码解读1

    EPANET代码解读1   

    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/

  • 相关阅读:
    LeetCode 258 Add Digits
    LeetCode 231 Power of Two
    LeetCode 28 Implement strStr()
    LeetCode 26 Remove Duplicates from Sorted Array
    LeetCode 21 Merge Two Sorted Lists
    LeetCode 20 Valid Parentheses
    图形处理函数库 ImageTTFBBox
    php一些函数
    func_get_arg(),func_get_args()和func_num_args()的用法
    人生不是故事,人生是世故,摸爬滚打才不会辜负功名尘土
  • 原文地址:https://www.cnblogs.com/KingOfFreedom/p/2805117.html
Copyright © 2011-2022 走看看