Skip to content

sharedImaginaryBoundaries places the boiling-point midpoints for most cuts, because the equal-mass quantile is clamped away #4307

Closed

Description

@PatrickMockridge

CharacterizationOptions.isSharedImaginaryBoundaries() documents its effect as placing the reference cut boundaries as equal-mass cut points on the reference fluid's imaginary (delumped) composition "instead of arithmetic boiling-point midpoints", and PseudoComponentCombiner.determineReferenceEqualMassBoundaries implements exactly that. In practice it returns the midpoints for most cuts.

The symptom, through the public API

Using the shared-slate test's own fluids (CharacterizeToReferenceSharedSlateTest.referenceFluid(6) and .gradedSource()):

SystemInterface reference = referenceFluid(6);
SystemInterface source = gradedSource();

CharacterizationOptions off = CharacterizationOptions.builder()
    .inheritReferenceProperties(false).sharedImaginaryBoundaries(false).build();
CharacterizationOptions on = CharacterizationOptions.builder()
    .inheritReferenceProperties(false).sharedImaginaryBoundaries(true).build();

SystemInterface a = PseudoComponentCombiner.characterizeToReference(source, reference, off);
SystemInterface b = PseudoComponentCombiner.characterizeToReference(source, reference, on);

a and b are identical — same pseudo-components, same moles, same molar mass, same normal-liquid density:

shared=false C8_PC    moles=0.200000000000 M=0.090000000 rho=0.710000000
shared=false C9_PC    moles=0.150000000000 M=0.140000000 rho=0.770000000
shared=false C11_PC   moles=0.100000000000 M=0.210000000 rho=0.820000000
shared=true  C8_PC    moles=0.200000000000 M=0.090000000 rho=0.710000000
shared=true  C9_PC    moles=0.150000000000 M=0.140000000 rho=0.770000000
shared=true  C11_PC   moles=0.100000000000 M=0.210000000 rho=0.820000000

The cause

Both boundary methods are private static, so comparing them directly needs reflection: build the contribution list the way extractComponents does — one PseudoComponentContribution per TBP/plus component with positive moles, sorted by sortingKey — and call both.

cut 0: midpoint 360.579   equal-mass 359.860   moved
cut 1: midpoint 403.897   equal-mass 403.897   the midpoint
cut 2: midpoint 451.702   equal-mass 451.702   the midpoint
cut 3: midpoint 501.090   equal-mass 501.090   the midpoint
cut 4: midpoint 550.830   equal-mass 550.830   the midpoint
-> 4 of 5 cuts are the midpoint

Across six reference slates at the default delumpResolution of 12, 3 of 23 cuts differ — and three of those six slates have every cut identical.

Why it is structural rather than a corner case

delumpContributions spreads each reference lump's normal boiling point over ± window (window = 0.4):

molarMasses[k] = parentMolarMass * (1.0 - window + 2.0 * window * t);
...
subTb = parentHasTb ? parent.normalBoilingPoint * ratio : parent.normalBoilingPoint;

So a lump's sub-fractions carry keys below and above the lump's own key, and the imaginary composition interleaves across neighbouring lumps instead of staying inside its own. determineQuantileBoundaries then answers with the key of the sub-fraction whose running mass crosses the target — and that key lands inside a lump's own spread far more often than in the narrow gap between two lumps. The clamp rejects any candidate not strictly inside the gap:

if (Double.isFinite(lower) && Double.isFinite(upper) && upper > lower) {
  double candidate = i < equalMass.size() ? equalMass.get(i) : Double.NaN;
  if (Double.isFinite(candidate) && candidate > lower && candidate < upper) {
    boundary = candidate;
  }
}

and boundary stays at the midpoint otherwise — the rule the option exists to replace.

The clamp's stated purpose is "preserving the strict one-to-one ordering on which the property inheritance in characterizeToReferenceCore relies". The midpoint is already strictly inside the same gap, so the fallback does not buy that property; it only removes the equal-mass placement.

No test can tell

CharacterizeToReferenceSharedSlateTest has five tests: the default is off, flag-off is backwards compatible, the total moles and mass are conserved, the one-to-one ordering survives for an identical source, and resolution <= 1 falls back to midpoints. None asserts that the equal-mass cuts differ from the midpoint cuts, so the suite passes unchanged with the flag doing nothing.

What would resolve it

If the equal-mass placement is meant to bite, the fallback should keep it inside the gap rather than discard it — clamp the quantile into (lower, upper) instead of replacing it, or take the nearest quantile that falls in the gap. If the midpoint fallback is the intent, the accessor's documentation and the class Javadoc should say so, since both currently describe a cut rule the code does not apply.

Complete runnable reproduction
import java.lang.reflect.Constructor;
import java.lang.reflect.Field;
import java.lang.reflect.Method;
import java.util.ArrayList;
import java.util.List;
import neqsim.thermo.characterization.CharacterizationOptions;
import neqsim.thermo.characterization.PseudoComponentCombiner;
import neqsim.thermo.component.ComponentInterface;
import neqsim.thermo.system.SystemInterface;
import neqsim.thermo.system.SystemPrEos;

public class SharedBoundariesRepro {

  private static final String COMBINER = "neqsim.thermo.characterization.PseudoComponentCombiner";
  private static final String CONTRIBUTION = COMBINER + "$PseudoComponentContribution";

  /** `CharacterizeToReferenceSharedSlateTest.referenceFluid(6)`, verbatim. */
  static SystemInterface referenceFluid(int cuts) {
    SystemInterface fluid = new SystemPrEos(298.15, 60.0);
    fluid.addComponent("methane", 0.6);
    double[] molarMass = {0.080, 0.100, 0.130, 0.160, 0.200, 0.240};
    double[] density = {0.70, 0.73, 0.76, 0.79, 0.81, 0.83};
    for (int i = 0; i < cuts; i++) {
      fluid.addTBPfraction("C" + (7 + i), 0.10, molarMass[i], density[i]);
    }
    return fluid;
  }

  /** `CharacterizeToReferenceSharedSlateTest.gradedSource()`, verbatim. */
  static SystemInterface gradedSource() {
    SystemInterface fluid = new SystemPrEos(298.15, 60.0);
    fluid.addComponent("methane", 0.5);
    fluid.addTBPfraction("S1", 0.20, 0.090, 0.71);
    fluid.addTBPfraction("S2", 0.15, 0.140, 0.77);
    fluid.addTBPfraction("S3", 0.10, 0.210, 0.82);
    return fluid;
  }

  static void print(String label, SystemInterface fluid) {
    for (int i = 0; i < fluid.getNumberOfComponents(); i++) {
      ComponentInterface c = fluid.getComponent(i);
      if ((c.isIsTBPfraction() || c.isIsPlusFraction()) && c.getNumberOfmoles() > 0.0) {
        System.out.printf("%s %-8s moles=%.12f M=%.9f rho=%.9f%n", label, c.getComponentName(),
            c.getNumberOfmoles(), c.getMolarMass(), c.getNormalLiquidDensity());
      }
    }
  }

  /** The two boundary rules, on one contribution list. Both methods are private static. */
  static void boundaries(SystemInterface reference) throws Exception {
    Class<?> row = Class.forName(CONTRIBUTION);
    Constructor<?> ctor = row.getDeclaredConstructor(String.class,
        Class.forName("neqsim.thermo.component.ComponentInterface"));
    ctor.setAccessible(true);
    Field tb = row.getDeclaredField("normalBoilingPoint");
    Field mass = row.getDeclaredField("molarMass");
    tb.setAccessible(true);
    mass.setAccessible(true);

    List<Object> rows = new ArrayList<>();
    for (int i = 0; i < reference.getNumberOfComponents(); i++) {
      ComponentInterface c = reference.getComponent(i);
      if (c.isIsTBPfraction() || c.isIsPlusFraction()) {
        rows.add(ctor.newInstance(c.getComponentName(), c));
      }
    }
    // `extractComponents` sorts by sorting key before either method is handed the list.
    for (int i = 1; i < rows.size(); i++) {
      for (int j = i; j > 0; j--) {
        Object a = rows.get(j - 1);
        Object b = rows.get(j);
        double ka = ((Number) tb.get(a)).doubleValue();
        double kb = ((Number) tb.get(b)).doubleValue();
        double va = Double.isFinite(ka) && ka > 0.0 ? ka : ((Number) mass.get(a)).doubleValue();
        double vb = Double.isFinite(kb) && kb > 0.0 ? kb : ((Number) mass.get(b)).doubleValue();
        if (va <= vb) {
          break;
        }
        rows.set(j - 1, b);
        rows.set(j, a);
      }
    }

    Method withResolution = Class.forName(COMBINER)
        .getDeclaredMethod("determineReferenceEqualMassBoundaries", List.class, int.class);
    Method midpoints =
        Class.forName(COMBINER).getDeclaredMethod("determineReferenceBoundaries", List.class);
    withResolution.setAccessible(true);
    midpoints.setAccessible(true);
    @SuppressWarnings("unchecked")
    List<Double> equal = (List<Double>) withResolution.invoke(null, rows, 12);
    @SuppressWarnings("unchecked")
    List<Double> mid = (List<Double>) midpoints.invoke(null, rows);

    int same = 0;
    for (int i = 0; i < mid.size(); i++) {
      boolean identical = Math.abs(equal.get(i) - mid.get(i)) < 1e-9;
      same += identical ? 1 : 0;
      System.out.printf("cut %d: midpoint %8.3f  equal-mass %8.3f  %s%n", i, mid.get(i), equal.get(i),
          identical ? "the midpoint" : "moved");
    }
    System.out.printf("-> %d of %d cuts are the midpoint%n", same, mid.size());
  }

  public static void main(String[] args) throws Exception {
    SystemInterface reference = referenceFluid(6);
    SystemInterface source = gradedSource();

    CharacterizationOptions off = CharacterizationOptions.builder()
        .inheritReferenceProperties(false).sharedImaginaryBoundaries(false).build();
    CharacterizationOptions on = CharacterizationOptions.builder()
        .inheritReferenceProperties(false).sharedImaginaryBoundaries(true).build();
    print("shared=false", PseudoComponentCombiner.characterizeToReference(source, reference, off));
    print("shared=true ", PseudoComponentCombiner.characterizeToReference(source, reference, on));

    boundaries(referenceFluid(6));
  }
}
javac -proc:none -cp neqsim-f0c7436.jar SharedBoundariesRepro.java
java -cp .:neqsim-f0c7436.jar SharedBoundariesRepro

Reproduced against f0c7436c6923766b1e22957b7075f650600457a7 (neqsim-f0c7436.jar).

Activity

  1. EvenSol commented on Oct 14, 2026

    @EvenSol
    Collaborator

    Reviewed against current master b0de5165abe8582dbfef1ea1b78fc87c11e6e043. The root cause is confirmed: finite quantiles outside their adjacent-reference gap are rejected, so the midpoint survives despite the documented clamping intent.

    Implemented the fix in #4330 (commit 6226b1fed968aab1caf4e61d359e4024fd2cde0f). Candidates are projected onto [Math.nextUp(lower), Math.nextDown(upper)]; already-valid interior candidates stay unchanged. Degenerate gaps, unavailable/non-finite candidates and resolution <= 1 retain the existing fallback.

    For the reported six-cut reference, all five corrected cuts now differ from the midpoints:

    Cut Midpoint (K) Corrected cut (K, rounded)
    0 360.579000 359.860359
    1 403.896500 429.023000
    2 451.701500 474.380000
    3 501.090000 527.800000
    4 550.830000 573.860000

    Cuts 1–4 are actually the floating-point values immediately below the upper reference keys, so they remain strictly inside the gaps.

    One nuance from the review: the original three coarse source lumps still stay in the same bins after this correction. None of their sorting keys crosses a moved boundary. That output alone cannot distinguish correct cut placement from a no-op. New public-API source probes demonstrate both upper-edge and lower-edge projection (and fail on the old code); the delumped version of the reported source also shows a composition change.

    Validation: 40 related tests passed, including 11 shared-boundary tests covering observable redistribution, interior candidates, exact boundary inclusion, reference identity, conservation, backwards compatibility and degenerate gaps. Spotless, both pre-commit stages, documentation search, engineering coverage, Checkstyle's enforced gate and focused Javadoc passed. Full commands and warning-level limitations are recorded in the PR.

    Documentation now explicitly describes constrained mass quantiles, rather than promising equal final cut masses after projection. The default path and flash calculations are unchanged. The PR is draft pending hosted CI and independent review; this issue remains open until the fix is merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions