Class PhaseSrkCPABroydenImplicit

java.lang.Object
neqsim.thermo.phase.Phase
neqsim.thermo.phase.PhaseEos
neqsim.thermo.phase.PhaseSrkEos
neqsim.thermo.phase.PhaseSrkCPA
neqsim.thermo.phase.PhaseSrkCPAs
neqsim.thermo.phase.PhaseSrkCPABroydenImplicit
All Implemented Interfaces:
Serializable, Cloneable, PhaseCPAInterface, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface
public class PhaseSrkCPABroydenImplicit extends PhaseSrkCPAs
Broyden quasi-Newton implicit CPA phase solver.

Improves on the fully implicit coupled Newton-Raphson approach of Igben et al. (2026) by using Broyden's rank-1 update of the inverse Jacobian after the first full Newton step. This reduces the per-iteration cost from O(n_s^3) (Gaussian elimination) to O(n_s^2) (rank-1 update via the Sherman-Morrison formula), while maintaining superlinear convergence.

Algorithm:

  1. Iteration 1: Compute full analytic Jacobian J, invert to get H = J^{-1}, solve delta_x = -H * R
  2. Iterations 2+: Update H via Broyden's "good" formula: H_{k+1} = H_k + (dx - H_k*df) * (dx^T * H_k) / (dx^T * H_k * df)
  3. Periodic refresh: Recompute full Jacobian when convergence stalls (||R_new||/||R_old|| > 0.5 for 2 consecutive steps)

References:

  • C.G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Comput. 19 (1965) 577-593.
  • O.N. Igben et al., Fully implicit algorithm for CPA EOS, Fluid Phase Equilib. 608 (2026) 114734.
Version:
1.0
Author:
Even Solbraa
See Also:

  • Field Details

    • serialVersionUID

      private static final long serialVersionUID
      Serialization version UID.
      See Also:

    • logger

      private static final org.apache.logging.log4j.Logger logger
      Logger object for class.

    • MAX_ITERATIONS

      private static final int MAX_ITERATIONS
      Maximum iterations for the coupled solver.
      See Also:

    • CONVERGENCE_TOL

      private static final double CONVERGENCE_TOL
      Convergence tolerance for the coupled system.
      See Also:

    • MAX_REL_STEP

      private static final double MAX_REL_STEP
      Maximum relative step size per iteration.
      See Also:

    • RESTART_ALPHA

      private static final double RESTART_ALPHA
      Restart threshold parameter alpha.
      See Also:

    • MIN_NEWTON_STEPS

      private static final int MIN_NEWTON_STEPS
      Minimum iterations before allowing Broyden updates.
      See Also:

    • BROYDEN_SWITCH_TOL

      private static final double BROYDEN_SWITCH_TOL
      Residual threshold below which Broyden updates are allowed.
      See Also:

    • STALL_RATIO

      private static final double STALL_RATIO
      Stall detection: refresh Jacobian if ratio exceeds this for consecutive steps.
      See Also:

    • callCount

      private static volatile long callCount
      Total molarVolume calls.

    • totalIters

      private static volatile long totalIters
      Total iterations across all calls.

    • jacobianEvals

      private static volatile long jacobianEvals
      Total Jacobian evaluations (full Newton steps).

    • broydenUpdates

      private static volatile long broydenUpdates
      Total Broyden updates (quasi-Newton steps).

    • fallbackCount

      private static volatile long fallbackCount
      Total fallbacks to nested solver.

    • workKlk

      private transient double[][] workKlk

    • workHess

      private transient double[][] workHess

    • workKsi

      private transient double[] workKsi

    • workM

      private transient double[] workM

    • workKlkKsi

      private transient double[] workKlkKsi

    • workXV

      private transient double[] workXV

    • workNs

      private transient int workNs

  • Constructor Details

    • PhaseSrkCPABroydenImplicit

      public PhaseSrkCPABroydenImplicit()
      Constructor for PhaseSrkCPABroydenImplicit.

  • Method Details

    • resetProfileCounters

      public static void resetProfileCounters()
      Reset profiling counters.

    • getProfileSummary

      public static String getProfileSummary()
      Get profiling summary string.
      Returns:
      a summary of profiling data

    • clone

      public PhaseSrkCPABroydenImplicit clone()

      clone.

      Specified by:
      clone in interface PhaseInterface
      Overrides:
      clone in class PhaseSrkCPAs
      Returns:
      a PhaseInterface object

    • addComponent

      public void addComponent(String name, double moles, double molesInPhase, int compNumber)

      Add component to component array and update moles variables.

      Specified by:
      addComponent in interface PhaseInterface
      Overrides:
      addComponent in class PhaseSrkCPAs
      Parameters:
      name - Name of component.
      moles - Total number of moles of component.
      molesInPhase - Number of moles in phase.
      compNumber - Index number of component in phase object component array.

    • molarVolume

      public double molarVolume(double pressure, double temperature, double A, double B, PhaseType pt) throws IsNaNException, TooManyIterationsException

      molarVolume.

      Broyden quasi-Newton implicit molar volume calculation. Uses the coupled (n_s+1)-dimensional system but replaces full Jacobian evaluation with Broyden rank-1 updates after the first Newton step. The inverse Jacobian H is maintained directly and updated via the Sherman-Morrison formula.

      Specified by:
      molarVolume in interface PhaseInterface
      Overrides:
      molarVolume in class PhaseSrkCPA
      Parameters:
      pressure - a double
      temperature - a double
      A - a double
      B - a double
      pt - the PhaseType of the phase
      Returns:
      a double
      Throws:
      IsNaNException - if any.
      TooManyIterationsException - if any.

    • buildJacobian

      private void buildJacobian(double[][] jac, double[] xSite, double[] siteMoles, double totalVol, double gdv1, double zeta, double btemp, int dim, int ns)
      Build the analytic Jacobian for the coupled (n_s+1)-dimensional system.
      Parameters:
      jac - the Jacobian matrix to populate (dim x dim)
      xSite - site fraction values
      siteMoles - moles per site
      totalVol - total volume V
      gdv1 - g'(V) - 1/V
      zeta - current B/(nV)
      btemp - co-volume B
      dim - system dimension (ns+1)
      ns - number of association sites

    • invertMatrix

      private static boolean invertMatrix(double[][] a, double[][] ainv, int n)
      Invert matrix A into Ainv using Gauss-Jordan elimination with partial pivoting.
      Parameters:
      a - input matrix (will be modified)
      ainv - output inverse matrix
      n - matrix dimension
      Returns:
      true if successful, false if singular

    • broydenUpdate

      private static void broydenUpdate(double[][] invJ, double[] dxVec, double[] dfVec, int n)
      Broyden's "good" rank-1 update of the inverse Jacobian H.

      Updates H in-place using: H_new = H + (dx - H*df) * (dx^T * H) / (dx^T * H * df)

      This is the Sherman-Morrison formula applied to the secant condition. Cost: O(n^2).

      Parameters:
      invJ - inverse Jacobian matrix H (updated in-place)
      dxVec - step vector: x_{k+1} - x_k
      dfVec - residual change vector: R(x_{k+1}) - R(x_k)
      n - system dimension

    • updateDeltaWithG

      private void updateDeltaWithG()
      Update delta = deltaNog * g.

    • setXsiteOnComponents

      private void setXsiteOnComponents(double[] xSite)
      Set XA site fractions from flat array onto component objects.
      Parameters:
      xSite - array of site fractions

    • readXsiteFromComponents

      private void readXsiteFromComponents(double[] xSite)
      Read XA site fractions from component objects into flat array.
      Parameters:
      xSite - array to fill

    • ensureWorkArrays

      private void ensureWorkArrays(int ns)
      Ensure work arrays are allocated for the given association site count.
      Parameters:
      ns - number of association sites

    • initCPAMatrix

      public void initCPAMatrix(int type)

      initCPAMatrix.

      Override CPA matrix initialization for efficient volume derivative computation.

      Overrides:
      initCPAMatrix in class PhaseSrkCPA
      Parameters:
      type - 1 for volume derivatives, 2+ for temperature/composition

    • molarVolumeChangePhase

      public double molarVolumeChangePhase(double pressure, double temperature, double A, double B, PhaseType pt) throws IsNaNException, TooManyIterationsException

      molarVolumeChangePhase.

      Overrides:
      molarVolumeChangePhase in class PhaseSrkCPA
      Parameters:
      pressure - a double
      temperature - a double
      A - a double
      B - a double
      pt - the PhaseType of the phase
      Returns:
      a double
      Throws:
      IsNaNException - if any.
      TooManyIterationsException - if any.

    • solveLinearSystem

      private static boolean solveLinearSystem(double[][] a, double[] b, int n)
      Solve a linear system Ax = b in-place using Gaussian elimination with partial pivoting.
      Parameters:
      a - coefficient matrix (modified in place)
      b - right-hand side (solution on return)
      n - system dimension
      Returns:
      true if solve succeeded