/*
 ******************************************************************************
 Project:      OWA EPANET
 Version:      2.3
 Module:       leakage.c
 Description:  models additional nodal demands due to pipe leaks.
 Authors:      see AUTHORS
 Copyright:    see AUTHORS
 License:      see LICENSE
 Last Updated: 06/14/2024
 ******************************************************************************
*/
/*
This module uses the FAVAD (Fixed and Variable Discharge) equation to model
leaky pipes:
  
  Q = Co * L * (Ao + m * H) * sqrt(H)
  
where Q = pipe leak flow rate, Co = an orifice coefficient (= 0.6*sqrt(2g)),
L = pipe length, Ao = initial area of leak per unit of pipe length,
m = change in leak area per unit of pressure head, and H = pressure head.

The inverted form of this equation is used to model the leakage demand from
a pipe's end node using a pair of equivalent emitters as follows:

  H = Cfa * Qfa^2
  H = Cva * Qva^(2/3)
  
where Qfa = fixed area node leakage rate, Qva = variable area node leakage rate,
Cfa = 1 / SUM(Co*(L/2)*Ao)^2, Cva = 1 / SUM(Co*(L/2)*m)^2/3, and
SUM(x) is the summation of x over all pipes connected to the node.

In implementing this model, the pipe property "LeakArea" represents Ao in
sq. mm per 100 units of pipe length and "LeakExpan" represents m in sq. mm
per unit of pressure head.

*/
#include <stdlib.h>
#include <math.h>

#include "types.h"
#include "funcs.h"

// Exported functions (declared in funcs.h)
//int     openleakage(Project *);
//void    closeleakage(Project *);
//double  findlinkleakage(Project *, int);
//void    leakagecoeffs(Project *);
//double  leakageflowchange(Project *, int);
//int     leakagehasconverged(Project *);

// Local functions
static int   check_for_leakage(Project *pr);
static int   create_leakage_objects(Project *pr);
static void  convert_pipe_to_node_leakage(Project *pr);
static void  init_node_leakage(Project *pr);
static int   leakage_headloss(Project* pr, int i, double *hfa,
             double *gfa, double *hva, double *gva);
static void  eval_leak_headloss(double RQtol, double q, double c,
             double n, double *hloss, double *hgrad);
static void  add_lower_barrier(double q, double *hloss, double *hgrad);


int openleakage(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  returns an error code
**   Purpose: opens the pipe leakage modeling system
**-------------------------------------------------------------
*/
{
    Hydraul *hyd = &pr->hydraul;

    int err;
    
    // Check if project includes leakage
    closeleakage(pr);
    hyd->HasLeakage = check_for_leakage(pr);
    if (!hyd->HasLeakage) return 0;
    
    // Allocate memory for leakage data objects
    err = create_leakage_objects(pr);
    if (err > 0) return err;
    
    // Convert pipe leakage coeffs. to node coeffs.
    convert_pipe_to_node_leakage(pr);
    init_node_leakage(pr);
    return 0;
}

    
int check_for_leakage(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  returns TRUE if any pipes can leak, FALSE otherwise
**   Purpose: checks if any pipes can leak.
**-------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    int i;
    Slink *link;
    
    for (i = 1; i <= net->Nlinks; i++)
    {
        // Only pipes have leakage
        link = &net->Link[i];
        if (link->Type > PIPE) continue;
        if (link->LeakArea > 0.0 || link->LeakExpan > 0.0) return TRUE;
    }
    return FALSE;
}    

    
int create_leakage_objects(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  returns an error code
**   Purpose: allocates an array of node leakage objects.
**-------------------------------------------------------------
*/
{    
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    int i;
    
    hyd->Leakage = (Sleakage *)calloc(net->Njuncs + 1, sizeof(Sleakage));
    if (hyd->Leakage == NULL) return 101;
    for (i = 1; i <= net->Njuncs; i++)
    {
        hyd->Leakage[i].cfa = 0.0;
        hyd->Leakage[i].cva = 0.0;
        hyd->Leakage[i].qfa = 0.0;
        hyd->Leakage[i].qva = 0.0;
    }
    return 0;
}

void convert_pipe_to_node_leakage(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  none
**   Purpose: converts pipe leakage parameters into node leakage
**            coefficents.
**-------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    int i;
    double c_area, c_expan, c_orif, len;
    Slink *link;
    Snode *node1;
    Snode *node2;

    // Orifice coeff. with conversion from sq. mm to sq. m
    c_orif = 4.8149866 * 1.e-6;
    
    // Examine each link
    for (i = 1; i <= net->Nlinks; i++)
    {
        // Only pipes have leakage
        link = &net->Link[i];
        if (link->Type > PIPE) continue;
        
        // Ignore leakage in a pipe connecting two tanks or
        // reservoirs (since those nodes don't have demands)
        node1 = &net->Node[link->N1];
        node2 = &net->Node[link->N2];
        if (node1->Type != JUNCTION && node2->Type != JUNCTION) continue;

        // Get pipe's fixed and variable area leak coeffs.
        if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) continue;
        c_area = c_orif * link->LeakArea / SQR(MperFT);
        c_expan = c_orif * link->LeakExpan;

        // Adjust for number of 100-ft pipe sections
        len = link->Len * pr->Ucf[LENGTH] / 100.;
        if (node1->Type == JUNCTION && node2->Type == JUNCTION)
        {
            len *= 0.5;
        }
        c_area *= len;
        c_expan *= len;
        
        // Add these coeffs. to pipe's end nodes
        if (node1->Type == JUNCTION)
        {
            hyd->Leakage[link->N1].cfa += c_area;
            hyd->Leakage[link->N1].cva += c_expan;
        }
        if (node2->Type == JUNCTION)
        {
            hyd->Leakage[link->N2].cfa += c_area;
            hyd->Leakage[link->N2].cva += c_expan;
        }
    }
}    
    
void init_node_leakage(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  none
**   Purpose: initializes node leakage coeffs. and flows.
**-------------------------------------------------------------
*/
{    
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    int i;
    double c_area, c_expan;

    for (i = 1; i <= net->Njuncs; i++)
    {
        // Coeff. for fixed area leakage
        c_area = hyd->Leakage[i].cfa;
        if (c_area > 0.0)
            hyd->Leakage[i].cfa = 1.0 / (c_area * c_area);
        else
            hyd->Leakage[i].cfa = 0.0;
        
        // Coeff. for variable area leakage
        c_expan = hyd->Leakage[i].cva;
        if (c_expan > 0.0)
            hyd->Leakage[i].cva = 1.0 / pow(c_expan, 2./3.);
        else
            hyd->Leakage[i].cva = 0.0;
        
        // Initialize leakage flow to a non-zero value (as required by
        // the hydraulic solver)
        if (hyd->Leakage[i].cfa > 0.0)
            hyd->Leakage[i].qfa = 0.001;
        if (hyd->Leakage[i].cva > 0.0)
            hyd->Leakage[i].qva = 0.001;
    }
}

void closeleakage(Project *pr)
/*-------------------------------------------------------------
**   Input:   none
**   Output:  none
**   Purpose: frees memory for nodal leakage objects.
**-------------------------------------------------------------
*/
{
    Hydraul *hyd = &pr->hydraul;
    if (hyd->Leakage) free(hyd->Leakage);
    hyd->Leakage = NULL;
    hyd->HasLeakage = FALSE;
}
   
double findlinkleakage(Project *pr, int i)
/*-------------------------------------------------------------
**   Input:   i = link index
**   Output:  returns link leakage flow (cfs)
**   Purpose: computes leakage flow from link i at current
**            hydraulic solution.
**-------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    Smatrix *sm = &hyd->smatrix;
    Slink *link = &net->Link[i];
    
    int    n1, n2;
    double h1, h2, hsqrt, a, m, c, len, q1, q2;

    // Only pipes can leak
    link = &net->Link[i];
    if (link->Type > PIPE) return 0.0;
    
    // No leakage if area & expansion are 0
    if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) return 0.0;
   
    // No leakage if link's end nodes are both fixed grade
    n1 = link->N1;
    n2 = link->N2;
    if (n1 > net->Njuncs && n2 > net->Njuncs) return 0.0;
    
    // Pressure head of end nodes
    h1 = hyd->NodeHead[n1] - net->Node[n1].El;
    h1 = MAX(h1, 0.0);
    h2 = hyd->NodeHead[n2] - net->Node[n2].El;
    h2 = MAX(h2, 0.0);
    
    // Pipe leak parameters converted to feet
    a = link->LeakArea / SQR(MperFT);
    m = link->LeakExpan;
    len = link->Len * pr->Ucf[LENGTH] / 100.; // # 100 ft pipe lengths
    c = 4.8149866 * len / 2.0 * 1.e-6;
    
    // Leakage from 1st half of pipe connected to node n1
    q1 = 0.0;
    if (n1 <= net->Njuncs)
    {
        hsqrt = sqrt(h1);    
        q1 = c * (a + m * h1) * hsqrt;
    }
    
    // Leakage from 2nd half of pipe connected to node n2
    q2 = 0.0;
    if (n2 <= net->Njuncs)
    {
        hsqrt = sqrt(h2);    
        q2 = c * (a + m * h2) * hsqrt;
    }
    
    // Adjust leakage flows to account for one node being fixed grade
    if (q2 == 0.0) q1 *= 2.0;
    if (q1 == 0.0) q2 *= 2.0;
    return q1 + q2;
}    

void leakagecoeffs(Project *pr)
/*
**--------------------------------------------------------------
**   Input:   none
**   Output:  none
**   Purpose: computes coeffs. of the linearized hydraulic eqns.
**            contributed by node leakages.
**--------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    Smatrix *sm = &hyd->smatrix;

    int    i, row;
    double hfa,    // head loss producing current fixed area leakage
           gfa,    // gradient of fixed area head loss
           hva,    // head loss producing current variable area leakage
           gva;    // gradient of variable area head loss
    
    Snode* node;
    
    for (i = 1; i <= net->Njuncs; i++)
    {
        // Skip junctions that don't leak
        node = &net->Node[i];
        if (!leakage_headloss(pr, i, &hfa, &gfa, &hva, &gva)) continue;

        // Addition to matrix diagonal & r.h.s
        row = sm->Row[i];
        if (gfa > 0.0)
        {
            sm->Aii[row] += 1.0 / gfa;
            sm->F[row] += (hfa + node->El) / gfa;
        }
        if (gva > 0.0)
        {
            sm->Aii[row] += 1.0 / gva;
            sm->F[row] += (hva + node->El) / gva;
        }

        // Update node's flow excess (inflow - outflow)
        hyd->Xflow[i] -= (hyd->Leakage[i].qfa + hyd->Leakage[i].qva);
    }
}

double leakageflowchange(Project *pr, int i)
/*
**--------------------------------------------------------------
**   Input:   i = node index
**   Output:  returns change in leakage flow rate
**   Purpose: finds new leakage flow rate at a node after new
**            heads are computed by the hydraulic solver.
**--------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    double  hfa, gfa, hva, gva;  // same as defined in leakage_solvercoeffs()
    double  h, dqfa, dqva;       // pressure head, change in leakage flows
    
    // Find the head loss and gradient of the inverted leakage
    // equation for both fixed and variable area leakage at the
    // current leakage flow rates
    if (!leakage_headloss(pr, i, &hfa, &gfa, &hva, &gva)) return 0.0;
    
    // Pressure head using latest head solution
    h = hyd->NodeHead[i] - net->Node[i].El;

    // GGA flow update formula for fixed area leakage
    dqfa = 0.0;
    if (gfa > 0.0)
    {
        dqfa = (hfa - h) / gfa * hyd->RelaxFactor;
        hyd->Leakage[i].qfa -= dqfa;
    }

    // GGA flow update formula for variable area leakage
    dqva = 0.0;
    if (gva > 0.0)
    {
        dqva = (hva - h) / gva * hyd->RelaxFactor;
        hyd->Leakage[i].qva -= dqva;
    }

    // New leakage flow at the node 
    hyd->LeakageFlow[i] = hyd->Leakage[i].qfa + hyd->Leakage[i].qva;
    return dqfa + dqva;
}

int leakagehasconverged(Project *pr)
/*
**--------------------------------------------------------------
**   Input:   none
**   Output:  returns TRUE if leakage calculations converged,
**            FALSE if not
**   Purpose: checks if leakage calculations have converged.
**--------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    
    int i;
    double h, qref, qtest;
    const double QTOL = 0.0001;  // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm)

    for (i = 1; i <= net->Njuncs; i++)
    {
        // Skip junctions that don't leak
        if (hyd->Leakage[i].cfa == 0 && hyd->Leakage[i].cva == 0) continue;
        
        // Evaluate node's pressure head
        h = hyd->NodeHead[i] - net->Node[i].El;
        
        // Directly compute a reference leakage at this pressure head
        qref = 0.0;
        // Contribution from pipes with fixed area leaks
        if (hyd->Leakage[i].cfa > 0.0)
            qref = sqrt(h / hyd->Leakage[i].cfa);
        // Contribution from pipes with variable area leaks
        if (hyd->Leakage[i].cva > 0.0)
            qref += pow((h / hyd->Leakage[i].cva), 1.5);
        
        // Compare reference leakage to solution leakage
        qtest = hyd->Leakage[i].qfa + hyd->Leakage[i].qva;
        if (fabs(qref - qtest) > QTOL) return FALSE;
    }
    return TRUE;
}

int leakage_headloss(Project* pr, int i, double *hfa, double *gfa,
    double *hva, double *gva)
/*
**--------------------------------------------------------------
**   Input:   i = node index
**   Output:  hfa = fixed area leak head loss (ft)
**            gfa = gradient of fixed area head loss (ft/cfs)
**            hva = variable area leak head loss (ft)
**            gva = gradient of variable area head loss (ft/cfs)
**            returns TRUE if node has leakage, FALSE otherwise
**   Purpose: finds head loss and its gradient for a node's
**            leakage as a function of leakage flow.
**--------------------------------------------------------------
*/
{     
    Hydraul *hyd = &pr->hydraul;
    
    if (hyd->Leakage[i].cfa == 0.0 && hyd->Leakage[i].cva == 0.0) return FALSE;
    if (hyd->Leakage[i].cfa == 0.0)
    {
        *hfa = 0.0;
        *gfa = 0.0;
    }
    else
        eval_leak_headloss(hyd->RQtol, hyd->Leakage[i].qfa, hyd->Leakage[i].cfa,
                           0.5, hfa, gfa);
    if (hyd->Leakage[i].cva == 0.0)
    {
        *hva = 0.0;
        *gva = 0.0;
    }
    else
        eval_leak_headloss(hyd->RQtol, hyd->Leakage[i].qva, hyd->Leakage[i].cva,
                           1.5, hva, gva);
    return TRUE;
}

void eval_leak_headloss(double RQtol, double q, double c, double n,
    double *hloss, double *hgrad)
/*
**--------------------------------------------------------------
**   Input:   RQtol = low gradient tolerance (ft/cfs)
**            q = leakage flow rate (cfs)
**            c = leakage head loss coefficient
**            n = leakage head loss exponent
**   Output:  hloss = leakage head loss (ft)
**            hgrad = gradient of leakage head loss (ft/cfs)
**   Purpose: evaluates inverted form of leakage equation to
**            compute head loss and its gradient as a function
**            flow.
**
**   Note:  Inverted leakage equation is:
**            hloss = c * q ^ (1/n)
**--------------------------------------------------------------
*/
{
    n = 1.0 / n;
    *hgrad = n * c * pow(fabs(q), n - 1.0);
    
    // Use linear head loss function for small gradient
/*    if (*hgrad < RQtol)
    {
        *hgrad = RQtol / n;
        *hloss = (*hgrad) * q;
    }            

    // Otherwise use normal leakage head loss function
    else  */
        *hloss = (*hgrad) * q / n;
    
    // Prevent leakage from going negative
    add_lower_barrier(q, hloss, hgrad);
}

void add_lower_barrier(double q, double* hloss, double* hgrad)
/*
**--------------------------------------------------------------------
**  Input:   q = current flow rate
**  Output:  hloss = head loss value
**           hgrad = head loss gradient value
**  Purpose: adds a head loss barrier to prevent flow from falling
**           below 0.
**--------------------------------------------------------------------
*/
{
    double a = 1.e9 * q;
    double b = sqrt(a*a + 1.e-6);
    *hloss += (a - b) / 2.;
    *hgrad += (1.e9 / 2.) * ( 1.0 - a / b);
}