update gradient flash (#1179)

* update gradient flash

* update

* update test

* fix for Java 8

* update

src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java

@@ -7,153 +7,175 @@
 package neqsim.thermodynamicoperations.flashops;
 
 import Jama.Matrix;
+import neqsim.thermo.component.ComponentInterface;
 import neqsim.thermo.system.SystemInterface;
 
 /**
- * <p>
- * TPgradientFlash class.
- * </p>
- *
- * @author even solbraa
- * @version $Id: $Id
+ * TPgradientFlash class handles thermodynamic calculations with temperature and pressure gradients.
  */
 public class TPgradientFlash extends Flash {
   private static final long serialVersionUID = 1000;
-
-  SystemInterface localSystem = null;
-  SystemInterface tempSystem = null;
-  double temperature = 0.0;
-  double height = 0.0;
-  double deltaHeight = 0.0;
-  double deltaT = 0.0;
-  Matrix Jac;
-  Matrix fvec;
-  Matrix u;
-  Matrix dx;
-  Matrix uold;
+  private SystemInterface system;
+  private double temperature;
+  private double height;
+  private Matrix Jac;
+  private Matrix fvec;
+  private Matrix dx;
+  private SystemInterface localSystem;
+  private SystemInterface tempSystem;
+  private double deltaHeight;
+  private double deltaT;
 
   /**
-   * <p>
-   * Constructor for TPgradientFlash.
-   * </p>
+   * Default constructor for TPgradientFlash.
    */
   public TPgradientFlash() {}
 
   /**
-   * <p>
-   * Constructor for TPgradientFlash.
-   * </p>
+   * Constructor for TPgradientFlash with system, height, and temperature.
    *
    * @param system a {@link neqsim.thermo.system.SystemInterface} object
-   * @param height a double
-   * @param temperature a double
+   * @param height a double representing the height
+   * @param temperature a double representing the temperature
    */
   public TPgradientFlash(SystemInterface system, double height, double temperature) {
     this.system = system;
     this.temperature = temperature;
     this.height = height;
-    Jac = new Matrix(system.getPhase(0).getNumberOfComponents(),
-        system.getPhase(0).getNumberOfComponents());
-    fvec = new Matrix(system.getPhase(0).getNumberOfComponents(), 1);
+    int numComponents = system.getPhase(0).getNumberOfComponents();
+    Jac = new Matrix(numComponents + 1, numComponents + 1);
+    fvec = new Matrix(numComponents + 1, 1);
   }
 
   /**
-   * <p>
-   * Setter for the field <code>fvec</code>.
-   * </p>
+   * Sets the fvec matrix values based on the thermodynamic calculations.
    */
   public void setfvec() {
-    for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
+    double sumx = 0.0;
+    int numComponents = system.getPhase(0).getNumberOfComponents();
+    double gravity = neqsim.thermo.ThermodynamicConstantsInterface.gravity;
+    double gasConstant = neqsim.thermo.ThermodynamicConstantsInterface.R;
+
+    for (int i = 0; i < numComponents; i++) {
+      ComponentInterface component = (ComponentInterface) localSystem.getPhase(0).getComponent(i);
+      double fugacityCoeff = component.getFugacityCoefficient();
+      double componentX = component.getx();
+      double pressure = localSystem.getPressure();
+
+      // Breakdown of terms for readability
+      double logTerm = Math.log(fugacityCoeff * componentX * pressure);
+      double molarMassTerm = component.getMolarMass() * gravity * deltaHeight
+          / (gasConstant * tempSystem.getPhase(0).getTemperature());
+
       fvec.set(i, 0,
-          Math.log(localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient()
-              * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure())
+          logTerm
               - Math.log(tempSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient()
                   * tempSystem.getPhases()[0].getComponents()[i].getx() * tempSystem.getPressure())
-              - tempSystem.getPhases()[0].getComponents()[i].getMolarMass()
-                  * neqsim.thermo.ThermodynamicConstantsInterface.gravity * deltaHeight
-                  / neqsim.thermo.ThermodynamicConstantsInterface.R
-                  / tempSystem.getPhase(0).getTemperature()
-              + tempSystem.getPhases()[0].getComponents()[i].getMolarMass()
-                  * (tempSystem.getPhases()[0].getEnthalpy()
-                      / tempSystem.getPhases()[0].getNumberOfMolesInPhase()
-                      / tempSystem.getPhase(0).getMolarMass()
-                      - tempSystem.getPhases()[0].getComponents()[i]
-                          .getEnthalpy(tempSystem.getPhase(0).getTemperature())
-                          / tempSystem.getPhases()[0].getComponent(i).getNumberOfMolesInPhase()
-                          / tempSystem.getPhase(0).getComponent(i).getMolarMass())
-                  * deltaT / tempSystem.getPhase(0).getTemperature()
-                  / neqsim.thermo.ThermodynamicConstantsInterface.R
-                  / tempSystem.getPhase(0).getTemperature());
+              - molarMassTerm);
+      sumx += componentX;
     }
+
+    fvec.set(numComponents, 0, sumx - 1.0);
   }
 
   /**
-   * <p>
-   * setJac.
-   * </p>
+   * Sets the Jacobian (Jac) matrix values based on the thermodynamic calculations.
    */
   public void setJac() {
     Jac.timesEquals(0.0);
-    double dij = 0.0;
-
-    double tempJ = 0.0;
-
-    for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
-      for (int j = 0; j < system.getPhase(0).getNumberOfComponents(); j++) {
-        dij = i == j ? 1.0 : 0.0; // Kroneckers delta
-        tempJ = 1.0
-            / (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient()
-                * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure())
-            * (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() * dij
-                * localSystem.getPressure()
-                + localSystem.getPhases()[0].getComponents()[i].getdfugdx(j)
-                    * localSystem.getPhases()[0].getComponents()[i].getx()
-                    * localSystem.getPressure());
+    int numComponents = system.getPhase(0).getNumberOfComponents();
+
+    for (int i = 0; i < numComponents; i++) {
+      for (int j = 0; j < numComponents; j++) {
+        double dij = (i == j) ? 1.0 : 0.0; // Kronecker delta
+        double fugacityCoeff =
+            localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient();
+        double componentX = localSystem.getPhases()[0].getComponents()[i].getx();
+        double pressure = localSystem.getPressure();
+
+        double tempJ = 1.0 / (fugacityCoeff * componentX * pressure)
+            * (fugacityCoeff * dij * pressure
+                + localSystem.getPhases()[0].getComponents()[i].getdfugdx(j) * componentX
+                    * pressure);
         Jac.set(i, j, tempJ);
       }
     }
+
+    // Set the last row of Jac
+    for (int j = 0; j < numComponents; j++) {
+      Jac.set(numComponents, j, 1.0);
+    }
+
+    for (int i = 0; i < numComponents; i++) {
+      double fugacityCoeff = localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient();
+      double componentX = localSystem.getPhases()[0].getComponents()[i].getx();
+      double pressure = localSystem.getPressure();
+
+      double tempJ = 1.0 / (fugacityCoeff * componentX * pressure)
+          * (localSystem.getPhases()[0].getComponents()[i].getdfugdp() * componentX * pressure
+              + fugacityCoeff * componentX);
+      Jac.set(i, numComponents, tempJ);
+    }
+
+    Jac.set(numComponents, numComponents, 0.0);
   }
 
   /**
-   * <p>
-   * setNewX.
-   * </p>
+   * Updates the composition and pressure in the local system based on calculated adjustments.
    */
   public void setNewX() {
-    for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
-      localSystem.getPhase(0).getComponent(i)
-          .setx(localSystem.getPhase(0).getComponent(i).getx() - 0.5 * dx.get(i, 0));
+    int numComponents = system.getPhase(0).getNumberOfComponents();
+    double relaxationFactor = 0.8; // Relaxation factor for numerical stability
+
+    for (int i = 0; i < numComponents; i++) {
+      ComponentInterface component = (ComponentInterface) localSystem.getPhase(0).getComponent(i);
+      double newX = component.getx() - relaxationFactor * dx.get(i, 0);
+      component.setx(newX);
     }
-    localSystem.getPhase(0).normalize();
+
+    // Update system pressure
+    double newPressure = localSystem.getPressure() - relaxationFactor * dx.get(numComponents, 0);
+    localSystem.setPressure(newPressure);
+
+    // Optional: Normalize phase composition
+    // localSystem.getPhase(0).normalize();
   }
 
   /** {@inheritDoc} */
   @Override
   public void run() {
+    // Clone the system for temporary calculations
     tempSystem = system.clone();
     tempSystem.init(0);
     tempSystem.init(3);
 
+    // Clone the system for local calculations
     localSystem = system.clone();
-    // localSystem.setPressure(height*9.81*height);
 
+    // Initialize delta values for temperature and height
     deltaT = (temperature - system.getTemperature()) / 20.0;
     deltaHeight = height / 20.0;
 
-    for (int i = 0; i < 20; i++) {
+    for (int step = 0; step < 20; step++) {
+      // Incrementally adjust the temperature in the local system
       localSystem.setTemperature(localSystem.getTemperature() + deltaT);
-      for (int o = 0; o < 15; o++) {
-        localSystem.init(3);
-        setfvec();
-        setJac();
-        dx = Jac.solve(fvec);
-        dx.print(10, 10);
-        setNewX();
-      }
-      // localSystem.display();
 
+      int iter = 0;
+      double tolerance = 1e-10; // Convergence criterion
+      int maxIterations = 50; // Maximum allowable iterations
+
+      do {
+        iter++;
+        localSystem.init(3); // Initialize system properties
+        setfvec(); // Calculate the function vector
+        setJac(); // Calculate the Jacobian matrix
+        dx = Jac.solve(fvec); // Solve for updates (dx)
+        setNewX(); // Update composition and pressure
+      } while (dx.norm2() > tolerance && iter < maxIterations);
+
+      // Clone the updated local system into tempSystem
       tempSystem = localSystem.clone();
-      tempSystem.init(3);
+      tempSystem.init(3); // Re-initialize tempSystem for the next step
     }
   }
 
@@ -162,10 +184,5 @@ public void run() {
   public SystemInterface getThermoSystem() {
     return localSystem;
   }
-
-  /** {@inheritDoc} */
-  @Override
-  public org.jfree.chart.JFreeChart getJFreeChart(String name) {
-    return null;
-  }
 }
+

src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java

@@ -0,0 +1,75 @@
+package neqsim.thermodynamicoperations.flashops;
+
+import static org.junit.jupiter.api.Assertions.assertEquals;
+import java.util.ArrayList;
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+import org.junit.jupiter.api.Test;
+import neqsim.thermo.system.SystemInterface;
+import neqsim.thermo.system.SystemSrkEos;
+import neqsim.thermodynamicoperations.ThermodynamicOperations;
+
+public class TPgradientFlashTest {
+  static Logger logger = LogManager.getLogger(TPgradientFlashTest.class);
+
+  @Test
+  void testRun() {
+    SystemInterface testSystem = new SystemSrkEos(345, 80.0);
+
+    ThermodynamicOperations testOps = new ThermodynamicOperations(testSystem);
+    testSystem.addComponent("hydrogen", 1.0);
+    testSystem.addComponent("methane", 8.0);
+    testSystem.addComponent("propane", 1.0);
+
+    testSystem.createDatabase(true);
+    testSystem.setMixingRule(2);
+
+    SystemInterface newSystem = null;
+    try {
+      newSystem = testOps.TPgradientFlash(1000, 355).phaseToSystem(0);
+      // newSystem.prettyPrint();
+    } catch (Exception ex) {
+      logger.error(ex.getMessage(), ex);
+    }
+
+    assertEquals(0.095513700959, newSystem.getComponent("hydrogen").getx(), 1e-4);
+  }
+
+  @Test
+  void testGradient() {
+
+    double depth = 1000.0;
+    double temperature = 273.15 + 70;
+    double pressure = 80.0;
+
+    ArrayList<Double> x_h2 = new ArrayList<>();
+    ArrayList<Double> p_depth = new ArrayList<>();
+
+    SystemInterface testSystem = new SystemSrkEos(temperature, pressure);
+    testSystem.addComponent("hydrogen", 1.0);
+    testSystem.addComponent("methane", 8.0);
+    testSystem.addComponent("propane", 1.0);
+
+    testSystem.createDatabase(true);
+    testSystem.setMixingRule(2);
+
+    ThermodynamicOperations testOps = new ThermodynamicOperations(testSystem);
+    SystemInterface newSystem = null;
+    double deltaHeight = 100.0;
+    double deltaT = 0.5;
+    for (int i = 0; i < 10; i++) {
+      try {
+        newSystem = testOps.TPgradientFlash(i * deltaHeight, temperature + i * deltaT);
+      } catch (Exception ex) {
+        logger.error(ex.getMessage(), ex);
+      }
+      x_h2.add(newSystem.getComponent("hydrogen").getx());
+      p_depth.add(newSystem.getPressure());
+      // System.out
+      // .println(newSystem.getComponent("hydrogen").getx() + " " + newSystem.getPressure());
+    }
+    assertEquals(0.0964169380341, x_h2.get(6), 1e-4);
+
+  }
+
+}