Class PhaseSrkCPAandersonReduced

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.PhaseSrkCPAandersonReduced
All Implemented Interfaces:
Serializable, Cloneable, PhaseCPAInterface, PhaseEosInterface, PhaseInterface, ThermodynamicConstantsInterface
public class PhaseSrkCPAandersonReduced extends PhaseSrkCPAs
Anderson-accelerated nested CPA phase solver with site symmetry reduction.

Combines two orthogonal acceleration strategies:

  1. Site symmetry reduction: groups equivalent association sites (same bonding pattern on the same component) into unique site types, reducing the inner loop dimension from n_s to p.
  2. Anderson acceleration: replaces successive substitution with Anderson mixing (depth m=3) on the reduced p-dimensional site fraction vector, achieving superlinear convergence.

The outer Halley iteration for molar volume and the volume derivative computation (initCPAMatrix(1)) also use the reduced dimension p, reducing the Hessian linear system from O(n_s^3) to O(p^3).

This is a nested-family solver: site fractions are fully converged at each volume step before computing exact volume derivatives via the implicit function theorem. It therefore avoids the coupled-family equilibrium sensitivity documented for the Broyden and fully implicit solvers.

Dimension reduction examples:

Dimension reduction for common CPA systems
System n_s (full) p (reduced) Reduction
Pure water (4C) 4 2 50%
Water + methanol (4C+2B) 6 4 33%
Water + MEG (4C+4C) 8 4 50%
NG + water + TEG (4C+4C) 8 4 50%
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.

    • ANDERSON_M

      private static final int ANDERSON_M
      Anderson mixing depth (number of history vectors).
      See Also:

    • MAX_INNER_ITERS

      private static final int MAX_INNER_ITERS
      Maximum inner (site fraction) iterations per outer step.
      See Also:

    • INNER_TOL

      private static final double INNER_TOL
      Inner loop convergence tolerance.
      See Also:

    • MAX_OUTER_ITERS

      private static final int MAX_OUTER_ITERS
      Maximum outer (volume) iterations.
      See Also:

    • OUTER_TOL

      private static final double OUTER_TOL
      Outer loop convergence tolerance.
      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.

    • 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

    • workXSiteFull

      private transient double[] workXSiteFull

    • workXTypeOuter

      private transient double[] workXTypeOuter

    • workTMolesOuter

      private transient double[] workTMolesOuter

    • workInnerKlk

      private transient double[][] workInnerKlk

    • workInnerXCurr

      private transient double[] workInnerXCurr

    • workInnerFX

      private transient double[] workInnerFX

    • workInnerGCurr

      private transient double[] workInnerGCurr

    • workInnerXNew

      private transient double[] workInnerXNew

    • workInnerGPrev

      private transient double[] workInnerGPrev

    • workInnerXPrev

      private transient double[] workInnerXPrev

    • workInnerGHist

      private transient double[][] workInnerGHist

    • workInnerXHist

      private transient double[][] workInnerXHist

    • workOuterP

      private transient int workOuterP

    • workOuterNs

      private transient int workOuterNs

    • skipDeltaUpdateInInitCPA

      private transient boolean skipDeltaUpdateInInitCPA
      When true, initCPAMatrix(1) skips the updateDeltaWithG(ns) call. Set by the Halley outer loop where delta has just been refreshed; cleared on exit.

    • callCount

      private static volatile long callCount
      Total molarVolume calls.

    • totalOuterIters

      private static volatile long totalOuterIters
      Total outer iterations across all calls.

    • totalInnerIters

      private static volatile long totalInnerIters
      Total inner iterations across all calls.

    • andersonConvergedCount

      private static volatile long andersonConvergedCount
      Calls where Anderson converged vs needed Newton fallback.

    • newtonFallbackCount

      private static volatile long newtonFallbackCount
      Fallbacks to Newton in inner loop.

  • Constructor Details

    • PhaseSrkCPAandersonReduced

      public PhaseSrkCPAandersonReduced()
      Constructor for PhaseSrkCPAandersonReduced.

  • 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

    • getCallCount

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

    • getAndersonConvergedCount

      public static long getAndersonConvergedCount()
      Get total Anderson converged count.
      Returns:
      number of inner loops converged by Anderson

    • getNewtonFallbackCount

      public static long getNewtonFallbackCount()
      Get total Newton fallback count.
      Returns:
      number of inner loops that fell back to Newton

    • clone

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

      Molar volume calculation using the Halley outer loop for volume and Anderson-accelerated successive substitution on reduced site type fractions. The Halley step uses the implicit function theorem to compute dF_CPA/dV derivatives in the reduced p-dimensional space.

      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.

    • solveXAndersonReduced

      private int solveXAndersonReduced(int p, double[] tMoles, int ns)
      Solve the reduced site fraction equations using Anderson-accelerated successive substitution.

      The fixed-point map operates on p site type fractions instead of n_s individual fractions: 

      f_alpha(X) = 1 / (1 + sum_beta m_beta * n_beta * Delta(rep_a,rep_b) * X_beta / V)

      Anderson acceleration (mixing depth m=3) is applied to this reduced map. After convergence, the p type fractions are expanded back to all n_s individual site fractions.

      Parameters:
      p - number of unique site types
      tMoles - moles per site type
      ns - total number of individual association sites
      Returns:
      number of inner iterations used

    • solveAndersonLeastSquares

      private static double[] solveAndersonLeastSquares(double[][] gMatrix, double[] gVec, int p, int histLen)
      Solve the Anderson least-squares problem: min ||g - G * gamma||^2.

      Uses the normal equations: (G^T G) gamma = G^T g. The system is at most m x m (typically 3x3) so a direct solve via Gaussian elimination is efficient and stable.

      Parameters:
      gMatrix - history of residual differences (m x p, using rows 0..histLen-1)
      gVec - current residual vector (length p)
      p - dimension of the vectors (number of site types)
      histLen - number of stored history vectors
      Returns:
      mixing coefficients gamma (length histLen)

    • 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

    • 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. Volume derivatives (FCPA, dFCPAdV, dFCPAdVdV, dFCPAdVdVdV) are computed using the type grouping, reducing the Hessian linear system from n_s to p dimensions.

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

    • solveLinearSystem

      private static void solveLinearSystem(double[][] a, double[] b, int n)
      Solve the 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

    • updateDeltaWithG

      private void updateDeltaWithG(int ns)
      Update delta = deltaNog * g for all individual sites.
      Parameters:
      ns - total number of individual sites

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