Class PhaseSrkCPAfullyImplicitReduced

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.PhaseSrkCPAfullyImplicitReduced
All Implemented Interfaces:
Serializable, Cloneable, PhaseCPAInterface, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface
public class PhaseSrkCPAfullyImplicitReduced extends PhaseSrkCPAs
Fully implicit CPA phase with site type reduction.

Combines two acceleration strategies:

  • Fully implicit coupled Newton-Raphson from Igben et al. (2026): simultaneous solution of molar volume and association site fractions, eliminating inner iterations.
  • Site type reduction: groups equivalent association sites (same deltaNog row on same component) into types with multiplicities, reducing the system dimension from (n_s + 1) to (p + 1).

The Newton Jacobian is built analytically on the reduced (p+1)-dimensional system at every iteration (no Broyden approximation), solved via Gaussian elimination O(p^3). This gives both the per-iteration cost reduction of dimension reduction AND the quadratic convergence of full Newton.

Dimension reduction examples:

Dimension reduction for common CPA systems
System n_s (full) p (reduced) Jacobian speedup
Pure water (4C) 4 2 4.6x
Water + methanol (4C+2B) 6 4 2.7x
Water + MEG (4C+4C) 8 4 5.8x
NG + water + TEG 8 4 5.8x
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

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

    • CONVERGENCE_TOL

      static final double CONVERGENCE_TOL
      Convergence tolerance for the max residual.
      See Also:

    • MAX_REL_STEP

      static final double MAX_REL_STEP
      Maximum relative step size per iteration to prevent divergence.
      See Also:

    • RESTART_ALPHA

      static final double RESTART_ALPHA
      Restart threshold parameter alpha (supercritical detection).
      See Also:

    • numTypes

      private transient int numTypes
      Number of unique site types for the current system.

    • cachedNs

      private transient int cachedNs
      Cached ns used to build the current map (-1 = invalid).

    • siteToType

      private transient int[] siteToType
      Maps individual site index to type index.

    • typeRepSite

      private transient int[] typeRepSite
      Representative individual site index for each type.

    • typeMult

      private transient int[] typeMult
      Multiplicity (count of equivalent sites) for each type.

    • typeCompIdx

      private transient int[] typeCompIdx
      Component index for each type.

    • workRedSum

      private transient double[] workRedSum

    • workXSiteFullRead

      private transient double[] workXSiteFullRead

    • workXType

      private transient double[] workXType

    • workTMoles

      private transient double[] workTMoles

    • workResidual

      private transient double[] workResidual

    • workJacobian

      private transient double[][] workJacobian

    • workDx

      private transient double[] workDx

    • workNewtonP

      private transient int workNewtonP

    • workNewtonNs

      private transient int workNewtonNs

    • callCount

      long callCount
      Number of molarVolume calls.

    • totalIters

      long totalIters
      Total iterations across all calls.

    • jacobianEvals

      long jacobianEvals
      Number of Jacobian evaluations (full Newton, no Broyden).

    • fallbackCount

      long fallbackCount
      Number of solver fallbacks to the base class.

    • lastNumTypes

      int lastNumTypes
      Number of site types (last call).

    • lastFullSites

      int lastFullSites
      Number of full sites (last call).

    • 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

    • workP

      private transient int workP

  • Constructor Details

    • PhaseSrkCPAfullyImplicitReduced

      public PhaseSrkCPAfullyImplicitReduced()
      Construct a PhaseSrkCPAfullyImplicitReduced phase.

  • Method Details

    • clone

      public PhaseSrkCPAfullyImplicitReduced 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.

      Fully implicit + reduced molar volume solver. Operates on the reduced (p+1)-dimensional system where p is the number of unique association site types. Uses full Newton Jacobian at every iteration (no Broyden approximation), solved via Gaussian elimination.

      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.

    • buildSiteTypeMap

      private void buildSiteTypeMap(int ns)
      Build the site type map by grouping equivalent association sites.

      Two individual sites are equivalent if they belong to the same component and have identical deltaNog rows (same bonding pattern to all other sites). This corresponds to sites with the same charge in the CPA association scheme.

      Parameters:
      ns - total number of individual association sites

    • expandAndSetSiteFractions

      private void expandAndSetSiteFractions(double[] xType, int ns)
      Expand reduced site fraction values to all individual sites on components.
      Parameters:
      xType - reduced site fraction array (length p)
      ns - total number of individual sites

    • readXsiteFromComponents

      private void readXsiteFromComponents(double[] xSite, int ns)
      Read individual site fractions from component objects into a flat array.
      Parameters:
      xSite - array to fill (length ns)
      ns - total number of individual sites

    • buildReducedJacobian

      private void buildReducedJacobian(double[][] jac, double[] xType, double[] tMoles, double[] redSum, double totalVol, double gdv1, double zeta, double btemp, int dim, int p)
      Build the analytic Jacobian for the reduced (p+1)-dimensional system.

      Block structure:

      • J_XX[a][b]: site fraction block with multiplicity weighting
      • j_X_zeta[a]: site fraction sensitivity to volume (via zeta)
      • j_zeta_X[a]: volume equation sensitivity to site fractions
      • j_zeta_zeta: volume equation self-sensitivity
      Parameters:
      jac - Jacobian matrix to populate (dim x dim)
      xType - reduced site fraction values (length p)
      tMoles - moles per type (length p)
      redSum - precomputed reduced summation for each type (length p)
      totalVol - total volume V
      gdv1 - g'(V) - 1/V
      zeta - current B/(nV)
      btemp - co-volume B
      dim - system dimension (p+1)
      p - number of unique site types

    • updateDeltaWithG

      private void updateDeltaWithG(int ns)
      Multiply delta by g-function for all individual sites.
      Parameters:
      ns - total number of individual sites

    • ensureWorkArrays

      private void ensureWorkArrays(int p)
      Ensure reduced work arrays are allocated.
      Parameters:
      p - number of unique site types

    • initCPAMatrix

      public void initCPAMatrix(int type)

      initCPAMatrix.

      Override type 1 initialization to use reduced-dimension site type computation. Higher-order volume derivatives (FCPA, dFCPAdV, dFCPAdVdV, dFCPAdVdVdV) are computed using the type grouping, reducing linear system size from n_s to p.

      Overrides:
      initCPAMatrix in class PhaseSrkCPA
      Parameters:
      type - a int

    • solveLinearSystem

      private static boolean solveLinearSystem(double[][] a, double[] b, int n)
      Solve a linear system A*x = b in-place (b overwritten with solution) using Gaussian elimination with partial pivoting.
      Parameters:
      a - coefficient matrix (modified in-place)
      b - right-hand side vector (overwritten with solution)
      n - system dimension
      Returns:
      true if successful, false if singular

    • 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.

    • getCallCount

      public long getCallCount()
      Get the total number of molarVolume calls.
      Returns:
      call count

    • getTotalIters

      public long getTotalIters()
      Get the total Newton iterations across all calls.
      Returns:
      total iteration count

    • getJacobianEvals

      public long getJacobianEvals()
      Get the number of Jacobian evaluations.
      Returns:
      Jacobian evaluation count

    • getFallbackCount

      public long getFallbackCount()
      Get the number of solver fallbacks to the base class.
      Returns:
      fallback count

    • getLastNumTypes

      public int getLastNumTypes()
      Get the number of unique site types from the last call.
      Returns:
      last number of unique types

    • getLastFullSites

      public int getLastFullSites()
      Get the number of full sites from the last call.
      Returns:
      last full site count