/*
 ******************************************************************************
 Project:      OWA HYDRAULIC
 Version:      2.2
 Module:       hydsolver.c
 Description:  computes flows and pressures throughout a pipe network using
               Todini's Global Gradient Algorithm
 Authors:      see AUTHORS
 Copyright:    see AUTHORS
 License:      see LICENSE
 Last Updated: 07/15/2019
 ******************************************************************************
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

// Hydraulic balance error for network being analyzed
typedef struct {
    double maxheaderror;
    double maxflowerror;
    double maxflowchange;
    int    maxheadlink;
    int    maxflownode;
    int    maxflowlink;
} Hydbalance;

// Exported functions
int  hydsolve(Project *, int *, double *);

// Imported functions
extern int  linsolve(Project *,Smatrix *, int);  //(see SMATRIX.C)
extern int  valvestatus(Project *);    //(see HYDSTATUS.C)
extern int  linkstatus(Project *);     //(see HYDSTATUS.C)

// Local functions
static int    badvalve(Project *, int);
static int    pswitch(Project *);

static double newflows(Project *, Hydbalance *);
static void   newlinkflows(Project *, Hydbalance *, double *, double *);
static void   newemitterflows(Project *, Hydbalance *, double *, double *);
static void   newdemandflows(Project *, Hydbalance *, double *, double *);

static void   checkhydbalance(Project *, Hydbalance *);
static int    hasconverged(Project *, double *, Hydbalance *);
static int    pdaconverged(Project *);
static void   reporthydbal(Project *, Hydbalance *);


int  hydsolve(Project *pr, int *iter, double *relerr)
/*
**-------------------------------------------------------------------
**  Input:   none
**  Output:  *iter   = # of iterations to reach solution
**           *relerr = convergence error in solution
**           returns error code
**  Purpose: solves network nodal equations for heads and flows
**           using Todini's Gradient algorithm
**
**  Notes:   Status checks on CVs, pumps and pipes to tanks are made
**           every CheckFreq iteration, up until MaxCheck iterations
**           are reached. Status checks on control valves are made
**           every iteration if DampLimit = 0 or only when the
**           convergence error is at or below DampLimit. If DampLimit
**           is > 0 then future computed flow changes are only 60% of
**           their full value. A complete status check on all links
**           is made when convergence is achieved. If convergence is
**           not achieved in MaxIter trials and ExtraIter > 0 then
**           another ExtraIter trials are made with no status changes
**           made to any links and a warning message is generated.
**
**   This procedure calls linsolve() which appears in SMATRIX.C.
**-------------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    Smatrix *sm = &hyd->smatrix;
    Report  *rpt = &pr->report;

    int    i;                     // Node index
    int    errcode = 0;           // Node causing solution error
    int    nextcheck;             // Next status check trial
    int    maxtrials;             // Max. trials for convergence
    double newerr;                // New convergence error
    int    valveChange;           // Valve status change flag
    int    statChange;            // Non-valve status change flag
    Hydbalance hydbal;            // Hydraulic balance errors
    double fullDemand;            // Full demand for a node (cfs)

    // Initialize status checking & relaxation factor
    nextcheck = hyd->CheckFreq;
    hyd->RelaxFactor = 1.0;
    
    // Initialize convergence criteria and PDA results
    hydbal.maxheaderror = 0.0;
    hydbal.maxflowchange = 0.0;
    hyd->DeficientNodes = 0;
    hyd->DemandReduction = 0.0;

    // Repeat iterations until convergence or trial limit is exceeded.
    // (ExtraIter used to increase trials in case of status cycling.)
    if (rpt->Statflag == FULL) writerelerr(pr, 0, 0);
    maxtrials = hyd->MaxIter;
    if (hyd->ExtraIter > 0) maxtrials += hyd->ExtraIter;
    *iter = 1;
    while (*iter <= maxtrials)
    {
        // Compute coefficient matrices A & F and solve A*H = F
        // where H = heads, A = Jacobian coeffs. derived from
        // head loss gradients, & F = flow correction terms.
        // Solution for H is returned in F from call to linsolve().

        headlosscoeffs(pr);
        matrixcoeffs(pr);
        errcode = linsolve(pr,sm, net->Njuncs);

        // Matrix ill-conditioning problem - if control valve causing problem,
        // fix its status & continue, otherwise quit with no solution.
        if (errcode > 0)
        {
            if (badvalve(pr, sm->Order[errcode])) continue;
            else break;
        }

        // Update current solution.
        // (Row[i] = row of solution matrix corresponding to node i)
        for (i = 1; i <= net->Njuncs; i++)
        {
            hyd->NodeHead[i] = sm->F[sm->Row[i]];   // Update heads
        }
        newerr = newflows(pr, &hydbal);             // Update flows
        *relerr = newerr;

        // Write convergence error to status report if called for
        if (rpt->Statflag == FULL)
        {
            writerelerr(pr, *iter, *relerr);
        }

        // Apply solution damping & check for change in valve status
        hyd->RelaxFactor = 1.0;
        valveChange = FALSE;
        if (hyd->DampLimit > 0.0)
        {
            if (*relerr <= hyd->DampLimit)
            {
                hyd->RelaxFactor = 0.6;
                valveChange = valvestatus(pr);
            }
        }
        else
        {
            valveChange = valvestatus(pr);
        }

        // Check for convergence
        if (hasconverged(pr, relerr, &hydbal))
        {
            // We have convergence - quit if we are into extra iterations
            if (*iter > hyd->MaxIter) break;

            // Quit if no status changes occur
            statChange = FALSE;
            if (valveChange)    statChange = TRUE;
            if (linkstatus(pr)) statChange = TRUE;
            if (pswitch(pr))    statChange = TRUE;
            if (!statChange)    break;

            // We have a status change so continue the iterations
            nextcheck = *iter + hyd->CheckFreq;
        }

        // No convergence yet - see if its time for a periodic status
        // check  on pumps, CV's, and pipes connected to tank
        else if (*iter <= hyd->MaxCheck && *iter == nextcheck)
        {
            linkstatus(pr);
            nextcheck += hyd->CheckFreq;
        }
        (*iter)++;
    }

    // Iterations ended - report any errors.
    if (errcode > 0)
    {
        writehyderr(pr, sm->Order[errcode]);    // Ill-conditioned matrix error
        errcode = 110;
    }

    // Store actual junction outflow in NodeDemand & full demand in DemandFlow
    for (i = 1; i <= net->Njuncs; i++)
    {
        fullDemand = hyd->NodeDemand[i];
        hyd->NodeDemand[i] = hyd->DemandFlow[i] + hyd->EmitterFlow[i];
        hyd->DemandFlow[i] = fullDemand;
    }

    // Save convergence info
    hyd->RelativeError = *relerr;
    hyd->MaxHeadError = hydbal.maxheaderror;
    hyd->MaxFlowChange = hydbal.maxflowchange;
    hyd->Iterations = *iter;
    return errcode;
}


int  badvalve(Project *pr, int n)
/*
**-----------------------------------------------------------------
**  Input:   n = node index
**  Output:  returns 1 if node n belongs to an active control valve,
**           0 otherwise
**  Purpose: determines if a node belongs to an active control valve
**           whose setting causes an inconsistent set of eqns. If so,
**           the valve status is fixed open and a warning condition
**           is generated.
**-----------------------------------------------------------------
*/
{
    Network *net  = &pr->network;
    Hydraul *hyd  = &pr->hydraul;
    Report  *rpt  = &pr->report;
    Times   *time = &pr->times;

    int i, k, n1, n2;
    Slink *link;
    LinkType t;

    for (i = 1; i <= net->Nvalves; i++)
    {
        k = net->Valve[i].Link;
        link = &net->Link[k];
        n1 = link->N1;
        n2 = link->N2;
        if (n == n1 || n == n2)
        {
            t = link->Type;
            if (t == PRV || t == PSV || t == FCV)
            {
                if (hyd->LinkStatus[k] == ACTIVE)
                {
                    if (rpt->Statflag == FULL)
                    {
                        sprintf(pr->Msg, FMT61,
                                clocktime(rpt->Atime, time->Htime), link->ID);
                        writeline(pr, pr->Msg);
                    }
                    if (link->Type == FCV) hyd->LinkStatus[k] = XFCV;
                    else                   hyd->LinkStatus[k] = XPRESSURE;
                    return 1;
                }
            }
            return 0;
        }
    }
    return 0;
}


int  pswitch(Project *pr)
/*
**--------------------------------------------------------------
**  Input:   none
**  Output:  returns 1 if status of any link changes, 0 if not
**  Purpose: adjusts settings of links controlled by junction
**           pressures after a hydraulic solution is found
**--------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;
    Report  *rpt = &pr->report;

    int   i,                 // Control statement index
          k,                 // Index of link being controlled
          n,                 // Node controlling link k
          reset,             // Flag on control conditions
          change,            // Flag for status or setting change
          anychange = 0;     // Flag for 1 or more control actions
    char  s;                 // Current link status
    Slink *link;

    // Check each control statement
    for (i = 1; i <= net->Ncontrols; i++)
    {
        reset = 0;
        k = net->Control[i].Link;
        if (k <= 0) continue;

        // Determine if control based on a junction, not a tank
        n = net->Control[i].Node;
        if (n > 0 && n <= net->Njuncs)
        {
            // Determine if control conditions are satisfied
            if (net->Control[i].Type == LOWLEVEL &&
                hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol)
            {
                reset = 1;
            }
            if (net->Control[i].Type == HILEVEL &&
                hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol)
            {
                reset = 1;
            }
        }

        // Determine if control forces a status or setting change
        if (reset == 1)
        {
            link = &net->Link[k];
            change = 0;
            s = hyd->LinkStatus[k];
            if (link->Type == PIPE)
            {
                if (s != net->Control[i].Status) change = 1;
            }
            if (link->Type == PUMP)
            {
                if (hyd->LinkSetting[k] != net->Control[i].Setting) change = 1;
            }
            if (link->Type >= PRV)
            {
                if (hyd->LinkSetting[k] != net->Control[i].Setting) change = 1;
                else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status)
                {
                    change = 1;
                }
            }

            // If a change occurs, update status & setting
            if (change)
            {
                hyd->LinkStatus[k] = net->Control[i].Status;
                if (link->Type > PIPE)
                {
                    hyd->LinkSetting[k] = net->Control[i].Setting;
                }
                if (rpt->Statflag == FULL)
                {
                    writestatchange(pr, k, s, hyd->LinkStatus[k]);
                }
                anychange = 1;
            }
        }
    }
    return anychange;
}


double newflows(Project *pr, Hydbalance *hbal)
/*
**----------------------------------------------------------------
**  Input:   hbal = ptr. to hydraulic balance information
**  Output:  returns solution convergence error
**  Purpose: updates link, emitter & demand flows after new
**           nodal heads are computed.
**----------------------------------------------------------------
*/
{
    Hydraul *hyd = &pr->hydraul;

    double  dqsum,                 // Network flow change
            qsum;                  // Network total flow

    // Initialize sum of flows & corrections
    qsum = 0.0;
    dqsum = 0.0;
    hbal->maxflowchange = 0.0;
    hbal->maxflowlink = 1;
    hbal->maxflownode = -1;

    // Update flows in all real and virtual links
    newlinkflows(pr, hbal, &qsum, &dqsum);
    newemitterflows(pr, hbal, &qsum, &dqsum);
    newdemandflows(pr, hbal, &qsum, &dqsum);

    // Return ratio of total flow corrections to total flow
    if (qsum > hyd->Hacc) return (dqsum / qsum);
    else return dqsum;
}


void  newlinkflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
/*
**----------------------------------------------------------------
**  Input:   hbal = ptr. to hydraulic balance information
**           qsum = sum of current system flows
**           dqsum = sum of system flow changes
**  Output:  updates hbal, qsum and dqsum
**  Purpose: updates link flows after new nodal heads computed
**----------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    double  dh,                    /* Link head loss       */
            dq;                    /* Link flow change     */
    int     k, n, n1, n2;
    Slink   *link;

    // Initialize net inflows (i.e., demands) at fixed grade nodes
    for (n = net->Njuncs + 1; n <= net->Nnodes; n++)
    {
        hyd->NodeDemand[n] = 0.0;
    }

    // Examine each link
    for (k = 1; k <= net->Nlinks; k++)
    {
        // Get link and its end nodes
        link = &net->Link[k];
        n1 = link->N1;
        n2 = link->N2;

        // Apply flow update formula:
        //   dq = Y - P * (new head loss)
        //    P = 1 / (previous head loss gradient)
        //    Y = P * (previous head loss)
        // where P & Y were computed in hlosscoeff() in hydcoeffs.c

        dh = hyd->NodeHead[n1] - hyd->NodeHead[n2];
        dq = hyd->Y[k] - hyd->P[k] * dh;

        // Adjust flow change by the relaxation factor
        dq *= hyd->RelaxFactor;

        // Prevent flow in constant HP pumps from going negative
        if (link->Type == PUMP)
        {
            n = findpump(net, k);
            if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlow[k])
            {
                dq = hyd->LinkFlow[k] / 2.0;
            }
        }

        // Update link flow and system flow summation
        hyd->LinkFlow[k] -= dq;
        *qsum += ABS(hyd->LinkFlow[k]);
        *dqsum += ABS(dq);

        // Update identity of element with max. flow change
        if (ABS(dq) > hbal->maxflowchange)
        {
            hbal->maxflowchange = ABS(dq);
            hbal->maxflowlink = k;
            hbal->maxflownode = -1;
        }

        // Update net flows to fixed grade nodes
        if (hyd->LinkStatus[k] > CLOSED)
        {
            if (n1 > net->Njuncs) hyd->NodeDemand[n1] -= hyd->LinkFlow[k];
            if (n2 > net->Njuncs) hyd->NodeDemand[n2] += hyd->LinkFlow[k];
        }
    }
}


void newemitterflows(Project *pr, Hydbalance *hbal, double *qsum,
                     double *dqsum)
/*
**----------------------------------------------------------------
**  Input:   hbal = ptr. to hydraulic balance information
**           qsum = sum of current system flows
**           dqsum = sum of system flow changes
**  Output:  updates hbal, qsum and dqsum
**  Purpose: updates nodal emitter flows after new nodal heads computed
**----------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    int     i;
    double  hloss, hgrad, dh, dq;

    // Examine each network junction
    for (i = 1; i <= net->Njuncs; i++)
    {
        // Skip junction if it does not have an emitter
        if (net->Node[i].Ke == 0.0) continue;

        // Find emitter head loss and gradient
        emitterheadloss(pr, i, &hloss, &hgrad);

        // Find emitter flow change
        dh = hyd->NodeHead[i] - net->Node[i].El;
        dq = (hloss - dh) / hgrad;
        dq *= hyd->RelaxFactor;
        hyd->EmitterFlow[i] -= dq;

        // Update system flow summation
        *qsum += ABS(hyd->EmitterFlow[i]);
        *dqsum += ABS(dq);

        // Update identity of element with max. flow change
        if (ABS(dq) > hbal->maxflowchange)
        {
            hbal->maxflowchange = ABS(dq);
            hbal->maxflownode = i;
            hbal->maxflowlink = -1;
        }
    }
}


void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
/*
**----------------------------------------------------------------
**  Input:   hbal = ptr. to hydraulic balance information
**           qsum = sum of current system flows
**           dqsum = sum of system flow changes
**  Output:  updates hbal, qsum and dqsum
**  Purpose: updates nodal pressure dependent demand flows after
**           new nodal heads computed
**----------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    double  dp,         // pressure range over which demand can vary (ft)
            dq,         // change in demand flow (cfs)
            n,          // exponent in head loss v. demand  function
            hloss,      // current head loss through outflow junction (ft)
            hgrad,      // head loss gradient with respect to flow (ft/cfs)
            dh;         // new head loss through outflow junction (ft)
    int     i;
    
    // Get demand function parameters
    if (hyd->DemandModel == DDA) return;
    dp = MAX((hyd->Preq - hyd->Pmin), MINPDIFF);
    n = 1.0 / hyd->Pexp;

    // Examine each junction
    for (i = 1; i <= net->Njuncs; i++)
    {
        // Skip junctions with no positive demand
        if (hyd->NodeDemand[i] <= 0.0) continue;
        
        // Find change in demand flow (see hydcoeffs.c)
        demandheadloss(pr, i, dp, n, &hloss, &hgrad);
        dh = hyd->NodeHead[i] - net->Node[i].El - hyd->Pmin;
        dq = (hloss - dh) / hgrad;
        dq *= hyd->RelaxFactor;
        hyd->DemandFlow[i] -= dq;

        // Update system flow summation
        *qsum += ABS(hyd->DemandFlow[i]);
        *dqsum += ABS(dq);

        // Update identity of element with max. flow change
        if (ABS(dq) > hbal->maxflowchange)
        {
            hbal->maxflowchange = ABS(dq);
            hbal->maxflownode = i;
            hbal->maxflowlink = -1;
        }
    }
}


void  checkhydbalance(Project *pr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
**   Input:   hbal = hydraulic balance errors
**   Output:  none
**   Purpose: finds the link with the largest head imbalance
**--------------------------------------------------------------
*/
{
    Network *net = &pr->network;
    Hydraul *hyd = &pr->hydraul;

    int k, n1, n2;
    double dh, headerror, headloss;
    Slink *link;

    hbal->maxheaderror = 0.0;
    hbal->maxheadlink = 1;
    headlosscoeffs(pr);
    for (k = 1; k <= net->Nlinks; k++)
    {
        if (hyd->LinkStatus[k] <= CLOSED) continue;
        if (hyd->P[k] == 0.0) continue;
        link = &net->Link[k];
        n1 = link->N1;
        n2 = link->N2;
        dh = hyd->NodeHead[n1] - hyd->NodeHead[n2];
        headloss = hyd->Y[k] / hyd->P[k];
        headerror = ABS(dh - headloss);
        if (headerror > hbal->maxheaderror)
        {
            hbal->maxheaderror = headerror;
            hbal->maxheadlink = k;
        }
    }
}


int  hasconverged(Project *pr, double *relerr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
**   Input:   relerr = current total relative flow change
**            hbal   = current hydraulic balance errors
**   Output:  returns 1 if system has converged or 0 if not
**   Purpose: checks various criteria to see if system has
**            become hydraulically balanced
**--------------------------------------------------------------
*/
{
    Hydraul *hyd = &pr->hydraul;
    
    // Check that total relative flow change is small enough
    if (*relerr > hyd->Hacc) return 0;

    // Find largest head loss error and absolute flow change
    checkhydbalance(pr, hbal);
    if (pr->report.Statflag == FULL)
    {
        reporthydbal(pr, hbal);
    }
    
    // Check that head loss error and flow change criteria are met
    if (hyd->HeadErrorLimit > 0.0 &&
        hbal->maxheaderror > hyd->HeadErrorLimit) return 0;
    if (hyd->FlowChangeLimit > 0.0 &&
        hbal->maxflowchange > hyd->FlowChangeLimit) return 0;
        
    // Check for pressure driven analysis convergence
    if (hyd->DemandModel == PDA) return pdaconverged(pr);
    return 1;
}

int pdaconverged(Project *pr)
/*
**--------------------------------------------------------------
**   Input:   none
**   Output:  returns 1 if PDA converged, 0 if not
**   Purpose: checks if pressure driven analysis has converged
**            and updates total demand deficit
**--------------------------------------------------------------
*/
{
    Hydraul *hyd = &pr->hydraul;

    const double TOL = 0.001;
    int i, converged = 1;
    double totalDemand = 0.0, totalReduction = 0.0;
    
    hyd->DeficientNodes = 0;
    hyd->DemandReduction = 0.0;
    
    // Add up number of junctions with demand deficits
    for (i = 1; i <= pr->network.Njuncs; i++)
    {
        // Skip nodes whose required demand is non-positive
        if (hyd->NodeDemand[i] <= 0.0) continue;
        
        // Check for negative demand flow or positive demand flow at negative pressure
        if (hyd->DemandFlow[i] < -TOL) converged = 0;
        if (hyd->DemandFlow[i] > TOL &&
            hyd->NodeHead[i] - pr->network.Node[i].El - hyd->Pmin < -TOL)
            converged = 0;
            
        // Accumulate total required demand and demand deficit
        if (hyd->DemandFlow[i] + 0.0001 < hyd->NodeDemand[i])
        {
            hyd->DeficientNodes++;
            totalDemand += hyd->NodeDemand[i];
            totalReduction += hyd->NodeDemand[i] - hyd->DemandFlow[i];
        }
    }
    if (totalDemand > 0.0)
        hyd->DemandReduction = totalReduction / totalDemand * 100.0;
    return converged;
}


void  reporthydbal(Project *pr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
**   Input:   hbal   = current hydraulic balance errors
**   Output:  none
**   Purpose: identifies links with largest flow change and
**            largest head loss error.
**--------------------------------------------------------------
*/
{
    double qchange = hbal->maxflowchange * pr->Ucf[FLOW];
    double herror = hbal->maxheaderror * pr->Ucf[HEAD];
    int    qlink = hbal->maxflowlink;
    int    qnode = hbal->maxflownode;
    int    hlink = hbal->maxheadlink;
    if (qlink >= 1)
    {
        sprintf(pr->Msg, FMT66, qchange, pr->network.Link[qlink].ID);
        writeline(pr, pr->Msg);
    }
    else if (qnode >= 1)
    {
        sprintf(pr->Msg, FMT67, qchange, pr->network.Node[qnode].ID);
        writeline(pr, pr->Msg);
    }
    if (hlink >= 1)
    {
        sprintf(pr->Msg, FMT68, herror, pr->network.Link[hlink].ID);
        writeline(pr, pr->Msg);
    }
}