import java.util.*;
/*
  Revision History:
  B: Output formatted, tested with 4 nodes, sing convention reversed: 
     (-) consumption (outflow), (+)injection (inflow), pressure drop value in
     output corrected.
  -----------------------------------------------------------------------------
  Picard-iteration pipe network solver with support for pressure-jump links.
  Unknowns: heads at non-fixed-head nodes (junctions). Fixed-head nodes are
  boundaries. Each link provides a linearized flow relation at the current
  iterate:

  h_ij = (h_i - h_j - PJMP) / C(Q),
   where PJMP is the constant head jump (+ for pump, - for valves),
   and C(Q) is the Picard slope for the head-loss law for the link

  Head loss model (power-law):
   dH = K* |Q|^{m-1} * Q, Picard slope: C(Q) = K * |Q|^{m-1}

  Heads in [meters] of water column
  Flows in [m^3/s], K such that K * |Q|^{m-1} * Q in [m]
  PJMP in [m]: convert from Pa by PJMP [m] = PJMP [Pa] / (rho*g)
*/

public class PicardPipeNetwork {
  // Define cut-off thresholds (e.g. to prevent division by ZERO)
  // to avoid zero slope when Q ~ 0
  private static final double EPS_Q_SLOPE = 1e-10;
  // for linear solver pivots
  private static final double EPS_PIVOT   = 1e-10;
  
  // Node Class (isfixedHead, head, inflow/outflow)
  static class Node {
    final int id; 
    boolean isFixedHead;
    double head;        // known if isFixedHead, otherwise unknown solved value
    double demand;      // (-) consumption (outflow), (+)injection (inflow)
    Node(int id) { this.id = id; }
  }
  
  static abstract class Link {
    final int s_node; final int e_node; 
    double Q;  // flow (positive along 's_node' to 'e_node' )
    Link(int s_node, int e_node) {
      this.s_node = s_node; 
      this.e_node = e_node;
    }
    //Constant head jump (gain or loss) term for this link, PJMP [m]
    abstract double PJMP();
    
    /*Picard slope C(Q_k) for linearized flow relation at current iterate Q_k
      Should be strictly positive (use small epsilon if needed).
    */
    abstract double slopeC(double Qk);

    /* Compute 'initial' flow from heads using Picard linearized relation with
      current slope of the link:
      Q0 = (h_from - h_to - PJMP) / C(Q_k)
    */
    double computeQ0(double hFrom, double hTo) {
      double C = Math.max(slopeC(Q), EPS_Q_SLOPE);
      return (hFrom - hTo - PJMP()) / C;
    }
  }
  
  static class PowerLawJumpLink extends Link {
    final double K;    // power-law coefficient (>=0, 2 for Darcy)
    final double m;    // exponent m >= 1
    final double PJMP; // constant head term (m)
    PowerLawJumpLink(int s_node, int e_node, double K, double m, double PJMP) {
      super(s_node, e_node); 
      this.K = K; 
      this.m = m; 
      this.PJMP = PJMP;
    }
    @Override
    double PJMP() { 
      return PJMP;
    }
      
    @Override
    double slopeC(double Qk) {
      // pure PJMP element: give tiny slope to avoid divergence
      if (K <= 0.0) return EPS_Q_SLOPE;
      double mag = Math.max(Math.abs(Qk), EPS_Q_SLOPE);
      double exp = Math.max(m - 1.0, 0.0);
      return K * Math.pow(mag, exp);
    }
  }
  
  // Calculate Pipe Information
  static PowerLawJumpLink makeQuadraticPipe(int s_node, int e_node, double Kpipe) {
    //h = K * |Q| * Q
    return new PowerLawJumpLink(s_node, e_node, Kpipe, 2.0, 0.0);
  }
  
  static PowerLawJumpLink makePressureJump (int s_node, int e_node, double PJMP) {
    // Pure constant-head jump, no quadratic loss
    return new PowerLawJumpLink(s_node, e_node, 0.0, 2.0, PJMP);
  }
  
  static PowerLawJumpLink makeJumpWithMinorLoss(int s_node, int e_node, double
    PJMP, double Kminor, double m) {
    // General jump: PJMP + K * |Q|^{m-1} * Q
    return new PowerLawJumpLink(s_node, e_node, Kminor, m, PJMP);
  }
  
  //Network container
  static class Network {
    final List<Node> nodes = new ArrayList<>();
    final List<Link> links = new ArrayList <> () ;
    Node addNode (boolean fixedHead, double head, double demand) {
      Node n = new Node (nodes.size());
      n.isFixedHead = fixedHead;
      n.head        = head; 
      n.demand      = demand; 
      nodes.add(n); 
      return n;
    }
    Link addLink (Link branch) {
      links.add(branch); 
      return branch;
    }
  }
  
  static class PicardSolver {
    final Network net;
    // Convergence options
    double tolKCL = 1e-8;   // max flow change between iterations [m^3/s]
    double tolQ   = 1e-8;   // nodal continuity residual [m^3/s]
    int maxIter   = 500;
    double alphaQ = 0.60;   // under-relaxation on flows (0 < alpha <=1)
   
    PicardSolver(Network net) { 
      this.net = net;
    }
    
    void solve() {
      // Map unknown node indices
      int nUnknown = 0;
      int[] nodeToUnknown = new int [net. nodes. size()];
      Arrays.fill(nodeToUnknown, -1);
      for (Node n : net.nodes) {
        if (!n.isFixedHead) nodeToUnknown[n.id] = nUnknown++;
      }
      // Initialize link flows from a crude estimate using current heads
      for (Link branch : net.links) {
        Node ni = net.nodes.get(branch.s_node);
        Node nj = net.nodes.get (branch.e_node);
        branch.Q = branch.computeQ0(ni.head, nj.head);
      }
      
      for (int it = 1; it <= maxIter; it++) {
        // Assemble linear system [A] {h} = {b} for unknown heads
        double[][] A = new double[nUnknown][nUnknown];
        double[]   b = new double[nUnknown];
       
        // Zero initialize
        for (int i = 0; i < nUnknown; i++) {
          Arrays.fill(A[i], 0.0);
          b[i] = 0.0;
        }
        // KCL at each unknown node: sum(outflows) = demand
        for (Node n : net.nodes) {
          if (n.isFixedHead) continue;
          int row     = nodeToUnknown[n.id];
          double rhs  = n.demand;
          double diag = 0.0;
          
          // demand is outflow (-)
          for (Link branch : net.links) {
            int i = branch.s_node, j = branch.e_node;
            
            // Determine if link is incident to node n
            if (i == n.id || j == n.id) {
              double C = Math.max(branch.slopeC(branch.Q), EPS_Q_SLOPE);
              double PJMP = branch.PJMP();
              
              if (i == n.id){
                // Contribution of +Qij = (h_i - h_j) - PJMP)/C
                diag = diag + 1.0 / C;
                Node other = net.nodes.get(j);
                if (other.isFixedHead) {
                  rhs = rhs + (other.head) / C; // move -h_j/C to RHS
                } else {
                  int col = nodeToUnknown [other.id];
                  A[row][col] += -1.0 / C;
                }
                rhs = rhs + PJMP / C; // move +PJMP/C to RHS
              } else {
                //n is the 'to' node j. KCL uses -Q_ij: -ij = h_j - h_i + PJMP)/C
                diag = diag + 1.0 / C;
                Node other = net.nodes.get(i);

                if (other.isFixedHead) {
                  rhs = rhs + (other.head) / C;
                  // move -h_i/C to RHS (note signs below)
                } else {
                  int col = nodeToUnknown [other.id];
                  A[row][col] += -1.0 / C;
                }
                rhs = rhs + (-PJMP) / C;
                // from (h_j - h_i + PJMP)/C -> move +PJMP/C to LHS gives -PJMP/C on RHS
              }
            }
          }
          A[row][row] += diag;
          b[row]      += rhs;
        }
        // Solve for unknown heads
        double[] x = solveLinearSystem(A, b);
        // Update heads (unknown nodes)
        for (Node n : net.nodes) {
           if (!n.isFixedHead) {
            int idx = nodeToUnknown[n.id];
            n.head = x[idx];
           }
        }
        // Compute new raw flows from updated heads, then under-relax
        double maxDeltaQ = 0.0;
        for (Link branch : net.links) {
          Node ni = net.nodes.get(branch.s_node);
          Node nj = net.nodes.get(branch.e_node);
          double Qraw = branch.computeQ0(ni.head, nj.head);
          double Qnew = (1.0 - alphaQ) * branch.Q + alphaQ * Qraw;
          maxDeltaQ = Math.max(maxDeltaQ, Math.abs(Qnew - branch.Q) );
          branch.Q = Qnew;
        }

        // Compute KCL residual (using current linearized Q)
        double maxKCL = 0.0;
        for (Node n : net.nodes) {
          if (n.isFixedHead) continue;
          double residual = 0.0;
          for (Link branch : net.links) {
            // (+) consumption (outflow), (-)injection (inflow)
            if (branch.s_node == n.id) {
              residual = residual + branch.Q;
            } else if (branch.e_node == n.id) {
              residual = residual - branch.Q;
            }
          }
          // (-) consumption (outflow), (+)injection (inflow)
          double res = Math.abs(residual - n.demand);
          maxKCL = Math.max (maxKCL, res);
        }
        // Convergence check
        if (maxDeltaQ < tolQ && maxKCL < tolKCL) {
          System.out.printf(
          "\nConverged in %d iterations. max|dQ|=%.3e, maxkCL=%.3e%n",
          it, maxDeltaQ, maxKCL);
          return;
        }
        if (it % 10 == 0) {
          System.out.printf("Iter %d, max |dQ|=%.3e maxKCL=%.3e%n", 
          it, maxDeltaQ, maxKCL);
        } 
      }
      throw new RuntimeException ("Picard did not converge within maxIter");
    } // solve function
  
    // Simple dense linear solver with partial pivoting: Gaussian elimination
    private static double[] solveLinearSystem(double[][] A, double[] b) {
      int n = b.length;                // Number of nodes (rows)
      double[][] M = new double[n][n]; // Coefficient matrix
      double[] rhs = new double[n];    // Residual {b}
      for (int i = 0; i < n; i++) {
        System. arraycopy(A[i], 0, M[i], 0, n);
        rhs[i] = b[i];
      }
    
      for (int k = 0; k < n; k++) {
        // Pivot
        int piv = k;
        double max = Math.abs (M[k][k]);
        for (int i = k + 1; i < n; i++) {
          double v = Math.abs (M[i][k]);
          if (v > max) { 
            max = v; 
            piv = i;
          }
        }
        if (max < EPS_PIVOT) 
          throw new RuntimeException("Singular system encountered.");
        if (piv != k) {
          double[] tmp = M[k]; 
          M[k] = M[piv]; 
          M[piv] = tmp;
          double t2 = rhs[k];
          rhs[k] = rhs[piv]; 
          rhs[piv] = t2;
        }
        // Eliminate
        double diag = M[k][k];
        for (int j = k; j < n; j++) M[k][j] /= diag;
        rhs[k] = rhs[k] / diag;

        for (int i = k + 1; i < n; i++) {
          double f = M[i][k];
          if (f == 0.0) continue;
          for (int j = k; j < n; j++) M[i][j] -= f* M[k][j];
          rhs[i] = rhs[i] - f* rhs[k];
        }
      }
      // Back substitution
      double [] x = new double[n];
      for (int i = n - 1; i >= 0; i--) {
        double s = rhs[i];
        for (int j = i+1; j < n; j++) s -= M[i][j] * x[j];
        x[i] = s;
      } return x;
    }
  }

  public static void main(String[] args) {
    Network net = new Network();
    
    //--- Nodes (isFixed, head [m], demand [m3/s]): 0-----1-----2-----3
    Node n0 = net.addNode (true, 120.0,  0.00);
    Node n1 = net.addNode (false, 0.00, -0.001);
    Node n2 = net.addNode (false, 0.00,  0.00);
    Node n3 = net.addNode (true,  0.00,  0.00);
    
    // Branches
    // Pipe 0 and 1: quadratic pipe, Pipe 2: Pressure jump
    // Kpipe chosen so that K * |Q|*Q yields head in [m] for Q in [m^3/s]
    double Kpipe  = 8.0e6;
    double dPJ    = +5.0;   // +5 m head gain
    double kMinor = 1.0e6;  // minor loss
    double mMinor = 2.0;    // quadratic minor loss |Q|Q
    
    net.addLink(makeQuadraticPipe(n0.id, n1.id, Kpipe));
    net.addLink(makeQuadraticPipe(n1.id, n2.id, Kpipe));
    net.addLink(makeJumpWithMinorLoss(n2.id, n3.id, dPJ, kMinor, mMinor));
    
    // Initialize flows with small value: solver re-estimates anyway
    for (Link branch : net.links) branch.Q = 0.001;
    PicardSolver solver = new PicardSolver (net) ;
    solver.solve();
    
    // Print Results
    System.out.println("\n---=== Pressures at Nodes ===---===---==.");
    for (Node n : net.nodes) {
      System.out.printf("Node %5d: head = %+9.3f m %s | demand=%+8.4f m^3/s%n",
        n.id, n.head, n.isFixedHead ? "(-fixed)" : "(solved)", n.demand);
    }
    System.out.println("\n---=== Flow through Branches ===---===--*");
    System.out.println("Branch ID    : Flow Rate          | Pressure Drop   | Picard slope");
    System.out.println("-------------:--------------------|-----------------|----------------");
    for (int k = 0; k < net.links.size(); k++) {
      Link branch = net.links.get(k);
      Node ni = net.nodes.get(branch.s_node);
      Node nj = net.nodes.get(branch.e_node);
      double dp = ni.head - nj.head;
      System.out.printf(
        "Link %d (%d->%d): Q = %+8.5f m^3/s | dP = %+8.3f   | C(Qk) =%8.1e%n",
        k, branch.s_node, branch.e_node, branch.Q, dp, 
        branch.slopeC(branch.Q));
    }
  }
}