Class GibbsReactor

java.lang.Object
neqsim.util.NamedBaseClass
neqsim.process.SimulationBaseClass
neqsim.process.equipment.ProcessEquipmentBaseClass
neqsim.process.equipment.TwoPortEquipment
neqsim.process.equipment.reactor.GibbsReactor
All Implemented Interfaces:
Serializable, Runnable, ProcessEquipmentInterface, TwoPortInterface, ProcessElementInterface, SimulationInterface, NamedInterface
public class GibbsReactor extends TwoPortEquipment
Gibbs reactor for chemical equilibrium calculations using Gibbs free energy minimization.

This reactor computes chemical equilibrium compositions by minimizing the total Gibbs free energy of the system subject to elemental mass balance constraints. The implementation uses the Newton-Raphson method with Lagrange multipliers to solve the constrained optimization problem.

Key Features

  • Supports both isothermal and adiabatic operation modes
  • Handles multi-component systems with element balance constraints
  • Allows marking specific components as inert (excluded from reactions)
  • Provides detailed convergence diagnostics and mass balance verification

Algorithm

The reactor minimizes the objective function: 

G = Σ nᵢ(μᵢ⁰ + RT ln(φᵢyᵢP)) - Σ λⱼ(Σ aᵢⱼnᵢ - bⱼ)

where nᵢ are molar amounts, μᵢ⁰ is standard chemical potential, φᵢ is fugacity coefficient, yᵢ is mole fraction, λⱼ are Lagrange multipliers, aᵢⱼ are stoichiometric coefficients, and bⱼ are element totals.

Usage Example

GibbsReactor reactor = new GibbsReactor("reactor", inletStream);
reactor.setEnergyMode(EnergyMode.ISOTHERMAL);
reactor.setMaxIterations(5000);
reactor.setConvergenceTolerance(1e-6);
reactor.setDampingComposition(0.01);
reactor.run();

if (reactor.hasConverged()) {
  Stream outlet = reactor.getOutletStream();
  double enthalpyChange = reactor.getEnthalpyOfReactions();
}
Version:
$Id: $Id
Author:
Sviatoslav Eroshkin
See Also:

  • Field Details

    • tempFugacitySystem

      private transient ThreadLocal<SystemInterface> tempFugacitySystem

    • serialVersionUID

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

    • REFERENCE_TEMPERATURE

      private static final double REFERENCE_TEMPERATURE
      Reference temperature in Kelvin for thermodynamic calculations (298.15 K = 25°C).
      See Also:

    • R_KJ

      private static final double R_KJ
      Universal gas constant in kJ/(mol·K).
      See Also:

    • MIN_MOLES

      private static final double MIN_MOLES
      Minimum mole amount to prevent numerical issues with log(0).
      See Also:

    • MIN_JACOBIAN_MOLES

      private static final double MIN_JACOBIAN_MOLES
      Minimum value for Jacobian calculations to avoid division by zero.
      See Also:

    • ELEMENT_ZERO_THRESHOLD

      private static final double ELEMENT_ZERO_THRESHOLD
      Small threshold for detecting zero element coefficients.
      See Also:

    • logger

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

    • ION_NAME_PATTERN

      private static final Pattern ION_NAME_PATTERN
      Pattern to identify ionic species (names ending with + or -).

    • energyMode

      private GibbsReactor.EnergyMode energyMode

    • componentAliases

      private static final Map<String,String> componentAliases

    • method

      private String method

    • useAllDatabaseSpecies

      private boolean useAllDatabaseSpecies

    • gibbsDatabase

      private transient List<GibbsReactor.GibbsComponent> gibbsDatabase

    • componentMap

      private transient Map<String, GibbsReactor.GibbsComponent> componentMap

    • lambda

      private double[] lambda

    • lagrangeContributions

      private Map<String,Double> lagrangeContributions

    • elementNames

      private String[] elementNames

    • processedComponents

      private List<String> processedComponents

    • objectiveFunctionValues

      private Map<String,Double> objectiveFunctionValues

    • inertComponents

      private final Set<String> inertComponents

    • initialMoles

      private Map<String,Double> initialMoles

    • finalMoles

      private Map<String,Double> finalMoles

    • elementMoleBalanceIn

      private double[] elementMoleBalanceIn

    • elementMoleBalanceOut

      private double[] elementMoleBalanceOut

    • elementMoleBalanceDiff

      private double[] elementMoleBalanceDiff

    • inlet_mole

      private List<Double> inlet_mole

    • outlet_mole

      private List<Double> outlet_mole

    • processedComponentIndexMap

      private Map<String,Integer> processedComponentIndexMap

    • objectiveMinimizationVector

      private double[] objectiveMinimizationVector

    • objectiveMinimizationVectorLabels

      private List<String> objectiveMinimizationVectorLabels

    • jacobianMatrix

      private double[][] jacobianMatrix

    • jacobianInverse

      private double[][] jacobianInverse

    • jacobianRowLabels

      private List<String> jacobianRowLabels

    • jacobianColLabels

      private List<String> jacobianColLabels

    • variableComponents

      private List<String> variableComponents

    • maxIterations

      private int maxIterations

    • convergenceTolerance

      private double convergenceTolerance

    • dampingComposition

      private double dampingComposition

    • actualIterations

      private int actualIterations

    • converged

      private boolean converged

    • finalConvergenceError

      private double finalConvergenceError

    • minIterations

      private int minIterations

    • useAdaptiveStepSize

      private boolean useAdaptiveStepSize

    • useArmijoLineSearch

      private boolean useArmijoLineSearch
      Enable Armijo backtracking line search for guaranteed Gibbs energy decrease.

      When enabled, the Newton step is scaled by alpha in (0, 1] such that the Armijo sufficient decrease condition is satisfied: G(n + alpha * dn) <= G(n) + c1 * alpha * grad^T * dn. This replaces fixed damping with an adaptive, globally convergent strategy.

      Reference: Nocedal & Wright, Numerical Optimization (2006), Algorithm 3.1.

    • armijoC1

      private double armijoC1
      Armijo sufficient decrease parameter c1 in (0, 1). Typical value 1e-4.

    • armijoRho

      private double armijoRho
      Backtracking contraction factor rho in (0, 1). Step is multiplied by rho each backtrack.

    • armijoMaxBacktracks

      private int armijoMaxBacktracks
      Maximum number of backtracking steps in Armijo line search.

    • useRegularization

      private boolean useRegularization
      Enable Tikhonov regularization for ill-conditioned Jacobians.

      When the condition number of the Jacobian exceeds regularizationThreshold, a regularization term tau * I is added to the Hessian block: H_reg = H + tau * I. This prevents divergence near phase boundaries where the Hessian becomes singular.

      Reference: Tikhonov, A. N. (1963). Dokl. Akad. Nauk SSSR, 151, 501-504.

    • regularizationThreshold

      private double regularizationThreshold
      Condition number threshold above which regularization is applied.

    • regularizationTau

      private double regularizationTau
      Tikhonov regularization parameter tau. Added to Hessian diagonal.

    • conditionNumberHistory

      private transient List<Double> conditionNumberHistory
      Record of condition numbers during iterations for diagnostics.

    • gibbsEnergyHistory

      private transient List<Double> gibbsEnergyHistory
      Record of Gibbs energy values during iterations for diagnostics.

    • stepSizeHistory

      private transient List<Double> stepSizeHistory
      Record of step sizes used during iterations for diagnostics.

    • elementBalanceErrorHistory

      private transient List<Double> elementBalanceErrorHistory
      Record of element balance errors during iterations for diagnostics.

    • columnScaleFactors

      private transient double[] columnScaleFactors
      Column scaling factors from the NASA CEA log-mole transformation applied to the Jacobian. Each entry j stores the factor by which Jacobian column j was multiplied. The Newton step deltaX_scaled must be divided by these factors to obtain the true deltaX in mole units.

    • MAX_STEP_LIMIT

      private static final double MAX_STEP_LIMIT
      Maximum relative change per step in moles, following NASA CEA step limiting (Gordon & McBride, 1994, NASA RP-1311). Prevents overshooting in Newton-Raphson iterations: max |alpha * deltaX[i]| / n_i &lt; MAX_STEP_LIMIT.

    • system

      private SystemInterface system

    • inletEnthalpy

      private double inletEnthalpy

    • outletEnthalpy

      private double outletEnthalpy

    • dT

      private double dT

    • tempUpdateIter

      private int tempUpdateIter

    • enthalpyOld

      double enthalpyOld

    • enthalpyOfReactions

      private double enthalpyOfReactions

    • temperatureChange

      private double temperatureChange

    • deltaNorm

      private double deltaNorm

    • GOLD

      private double GOLD

    • G

      private double G

    • dG

      private double dG

    • formulaToComponent

      private static final Map<String,String> formulaToComponent

  • Constructor Details

    • GibbsReactor

      public GibbsReactor(String name)
      Constructor for GibbsReactor.
      Parameters:
      name - Name of GibbsReactor

    • GibbsReactor

      public GibbsReactor(String name, StreamInterface stream)
      Constructor for GibbsReactor.
      Parameters:
      name - Name of GibbsReactor
      stream - Stream to set as inlet Stream. A clone of stream is set as outlet stream.

  • Method Details

    • getMassBalanceError

      public double getMassBalanceError()
      Get the relative mass balance error as a percentage. Computed as 100 * |massIn - massOut| / massIn.
      Returns:
      relative mass balance error in percent (e.g., 0.001 means 0.001%)

    • getMassBalanceConverged

      public boolean getMassBalanceConverged()
      Returns true if the relative mass balance error is less than 0.001%.
      Returns:
      true if mass balance is converged, false otherwise

    • ensureTempFugacitySystemInitialized

      private void ensureTempFugacitySystemInitialized()
      Ensures tempFugacitySystem is initialized (handles deserialization case).

    • getEnthalpyOfReactions

      public double getEnthalpyOfReactions()
      Get the cumulative enthalpy of reaction (sum of dH for all iterations).
      Returns:
      enthalpyOfReactions in kJ

    • getTemperatureChange

      public double getTemperatureChange()
      Get the cumulative temperature change during the reaction (sum of dT for all iterations).
      Returns:
      temperatureChange in K

    • getPower

      public double getPower()
      Get the reactor power (default: W, negative enthalpyOfReactions*1000).
      Returns:
      Power in Watts (W)

    • getPower

      public double getPower(String unit)
      Get the reactor power in the specified unit ("W", "kW", or "MW").
      Parameters:
      unit - Power unit: "W", "kW", or "MW" (case-insensitive, default is W)
      Returns:
      Power in the specified unit

    • calculateMixtureEnthalpy

      public double calculateMixtureEnthalpy(List<String> componentNames, List<Double> n, double T, Map<String, GibbsReactor.GibbsComponent> componentMap)
      Calculate the total enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      T - Temperature in K
      componentMap - Map from component name (lowercase) to GibbsComponent
      Returns:
      Total enthalpy (kJ)

    • calculateMixtureGibbsEnergy

      public double calculateMixtureGibbsEnergy(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap, double T)
      Calculate the total Gibbs energy of a mixture: sum_i n_i * gibbs_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      T - Temperature in K
      Returns:
      Total Gibbs energy (kJ)

    • getMixtureEnthalpy

      public double getMixtureEnthalpy()
      Get the last calculated mixture enthalpy (kJ).
      Returns:
      a double

    • getMixtureGibbsEnergy

      public double getMixtureGibbsEnergy()
      Get the last calculated mixture Gibbs energy (kJ).
      Returns:
      a double

    • calculateMixtureEnthalpyStandard

      public double calculateMixtureEnthalpyStandard(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap)
      Calculate the total standard enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      Returns:
      Total enthalpy (kJ)

    • calculateMixtureEnthalpy

      public double calculateMixtureEnthalpy(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap, double T)
      Calculate the total standard enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      T - a double
      Returns:
      Total enthalpy (kJ)

    • setEnergyMode

      public void setEnergyMode(GibbsReactor.EnergyMode mode)
      Set the energy mode of the reactor (isothermal or adiabatic).
      Parameters:
      mode - EnergyMode.ISOTHERMAL or EnergyMode.ADIABATIC

    • setEnergyMode

      public void setEnergyMode(String mode)
      Set the energy mode of the reactor using a string (case-insensitive). Accepts "adiabatic" or "isothermal" (case-insensitive).
      Parameters:
      mode - String representing the energy mode
      Throws:
      IllegalArgumentException - if the mode is not recognized

    • getEnergyMode

      public GibbsReactor.EnergyMode getEnergyMode()
      Gets the current energy mode of the reactor.
      Returns:
      the current EnergyMode (ISOTHERMAL or ADIABATIC)

    • setComponentAsInert

      public void setComponentAsInert(String componentName)
      Mark a component as inert by name.
      Parameters:
      componentName - component name as in the thermo system

    • setComponentAsInert

      public void setComponentAsInert(int index)
      Mark a component as inert by index in inlet system.
      Parameters:
      index - component index in inlet system

    • isComponentInert

      public boolean isComponentInert(String componentName)
      Check if a component is inert.
      Parameters:
      componentName - the name of the component to check
      Returns:
      true if the component is inert, false otherwise

    • isIonicComponent

      private boolean isIonicComponent(String moleculeName)

    • loadGibbsDatabase

      private void loadGibbsDatabase()
      Load the Gibbs reaction database from resources.

    • setUseAllDatabaseSpecies

      public void setUseAllDatabaseSpecies(boolean useAllDatabaseSpecies)
      Set whether to use all database species or only species in the system.
      Parameters:
      useAllDatabaseSpecies - true to use all database species, false to use only system species

    • getUseAllDatabaseSpecies

      public boolean getUseAllDatabaseSpecies()
      Get whether using all database species or only species in the system.
      Returns:
      true if using all database species, false if using only system species

    • getMethod

      public String getMethod()
      Get the method used for Gibbs minimization.
      Returns:
      The method name

    • setMethod

      public void setMethod(String method)
      Set the method used for Gibbs minimization.
      Parameters:
      method - The method name

    • getLagrangeContributions

      public Map<String,Double> getLagrangeContributions()
      Get the Lagrange contributions (legacy method).
      Returns:
      Map of Lagrange contributions

    • run

      public void run(UUID id)

      In this method all thermodynamic and unit operations will be calculated in a steady state calculation.

      Parameters:
      id - UUID

    • performGibbsMinimization

      private void performGibbsMinimization(SystemInterface system)
      Perform Gibbs free energy minimization.
      Parameters:
      system - The thermodynamic system

    • calculateElementMoleBalance

      private void calculateElementMoleBalance(SystemInterface system, double[] elementBalance, boolean isInput)
      Calculate element mole balance for a system.
      Parameters:
      system - The thermodynamic system
      elementBalance - Array to store the element balance
      isInput - true if this is input balance, false if output balance

    • calculateObjectiveFunctionValues

      private void calculateObjectiveFunctionValues(SystemInterface system)
      Calculate objective function values for each component.
      Parameters:
      system - The thermodynamic system

    • setLagrangeMultiplier

      public void setLagrangeMultiplier(int index, double value)
      Set a Lagrange multiplier value.
      Parameters:
      index - The element index (0=O, 1=N, 2=C, 3=H, 4=S, 5=Ar)
      value - The Lagrange multiplier value

    • getLagrangianMultipliers

      public double[] getLagrangianMultipliers()
      Get the Lagrange multiplier values.
      Returns:
      Array of Lagrange multiplier values [O, N, C, H, S, Ar]

    • getElementNames

      public String[] getElementNames()
      Get the element names.
      Returns:
      Array of element names

    • getLagrangeMultiplierContributions

      public Map<String, Map<String,Double>> getLagrangeMultiplierContributions()
      Get the Lagrange multiplier contributions for each component.
      Returns:
      Map with component names and their Lagrange multiplier contributions

    • getObjectiveFunctionValues

      public Map<String,Double> getObjectiveFunctionValues()
      Get the objective function values for each component.
      Returns:
      Map with component names and their objective function values

    • getElementMoleBalanceIn

      public double[] getElementMoleBalanceIn()
      Get the element mole balance for input streams.
      Returns:
      Array of total moles for each element in input [O, N, C, H, S, Ar]

    • getElementMoleBalanceOut

      public double[] getElementMoleBalanceOut()
      Get the element mole balance for output streams.
      Returns:
      Array of total moles for each element in output [O, N, C, H, S, Ar]

    • getElementMoleBalanceDiff

      public double[] getElementMoleBalanceDiff()
      Get the element mole balance difference (in - out).
      Returns:
      Array of mole differences for each element [O, N, C, H, S, Ar]

    • getDetailedMoleBalance

      public Map<String, Map<String,Double>> getDetailedMoleBalance()
      Get detailed mole balance information for each component.
      Returns:
      Map with component names and their element contributions to mole balance

    • calculateObjectiveMinimizationVector

      private void calculateObjectiveMinimizationVector()
      Calculate the objective minimization vector. This vector contains the F values for each component and the mass balance constraints. The system is in equilibrium when this vector is zero. Only includes elements that are actually present in the system.

    • getObjectiveMinimizationVector

      public double[] getObjectiveMinimizationVector()
      Get the objective minimization vector.
      Returns:
      The objective minimization vector

    • getObjectiveVectorForVariables

      private double[] getObjectiveVectorForVariables()
      Build objective vector matching the optimization variables (variableComponents) and active element balances. This excludes inert components from the variable part but keeps balances.
      Returns:
      the objective vector for optimization variables

    • getObjectiveMinimizationVectorLabels

      public List<String> getObjectiveMinimizationVectorLabels()
      Get the labels for the objective minimization vector.
      Returns:
      List of labels for each element in the objective minimization vector

    • calculateJacobian

      private void calculateJacobian()
      Calculate the Jacobian matrix for the Newton-Raphson method. The Jacobian represents the derivatives of the objective function with respect to the variables. Only includes elements that are actually present in the system to avoid singular matrices.

    • getJacobianMatrix

      public double[][] getJacobianMatrix()
      Get the Jacobian matrix.
      Returns:
      The Jacobian matrix

    • getJacobianRowLabels

      public List<String> getJacobianRowLabels()
      Get the Jacobian row labels.
      Returns:
      List of row labels for the Jacobian matrix

    • getJacobianColLabels

      public List<String> getJacobianColLabels()
      Get the Jacobian column labels.
      Returns:
      List of column labels for the Jacobian matrix

    • getJacobianInverse

      public double[][] getJacobianInverse()
      Get the Jacobian inverse matrix.
      Returns:
      The Jacobian inverse matrix, or null if it couldn't be calculated

    • calculateJacobianInverse

      private double[][] calculateJacobianInverse()
      Calculate the inverse of the Jacobian matrix using JAMA Matrix library.
      Returns:
      Inverse matrix, or null if matrix is singular

    • solveNewtonSystem

      private double[] solveNewtonSystem(double[] objectiveVector)
      Solve the Newton-Raphson linear system J_scaled * deltaX_scaled = -F using LU decomposition, then unscale using the column scaling factors from the NASA CEA log-mole transformation.

      The Jacobian was column-scaled: J_scaled[:,j] = J[:,j] * n_j. Solving gives deltaX_scaled, and the true step is deltaX[j] = deltaX_scaled[j] / n_j. Falls back to pseudo-inverse if LU solve fails (e.g., singular matrix).

      Parameters:
      objectiveVector - the right-hand side vector F
      Returns:
      the Newton step deltaX = -J^{-1} F (unscaled), or null if solve fails

    • unscaleNewtonStep

      private void unscaleNewtonStep(double[] deltaX)
      Unscale a Newton step vector using the column scaling factors from the log-mole transformation. The scaled system solves J_s * Δξ = -F where J_s = J * D (D = diag(n_j)). Recovering the original step: Δn = D * Δξ, so Δn_j = n_j * Δξ_j. Lagrange multiplier entries have scale factor 1.0 (no change).
      Parameters:
      deltaX - the Newton step vector to unscale in place

    • enforceMinimumConcentrations

      private void enforceMinimumConcentrations(SystemInterface system)
      Enforce minimum concentration threshold to prevent numerical issues. If any component has moles less than MIN_MOLES, it will be set to MIN_MOLES.
      Parameters:
      system - The thermodynamic system to check and modify

    • getInletMole

      public List<Double> getInletMole()
      Get the inlet mole list.
      Returns:
      List of inlet moles for each component

    • getOutletMole

      public List<Double> getOutletMole()
      Get the outlet mole list.
      Returns:
      List of outlet moles for each component (with 0 values replaced by 1E-6)

    • getInletMoles

      public List<Double> getInletMoles()
      Get the inlet mole list for debugging.
      Returns:
      List of inlet moles

    • getOutletMoles

      public List<Double> getOutletMoles()
      Get the outlet mole list for debugging.
      Returns:
      List of outlet moles

    • printDatabaseComponents

      public void printDatabaseComponents()
      Print loaded database components for debugging.

    • findActiveElements

      private List<Integer> findActiveElements()
      Find which elements are actually present in the system (have non-zero coefficients).
      Returns:
      List of indices of active elements

    • getActiveElementIndices

      private List<Integer> getActiveElementIndices()
      Get the indices of active elements (elements that have non-zero coefficients in any component).
      Returns:
      List of active element indices

    • verifyJacobianInverse

      public boolean verifyJacobianInverse()
      Verify that the Jacobian inverse is correct by multiplying J * J^-1. Should return the identity matrix if the inverse is correct.
      Returns:
      True if the inverse is correct (within tolerance)

    • performNewtonRaphsonIteration

      public double[] performNewtonRaphsonIteration()
      Perform one Newton-Raphson iteration step to calculate the delta vector (dX). Uses LU decomposition to solve J * dX = -F directly (Nocedal & Wright, 2000), falling back to J^{-1} * F if LU solve is not available. The LU approach is ~3x faster and more numerically stable than explicit matrix inversion.
      Returns:
      The delta vector (dX) for updating variables, or null if calculation fails

    • performIterationUpdate

      public boolean performIterationUpdate(double[] deltaX, double alphaComposition)
      Perform a Newton-Raphson iteration update. Updates outlet compositions with damping factor and Lagrange multipliers directly.
      Parameters:
      deltaX - The delta vector from Newton-Raphson iteration
      alphaComposition - Damping factor for composition updates (e.g., 0.0001)
      Returns:
      True if update was successful, false otherwise

    • updateSystemWithNewCompositions

      private boolean updateSystemWithNewCompositions()
      Update the thermodynamic system with the new outlet compositions.
      Returns:
      True if update was successful, false otherwise

    • getFugacityCoefficient

      public double[] getFugacityCoefficient(Object phaseNameOrIndex)
      Get the fugacity coefficient array for all components in a specified phase using the current outlet composition. Uses direct phase composition assignment for efficiency.
      Parameters:
      phaseNameOrIndex - Name or index of the phase (e.g., "gas", "oil", "aqueous", or 0/1/2)
      Returns:
      Fugacity coefficient (phi) array for all components in the specified phase, or Double.NaN if not found

    • setMaxIterations

      public void setMaxIterations(int maxIterations)
      Set maximum number of Newton-Raphson iterations.
      Parameters:
      maxIterations - Maximum number of iterations

    • getMaxIterations

      public int getMaxIterations()
      Get maximum number of Newton-Raphson iterations.
      Returns:
      Maximum number of iterations

    • setConvergenceTolerance

      public void setConvergenceTolerance(double convergenceTolerance)
      Set convergence tolerance for Newton-Raphson iterations.
      Parameters:
      convergenceTolerance - Convergence tolerance

    • getConvergenceTolerance

      public double getConvergenceTolerance()
      Get convergence tolerance for Newton-Raphson iterations.
      Returns:
      Convergence tolerance

    • setDampingComposition

      public void setDampingComposition(double dampingComposition)
      Set damping factor for composition updates in Newton-Raphson iterations.
      Parameters:
      dampingComposition - Damping factor (typically 0.0001 to 0.01)

    • getDampingComposition

      public double getDampingComposition()
      Get damping factor for composition updates.
      Returns:
      Damping factor for composition updates

    • setUseConsistentOffDiagonal

      @Deprecated public void setUseConsistentOffDiagonal(boolean useConsistent)
      Deprecated.
      The consistent off-diagonal formulation is now always enabled. This setter is a no-op.
      Enable the mathematically consistent off-diagonal Jacobian formulation. This is now always enabled (the RT-corrected formulation is the only implementation). This method is retained for backward compatibility but has no effect.
      Parameters:
      useConsistent - ignored — consistent formulation is always active

    • isUseConsistentOffDiagonal

      @Deprecated public boolean isUseConsistentOffDiagonal()
      Deprecated.
      The consistent off-diagonal formulation is now always enabled.
      Check if the consistent off-diagonal Jacobian formulation is enabled. Always returns true.
      Returns:
      always true

    • getActualIterations

      public int getActualIterations()
      Get actual number of iterations performed.
      Returns:
      Actual iterations performed

    • hasConverged

      public boolean hasConverged()
      Check if iterations converged.
      Returns:
      true if converged, false otherwise

    • getFinalConvergenceError

      public double getFinalConvergenceError()
      Get final convergence error.
      Returns:
      Final convergence error

    • setMinIterations

      public void setMinIterations(int minIterations)
      Set minimum number of iterations before convergence is checked. Default is 100. The solver will not declare convergence before completing this many iterations, even if the convergence criterion is satisfied. Setting this too high wastes iterations; too low may cause premature termination.
      Parameters:
      minIterations - Minimum iterations before convergence check (must be at least 1)

    • getMinIterations

      public int getMinIterations()
      Get minimum number of iterations before convergence is checked.
      Returns:
      Minimum iterations before convergence check

    • setUseAdaptiveStepSize

      public void setUseAdaptiveStepSize(boolean useAdaptive)
      Enable or disable NASA CEA-style adaptive step sizing (Gordon & McBride, 1994, NASA RP-1311). When enabled, the step size is automatically computed each iteration to limit the maximum relative change in component moles, allowing larger steps when safe and smaller steps near steep gradients. When disabled, the fixed dampingComposition factor is used for all iterations.
      Parameters:
      useAdaptive - true to enable adaptive step sizing

    • isUseAdaptiveStepSize

      public boolean isUseAdaptiveStepSize()
      Check if adaptive step sizing is enabled.
      Returns:
      true if adaptive step sizing is active

    • setUseArmijoLineSearch

      public void setUseArmijoLineSearch(boolean useArmijo)
      Enable or disable Armijo backtracking line search.

      When enabled, the Newton step is adaptively scaled to guarantee sufficient decrease in the total Gibbs free energy at each iteration. This provides a globally convergent algorithm replacing fixed damping.

      Parameters:
      useArmijo - true to enable Armijo line search

    • isUseArmijoLineSearch

      public boolean isUseArmijoLineSearch()
      Check if Armijo backtracking line search is enabled.
      Returns:
      true if Armijo line search is active

    • setArmijoC1

      public void setArmijoC1(double c1)
      Set the Armijo sufficient decrease parameter c1.
      Parameters:
      c1 - value in (0, 1), typically 1e-4

    • setArmijoRho

      public void setArmijoRho(double rho)
      Set the Armijo backtracking contraction factor.
      Parameters:
      rho - value in (0, 1), typically 0.5

    • setArmijoMaxBacktracks

      public void setArmijoMaxBacktracks(int maxBacktracks)
      Set maximum number of backtracking steps.
      Parameters:
      maxBacktracks - maximum number of backtracks

    • setUseRegularization

      public void setUseRegularization(boolean useReg)
      Enable or disable Tikhonov regularization for ill-conditioned Jacobians.
      Parameters:
      useReg - true to enable regularization

    • isUseRegularization

      public boolean isUseRegularization()
      Check if Tikhonov regularization is enabled.
      Returns:
      true if regularization is active

    • setRegularizationThreshold

      public void setRegularizationThreshold(double threshold)
      Set the condition number threshold above which regularization is applied.
      Parameters:
      threshold - condition number threshold

    • setRegularizationTau

      public void setRegularizationTau(double tau)
      Set the Tikhonov regularization parameter tau.
      Parameters:
      tau - regularization parameter added to Hessian diagonal

    • getConditionNumberHistory

      public List<Double> getConditionNumberHistory()
      Get the condition number history from the last solve (one entry per iteration).
      Returns:
      list of Jacobian condition numbers

    • getGibbsEnergyHistory

      public List<Double> getGibbsEnergyHistory()
      Get the Gibbs energy history from the last solve (one entry per iteration).
      Returns:
      list of total Gibbs energy values (kJ)

    • getStepSizeHistory

      public List<Double> getStepSizeHistory()
      Get the step size history from the last solve (one entry per iteration).
      Returns:
      list of effective step sizes used

    • getElementBalanceErrorHistory

      public List<Double> getElementBalanceErrorHistory()
      Get the element balance error history from the last solve.
      Returns:
      list of element balance error norms

    • evaluateTotalGibbsEnergy

      private double evaluateTotalGibbsEnergy()
      Compute the total Gibbs free energy including mixing terms for the current outlet composition.

      G_total = sum_i n_i * [Gf0_i(T) + RT*ln(phi_i) + RT*ln(y_i) + RT*ln(P/Pref)]

      This method evaluates the actual objective function for line search, distinct from calculateMixtureGibbsEnergy which only uses standard Gibbs energies.

      Returns:
      total Gibbs free energy (kJ) including mixing and non-ideality contributions

    • armijoLineSearch

      private double armijoLineSearch(double[] deltaX, double alphaMax, double currentG)
      Perform Armijo backtracking line search to find a step size that guarantees sufficient decrease in the total Gibbs free energy.

      Starting from alpha = alpha_max, the step is contracted by factor rho until the Armijo condition is satisfied: G(n + alpha*dn) <= G(n) + c1*alpha*grad^T*dn.

      Reference: Nocedal & Wright, Numerical Optimization (2006), Algorithm 3.1.

      Parameters:
      deltaX - the full Newton step vector (composition components only, not Lagrange portion)
      alphaMax - the maximum step size (from adaptive step or fixed damping)
      currentG - current total Gibbs energy
      Returns:
      the accepted step size satisfying the Armijo condition

    • applyRegularization

      private double applyRegularization()
      Apply Tikhonov regularization to the Jacobian matrix if the condition number exceeds the threshold. Adds tau*I to the Hessian (composition-composition) block of the Jacobian.

      This converts the saddle-point system into a positive-definite system when the Hessian is nearly singular, ensuring the Newton direction remains well-defined.

      Returns:
      the condition number of the (possibly regularized) Jacobian

    • calculateAdaptiveAlpha

      private double calculateAdaptiveAlpha(double[] deltaX, double requestedAlpha)
      Calculate adaptive step size using NASA CEA-style step limiting (Gordon & McBride, 1994). Limits the maximum relative mole change so that no major component changes by more than a factor of e^MAX_STEP_LIMIT (~5x) in a single step. Components with moles below a significance threshold are excluded from the limiting (they are growing from near-zero and need large relative steps).
      Parameters:
      deltaX - The raw Newton step vector
      requestedAlpha - The starting step size (typically 1.0)
      Returns:
      The step size bounded to prevent overshooting (clamped to [1e-4, requestedAlpha])

    • solveGibbsEquilibrium

      public boolean solveGibbsEquilibrium(double alphaComposition)
      Solve Gibbs equilibrium using Newton-Raphson iterations with specified step size.
      Parameters:
      alphaComposition - Step size for composition updates
      Returns:
      true if converged, false otherwise

    • solveGibbsEquilibrium

      public boolean solveGibbsEquilibrium()
      Solve Gibbs equilibrium using Newton-Raphson iterations with default damping factor.
      Returns:
      true if converged, false otherwise