diff --git a/CHANGELOG_AGENT_NOTES.md b/CHANGELOG_AGENT_NOTES.md index 8fcc249065..57afef2046 100644 --- a/CHANGELOG_AGENT_NOTES.md +++ b/CHANGELOG_AGENT_NOTES.md @@ -9,6 +9,87 @@ --- +## 2026-06-24 — Span-aware reference cut boundaries (`ReferenceBoundaryMode.CENTROID_SPAN`) + +### Summary +`characterizeToReference` gains an opt-in boundary-placement mode that fixes +mis-binning when the reference slate mixes a narrow cut with a much wider +neighbour. The default `MIDPOINT` rule places each cut edge at the arithmetic +midpoint of two adjacent cut keys — correct only for equal-width cuts. The new +`CENTROID_SPAN` rule requires each cut key to be the **centroid of its own +span**, walking light→heavy from `b₀ = key₀ − ½(key₁ − key₀)` via the +recurrence `bᵢ = 2·keyᵢ − bᵢ₋₁`, with each boundary clamped strictly between +its adjacent keys (midpoint fallback on overshoot). The recurrence is linear in +molar mass, so it is applied only on the `MOLAR_MASS` basis; the boiling-point +basis falls back to `MIDPOINT`. Additive, no breaking change — the default stays +`MIDPOINT` so all existing results are unchanged. + +### What's new (all additive) +- `CharacterizationOptions`: + - nested enum `ReferenceBoundaryMode {MIDPOINT, CENTROID_SPAN}` (default + `MIDPOINT`) + - getter `getReferenceBoundaryMode()` and builder setter + `referenceBoundaryMode(ReferenceBoundaryMode)` (null-guarded) +- `PseudoComponentCombiner`: `determineReferenceBoundaries(...)` gained a + `ReferenceBoundaryMode` argument and a new private + `determineReferenceCentroidSpanBoundaries(...)` helper; the boundary mode is + threaded through `characterizeToReferenceCore(...)`. + +### Effect +On a Grane→Oseberg per-series case the contested heavy cuts re-align toward the +reference: F4 moves 33.94 → 26.38 (reference 25.21) and F3 5.99 → 1.52 when +`CENTROID_SPAN` is enabled on the molar-mass basis. + +### Verification +```bash +./mvnw test -Dtest=CharacterizeToReferenceCentroidBoundaryTest,CharacterizeToReferenceCommonSlateTest,CharacterizeToReferenceSharedSlateTest,CharacterizeToReferenceDelumpTest,CharacterizationOptionsTest +``` + +### Docs touched +- `docs/pvtsimulation/fluid_characterization_mathematics.md` — `referenceBoundaryMode` + options-table row and "Span-aware reference boundaries" subsection. + +--- + +## 2026-06-24 — Pedersen-faithful gamma delumping for `characterizeToReference` + +### Summary +The delumping stage of `PseudoComponentCombiner` (`delumpContributions`) was +rewritten to follow the Pedersen et al. (Chapter 5) gamma-distribution +treatment instead of the previous fixed ±0.4 carbon-window exponential-decay +heuristic. Each coarse lump is now split into single-carbon-number +sub-fractions over a **neighbour-bounded** carbon range, with sub-fraction +moles drawn from a **Whitson gamma molar distribution on the molar-mass axis**, +molar mass set by Eq. 5.27 (`M = 14·C − 4`), and normal boiling point from the +Katz–Firoozabadi correlation (Eq. 5.28). Parent moles and mass are conserved +according to a configurable mode. Additive, no breaking change. + +### What's new (all additive) +- `CharacterizationOptions` enums + builder methods: + - `DelumpBinningBasis {MOLAR_MASS, BOILING_POINT}` (default `MOLAR_MASS`) + - `DelumpGammaScope {NEIGHBOURS, GLOBAL}` (default `NEIGHBOURS`) + - `DelumpConservation {BOTH, MOLES, MASS}` (default `BOTH`) + - getters `getDelumpBinningBasis()` / `getDelumpGammaScope()` / + `getDelumpConservation()` and matching builder setters. +- `PseudoComponentCombiner`: 4-arg + `delumpContributions(List, int, DelumpGammaScope, DelumpConservation)` gamma + implementation plus molar-mass-basis overloads of + `determineQuantileBoundaries`, `determineReferenceBoundaries`, and + `distributeToProfiles`, and gamma helpers (`fitGlobalGamma`, + `fitGammaByMoments`, `katzFiroozabadiBoilingPoint`, carbon-number/molar-mass + conversions). + +### Verification +```bash +./mvnw test -Dtest=PseudoComponentDelumpDistributionTest,CharacterizeToReferenceDelumpTest,CharacterizeToCommonSlateTest +``` + +### Docs touched +- `docs/pvtsimulation/fluid_characterization_mathematics.md` — gamma delumping + subsection, options table, and builder example. + +--- + ## 2026-06-23 — Process-safety rollout: NOG 070 / STS-0131 / API 14C / NORSOK P-002 / ISO 17776 / EI AVIFF / API 521 §7 / IEC 60079 / API 537 PFP ### Summary diff --git a/docs/pvtsimulation/fluid_characterization_mathematics.md b/docs/pvtsimulation/fluid_characterization_mathematics.md index b96af62738..8e7ace2b28 100644 --- a/docs/pvtsimulation/fluid_characterization_mathematics.md +++ b/docs/pvtsimulation/fluid_characterization_mathematics.md @@ -943,6 +943,9 @@ For advanced control, use the `CharacterizationOptions` builder: ```java import neqsim.thermo.characterization.CharacterizationOptions; import neqsim.thermo.characterization.CharacterizationOptions.NamingScheme; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpBinningBasis; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpGammaScope; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpConservation; CharacterizationOptions options = CharacterizationOptions.builder() .transferBinaryInteractionParameters(true) // Copy BIPs from reference @@ -952,6 +955,9 @@ CharacterizationOptions options = CharacterizationOptions.builder() .inheritReferenceProperties(true) // Inherit reference lump properties .delumpBeforeRecharacterization(false) // Split source lumps before re-binning .delumpResolution(12) // Sub-fractions per source lump + .delumpBinningBasis(DelumpBinningBasis.MOLAR_MASS) // Bin sub-fractions on molar mass + .delumpGammaScope(DelumpGammaScope.NEIGHBOURS) // Local Whitson gamma per lump + .delumpConservation(DelumpConservation.BOTH) // Conserve moles and mass .sharedImaginaryBoundaries(false) // Equal-mass reference cut points .build(); @@ -966,9 +972,13 @@ SystemInterface matched = PseudoComponentCombiner.characterizeToReference( | `namingScheme` | Use SOURCE, REFERENCE, or MERGED names | `REFERENCE` | | `generateValidationReport` | Generate validation report | `false` | | `inheritReferenceProperties` | Inherit the reference lump properties (molar mass, density, critical constants). When `false`, lump properties are recomputed from the source mass on the reference cut grid | `true` | -| `delumpBeforeRecharacterization` | Split each coarse source lump into finer single-carbon-number sub-fractions before re-binning onto the reference cuts (Pedersen Ch. 5 delumping). Conserves source moles and mass exactly | `false` | +| `delumpBeforeRecharacterization` | Split each coarse source lump into finer single-carbon-number sub-fractions before re-binning onto the reference cuts (Pedersen Ch. 5 delumping). Each parent lump is delumped only inside its own neighbour-bounded carbon range using a fitted Whitson gamma molar distribution | `false` | | `delumpResolution` | Number of sub-fractions per source lump when delumping (also the fine grid for equal-mass reference boundaries). Values ≤ 1 disable splitting | `12` | +| `delumpBinningBasis` | Basis on which delumped sub-fractions are binned onto the reference cuts: `MOLAR_MASS` (the conserved, monotonic quantity on which the Pedersen molar distribution is defined) or `BOILING_POINT` (legacy sorting-key behaviour). Applied only when delumping is on | `MOLAR_MASS` | +| `delumpGammaScope` | Scope of the Whitson gamma molar-distribution fit used to shape the delumping: `NEIGHBOURS` (a local gamma per lump fitted from its immediate neighbours) or `GLOBAL` (one gamma fitted to the whole C7+ slate) | `NEIGHBOURS` | +| `delumpConservation` | Quantity conserved exactly per parent lump when delumping: `BOTH` (moles and mass), `MOLES` (moles only, mass floats), or `MASS` (mass only, moles float) | `BOTH` | | `sharedImaginaryBoundaries` | Place reference cut boundaries as carbon-number **equal-mass** cut points on the reference's imaginary (delumped) composition (Pedersen Ch. 5.6, Eqs. 5.58–5.59) instead of boiling-point midpoints | `false` | +| `referenceBoundaryMode` | How each reference cut boundary is placed between adjacent cut keys: `MIDPOINT` (arithmetic midpoint of the two cut keys) or `CENTROID_SPAN` (each cut key is required to be the centroid of its own span). `CENTROID_SPAN` is linear in molar mass and therefore applied only on the `MOLAR_MASS` basis; the boiling-point basis falls back to `MIDPOINT` | `MIDPOINT` | #### Property inheritance vs. recomputation (`inheritReferenceProperties`) @@ -976,7 +986,19 @@ The **Common EoS** slate of Pedersen Ch. 5.6 requires every fluid characterized #### Delumping before re-characterization (`delumpBeforeRecharacterization`) -When a field's native lumps already sit close to the reference grid, the naive source→reference mapping is effectively the identity: lump mole fractions are frozen and per-cut mass is not conserved against the reference molar masses. Enabling `delumpBeforeRecharacterization(true)` first splits each parent lump into `delumpResolution` single-carbon-number sub-fractions whose moles and mass **exactly** reproduce the parent (a single linear molar-mass rescale enforces $\sum_k n_k M_k = n_\text{parent} M_\text{parent}$). The boiling point spreads monotonically with molar mass so sub-fractions can cross reference cut boundaries; density and critical constants stay at the parent values. The sub-fractions are then re-lumped onto the reference grid, so each cut's molar mass = mass / moles is recomputed self-consistently. This is most effective with `inheritReferenceProperties(false)`; combining it with `inheritReferenceProperties(true)` logs a warning because the reference values still overwrite the redistributed lump properties. +When a field's native lumps already sit close to the reference grid, the naive source→reference mapping is effectively the identity: lump mole fractions are frozen and per-cut mass is not conserved against the reference molar masses. Enabling `delumpBeforeRecharacterization(true)` first splits each parent lump into `delumpResolution` single-carbon-number sub-fractions, then re-lumps them onto the reference grid so each cut's molar mass = mass / moles is recomputed self-consistently. The delumping follows the Pedersen et al. (Chapter 5) scheme faithfully: + +1. **Neighbour-bounded carbon range, not a fixed window.** Each parent lump is mapped to a carbon number with the Pedersen molar-mass relation (Eq. 5.27 inverted), $C = (M_{[\text{g/mol}]} + 4)/14$. The sub-fraction interval runs from the carbon-number **midpoint to the lower neighbour** to the **midpoint to the upper neighbour** (first/last lumps mirror their one available gap to the open side). Sub-fractions therefore never cross into a neighbour's range and never fabricate material outside the lump — replacing the earlier arbitrary fixed ±40 % molar-mass window that smeared a lump across carbon numbers it does not contain. + +2. **Whitson gamma molar distribution on the molar-mass axis.** The mole split is shaped by a shifted gamma distribution (Whitson 1983; Pedersen Eq. 5.27) fitted by the method of moments — locally per lump from its immediate neighbours (`delumpGammaScope = NEIGHBOURS`, default) or once over the whole C7+ slate (`GLOBAL`). The gamma is sliced between the single-carbon-number cell molar-mass edges: each cell's mole weight is the gamma probability mass and its molar mass is the gamma conditional mean over the cell (closed form via the regularized lower incomplete gamma function). A fitted gamma is **unimodal** (rises to a mode, then decays) rather than strictly monotonic. When a gamma cannot be fitted (degenerate variance, isolated lump) the method falls back to Pedersen's exponential molar distribution (Eq. 5.15, $\ln z = A + B\,C$) on the cell-centred grid. + +3. **Molar mass and boiling point from non-linear correlations.** Sub-fraction molar masses come from the gamma fit (Eq. 5.27), and normal boiling points from the non-linear Katz–Firoozabadi correlation (Eq. 5.28), $T_b = 97.58\,M^{0.3323}\rho^{0.04609}$ — **not** a $T_b \propto M$ linear spread. When the parent has a boiling point its Katz–Firoozabadi shape is anchored to the parent value, $T_{b,k} = T_{b,\text{parent}}\,\kappa(M_k)/\kappa(M_\text{parent})$, keeping the sub-fractions on the parent/reference boiling-point scale while preserving the non-linear carbon-number dependence. Density and critical constants stay at the parent values. + +4. **User-selectable conservation closure (Eqs. 5.35–5.37).** `delumpConservation` selects what is conserved exactly per parent lump: `BOTH` (default) normalizes the mole weights ($\sum_k n_k = n_\text{parent}$) and rescales the gamma molar masses so $\sum_k n_k M_k = n_\text{parent} M_\text{parent}$ to $10^{-9}$; `MOLES` keeps the gamma molar masses unscaled and lets the mass float; `MASS` rescales the moles so the mass is conserved exactly and lets the moles float. + +The re-binning uses `delumpBinningBasis` (default `MOLAR_MASS`, the conserved monotonic quantity) to place sub-fractions onto the reference cuts. Delumping is most effective with `inheritReferenceProperties(false)`; combining it with `inheritReferenceProperties(true)` logs a warning because the reference values still overwrite the redistributed lump properties. + +The same neighbour-aware delumping also drives the imaginary-composition cut placement used by `sharedImaginaryBoundaries(true)`, so that path benefits from the gamma shape too. #### Equal-mass reference boundaries (`sharedImaginaryBoundaries`) @@ -984,6 +1006,16 @@ By default the reference cut boundaries are the arithmetic boiling-point midpoin > **Why reference-only here?** `characterizeToReference` inherits the trusted reference's lump properties one-to-one, so the reference fluid is the correct authority for the grid. The pooled multi-fluid imaginary composition of Eqs. 5.58–5.59 remains used by `characterizeToCommonSlate` (free new slate); applying it to the inherit-from-fixed-reference path would break the inherit alignment. +#### Span-aware reference boundaries (`referenceBoundaryMode`) + +The default `MIDPOINT` rule places each cut edge at the arithmetic midpoint of the two adjacent reference cut keys. That midpoint of the means is the true span edge only when neighbouring cuts have *equal width*; when a narrow cut sits next to a much wider one, the midpoint is biased toward the wider cut and material that physically belongs to the heavy cut is mis-binned into the narrow neighbour. With `referenceBoundaryMode(ReferenceBoundaryMode.CENTROID_SPAN)` each cut key is instead required to be the **centroid of its own span**: walking light→heavy from the anchor $b_0 = \text{key}_0 - \tfrac{1}{2}(\text{key}_1 - \text{key}_0)$, each successive boundary follows the recurrence + +$$ +b_i = 2\,\text{key}_i - b_{i-1}, +$$ + +so that $\text{key}_i = \tfrac{1}{2}(b_{i-1} + b_i)$ holds exactly. Every boundary is then *clamped* strictly between its two adjacent cut keys (falling back to the midpoint if the recurrence would overshoot a neighbour), preserving the one-to-one inheritance guarantee. The recurrence is linear in molar mass, so `CENTROID_SPAN` is only applied on the `MOLAR_MASS` basis (i.e. with delumping on and `delumpBinningBasis(MOLAR_MASS)`); the boiling-point basis falls back to the midpoint rule. The default stays `MIDPOINT`, so existing results are unchanged unless the mode is set explicitly. + ### Common Slate for Multiple Fluids (`characterizeToCommonSlate`) `characterizeToCommonSlate` implements the Pedersen et al. (Chapter 5.6, Eqs. 5.55–5.60) common-slate procedure for keeping several fluids as **separate** systems while forcing them onto an **identical** pseudo-component slate (so they are mutually EoS-compatible without merging compositions). diff --git a/src/main/java/neqsim/thermo/characterization/CharacterizationOptions.java b/src/main/java/neqsim/thermo/characterization/CharacterizationOptions.java index e4a049ae36..2e00ed6f06 100644 --- a/src/main/java/neqsim/thermo/characterization/CharacterizationOptions.java +++ b/src/main/java/neqsim/thermo/characterization/CharacterizationOptions.java @@ -39,6 +39,66 @@ public enum NamingScheme { CARBON_NUMBER } + /** + * Basis on which delumped single-carbon-number sub-fractions are binned onto the reference cuts. + */ + public enum DelumpBinningBasis { + /** + * Bin on molar mass (the conserved, monotonic quantity; the default). Preferred because molar mass is the basis on + * which the Pedersen molar distribution (Eq. 5.27) and the conservation closures (Eqs. 5.35-5.37) are defined. + */ + MOLAR_MASS, + /** Bin on normal boiling point (the legacy sorting-key behaviour, with molar-mass fallback). */ + BOILING_POINT + } + + /** + * Scope of the Whitson gamma molar-distribution fit used to shape the delumping of a lump into sub-fractions. + */ + public enum DelumpGammaScope { + /** Fit one global gamma to the whole C7+ lump set (smoother, more Whitson-faithful). */ + GLOBAL, + /** Fit a local gamma to each lump from its immediate neighbours (the default). */ + NEIGHBOURS + } + + /** + * Quantity conserved exactly when a coarse lump is delumped into single-carbon-number sub-fractions. + */ + public enum DelumpConservation { + /** + * Conserve both the parent moles and the parent mass (the default). Achieved by normalizing the sub-fraction mole + * weights (moles) and rescaling the gamma conditional-mean molar masses to the parent mass (mass). + */ + BOTH, + /** Conserve the parent moles exactly and keep the sub-fraction molar masses unscaled; the parent mass may drift. */ + MOLES, + /** Conserve the parent mass exactly by rescaling the sub-fraction moles; the parent moles may drift. */ + MASS + } + + /** + * Rule used to place the cut edges between adjacent reference pseudo-components when a source fluid is + * re-characterized onto a fixed reference slate. + */ + public enum ReferenceBoundaryMode { + /** + * Place each cut edge at the arithmetic midpoint of the two adjacent cut keys (the default, legacy behaviour). This + * implicitly assumes the reference cuts are equally wide. + */ + MIDPOINT, + /** + * Require each reference cut key to be the centroid of its own span, giving the recurrence + * {@code b_i = 2*key_i - b_(i-1)} anchored at {@code b_0 = key_0 - 0.5*(key_1 - key_0)} and walked from the + * lightest to the heaviest cut. Each candidate boundary is clamped strictly between its two adjacent keys (with a + * midpoint fallback when it would fall out of range or break monotonicity). The recurrence is linear in the cut + * key, so it is only applied on the {@link DelumpBinningBasis#MOLAR_MASS} basis; the boiling-point basis falls back + * to {@link #MIDPOINT}. Use this when the reference slate has very unequal cut widths, where the midpoint of the + * means mis-bins material into the narrower neighbour. + */ + CENTROID_SPAN + } + private final boolean transferBinaryInteractionParameters; private final boolean normalizeComposition; private final NamingScheme namingScheme; @@ -48,6 +108,10 @@ public enum NamingScheme { private final boolean delumpBeforeRecharacterization; private final int delumpResolution; private final boolean sharedImaginaryBoundaries; + private final DelumpBinningBasis delumpBinningBasis; + private final DelumpGammaScope delumpGammaScope; + private final DelumpConservation delumpConservation; + private final ReferenceBoundaryMode referenceBoundaryMode; private CharacterizationOptions(Builder builder) { this.transferBinaryInteractionParameters = builder.transferBinaryInteractionParameters; @@ -59,6 +123,10 @@ private CharacterizationOptions(Builder builder) { this.delumpBeforeRecharacterization = builder.delumpBeforeRecharacterization; this.delumpResolution = builder.delumpResolution; this.sharedImaginaryBoundaries = builder.sharedImaginaryBoundaries; + this.delumpBinningBasis = builder.delumpBinningBasis; + this.delumpGammaScope = builder.delumpGammaScope; + this.delumpConservation = builder.delumpConservation; + this.referenceBoundaryMode = builder.referenceBoundaryMode; } /** @@ -176,6 +244,61 @@ public boolean isSharedImaginaryBoundaries() { return sharedImaginaryBoundaries; } + /** + * Returns the basis on which delumped single-carbon-number sub-fractions are binned onto the reference cuts. + * + *

+ * Defaults to {@link DelumpBinningBasis#MOLAR_MASS}. Molar mass is the quantity on which the Pedersen molar + * distribution (Eq. 5.27) and the conservation closures are defined, so it gives a monotonic, conservation-faithful + * binning. {@link DelumpBinningBasis#BOILING_POINT} reproduces the legacy sorting-key behaviour. + * + * @return the delumping binning basis (never null) + */ + public DelumpBinningBasis getDelumpBinningBasis() { + return delumpBinningBasis; + } + + /** + * Returns the scope of the Whitson gamma molar-distribution fit used to shape the delumping. + * + *

+ * Defaults to {@link DelumpGammaScope#NEIGHBOURS}, which fits a local gamma to each lump from its immediate + * neighbours. {@link DelumpGammaScope#GLOBAL} fits a single gamma to the whole C7+ lump set. + * + * @return the delumping gamma scope (never null) + */ + public DelumpGammaScope getDelumpGammaScope() { + return delumpGammaScope; + } + + /** + * Returns the quantity conserved exactly when a coarse lump is delumped into single-carbon-number sub-fractions. + * + *

+ * Defaults to {@link DelumpConservation#BOTH}, which conserves both the parent moles and the parent mass. + * {@link DelumpConservation#MOLES} conserves moles only; {@link DelumpConservation#MASS} conserves mass only. + * + * @return the delumping conservation mode (never null) + */ + public DelumpConservation getDelumpConservation() { + return delumpConservation; + } + + /** + * Returns the rule used to place the cut edges between adjacent reference pseudo-components on the reference-only + * re-characterization path. + * + *

+ * Defaults to {@link ReferenceBoundaryMode#MIDPOINT} (the legacy arithmetic-midpoint behaviour). + * {@link ReferenceBoundaryMode#CENTROID_SPAN} requires each reference cut key to be the centroid of its own span and + * is only applied on the {@link DelumpBinningBasis#MOLAR_MASS} basis. + * + * @return the reference boundary placement mode (never null) + */ + public ReferenceBoundaryMode getReferenceBoundaryMode() { + return referenceBoundaryMode; + } + /** * Creates a new builder with default options. * @@ -216,6 +339,10 @@ public static class Builder { private boolean delumpBeforeRecharacterization = false; private int delumpResolution = 12; private boolean sharedImaginaryBoundaries = false; + private DelumpBinningBasis delumpBinningBasis = DelumpBinningBasis.MOLAR_MASS; + private DelumpGammaScope delumpGammaScope = DelumpGammaScope.NEIGHBOURS; + private DelumpConservation delumpConservation = DelumpConservation.BOTH; + private ReferenceBoundaryMode referenceBoundaryMode = ReferenceBoundaryMode.MIDPOINT; /** * Set whether to transfer binary interaction parameters from the reference fluid. @@ -351,6 +478,79 @@ public Builder sharedImaginaryBoundaries(boolean shared) { return this; } + /** + * Set the basis on which delumped single-carbon-number sub-fractions are binned onto the reference cuts. + * + *

+ * Defaults to {@link DelumpBinningBasis#MOLAR_MASS}. Molar mass is the quantity on which the Pedersen molar + * distribution and the conservation closures are defined, giving a monotonic, conservation-faithful binning; + * {@link DelumpBinningBasis#BOILING_POINT} reproduces the legacy boiling-point sorting-key behaviour. + * + * @param basis the delumping binning basis (must not be null) + * @return this builder + */ + public Builder delumpBinningBasis(DelumpBinningBasis basis) { + if (basis != null) { + this.delumpBinningBasis = basis; + } + return this; + } + + /** + * Set the scope of the Whitson gamma molar-distribution fit used to shape the delumping. + * + *

+ * Defaults to {@link DelumpGammaScope#NEIGHBOURS} (a local gamma fitted per lump from its immediate neighbours). + * {@link DelumpGammaScope#GLOBAL} fits a single gamma to the whole C7+ lump set. + * + * @param scope the delumping gamma scope (must not be null) + * @return this builder + */ + public Builder delumpGammaScope(DelumpGammaScope scope) { + if (scope != null) { + this.delumpGammaScope = scope; + } + return this; + } + + /** + * Set the quantity conserved exactly when a coarse lump is delumped into single-carbon-number sub-fractions. + * + *

+ * Defaults to {@link DelumpConservation#BOTH} (parent moles and mass both conserved). + * {@link DelumpConservation#MOLES} conserves moles only; {@link DelumpConservation#MASS} conserves mass only. + * + * @param conservation the delumping conservation mode (must not be null) + * @return this builder + */ + public Builder delumpConservation(DelumpConservation conservation) { + if (conservation != null) { + this.delumpConservation = conservation; + } + return this; + } + + /** + * Set the rule used to place the cut edges between adjacent reference pseudo-components on the reference-only + * re-characterization path. + * + *

+ * Defaults to {@link ReferenceBoundaryMode#MIDPOINT} (legacy arithmetic midpoint of the adjacent cut keys). Set to + * {@link ReferenceBoundaryMode#CENTROID_SPAN} to require each reference cut key to be the centroid of its own span, + * which avoids mis-binning material into a narrower neighbour when the reference cut widths are very unequal. The + * centroid-span rule is only applied on the {@link DelumpBinningBasis#MOLAR_MASS} basis; the boiling-point basis + * falls back to the midpoint rule. + * + * @param mode the reference boundary placement mode (must not be null) + * @return this builder + */ + public Builder referenceBoundaryMode(ReferenceBoundaryMode mode) { + if (mode != null) { + this.referenceBoundaryMode = mode; + } + return this; + } + /** * Build the CharacterizationOptions instance. * diff --git a/src/main/java/neqsim/thermo/characterization/PseudoComponentCombiner.java b/src/main/java/neqsim/thermo/characterization/PseudoComponentCombiner.java index 133410c210..88916c6ab6 100644 --- a/src/main/java/neqsim/thermo/characterization/PseudoComponentCombiner.java +++ b/src/main/java/neqsim/thermo/characterization/PseudoComponentCombiner.java @@ -9,8 +9,13 @@ import java.util.List; import java.util.Map; import java.util.Objects; +import org.apache.commons.math3.special.Gamma; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpBinningBasis; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpConservation; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpGammaScope; +import neqsim.thermo.characterization.CharacterizationOptions.ReferenceBoundaryMode; import neqsim.thermo.component.ComponentInterface; import neqsim.thermo.component.attractiveeosterm.AttractiveTermInterface; import neqsim.thermo.phase.PhaseEos; @@ -162,7 +167,8 @@ public static SystemInterface combineReservoirFluids(int targetPseudoComponents, * @return characterized fluid containing pseudo components compatible with the reference fluid */ public static SystemInterface characterizeToReference(SystemInterface source, SystemInterface reference) { - return characterizeToReferenceCore(source, reference, false, false, 0, false); + return characterizeToReferenceCore(source, reference, false, false, 0, false, DelumpBinningBasis.MOLAR_MASS, + DelumpGammaScope.NEIGHBOURS, DelumpConservation.BOTH, ReferenceBoundaryMode.MIDPOINT); } /** @@ -184,10 +190,19 @@ public static SystemInterface characterizeToReference(SystemInterface source, Sy * less disable splitting * @param sharedImaginaryBoundaries whether the reference cut boundaries are placed as equal-mass cuts on the * reference's imaginary (delumped) composition instead of boiling-point midpoints + * @param binningBasis basis on which delumped sub-fractions are binned onto the reference cuts; only applied when + * {@code delump} is {@code true} (otherwise the legacy boiling-point key is used to preserve backward compatibility) + * @param gammaScope scope of the Whitson gamma molar-distribution fit used to shape the delumping + * @param conservation quantity (moles, mass, or both) conserved exactly when a lump is delumped + * @param boundaryMode rule used to place the cut edges between adjacent reference pseudo-components; only the + * {@link ReferenceBoundaryMode#CENTROID_SPAN} mode on the {@link DelumpBinningBasis#MOLAR_MASS} basis changes the + * legacy midpoint placement * @return characterized fluid containing pseudo components compatible with the reference fluid */ private static SystemInterface characterizeToReferenceCore(SystemInterface source, SystemInterface reference, - boolean inheritReferenceProperties, boolean delump, int delumpResolution, boolean sharedImaginaryBoundaries) { + boolean inheritReferenceProperties, boolean delump, int delumpResolution, boolean sharedImaginaryBoundaries, + DelumpBinningBasis binningBasis, DelumpGammaScope gammaScope, DelumpConservation conservation, + ReferenceBoundaryMode boundaryMode) { Objects.requireNonNull(source, "source"); Objects.requireNonNull(reference, "reference"); @@ -212,16 +227,21 @@ private static SystemInterface characterizeToReferenceCore(SystemInterface sourc return characterized; } + // The molar-mass binning basis is applied only when delumping. With delumping off the legacy boiling-point sorting + // key is retained so the non-delumped and shared-imaginary paths reproduce their previous results exactly. + DelumpBinningBasis basis = delump ? binningBasis : DelumpBinningBasis.BOILING_POINT; + List sourcePseudoComponents = sourceExtraction.pseudoComponents; if (delump) { - sourcePseudoComponents = delumpContributions(sourcePseudoComponents, delumpResolution); + sourcePseudoComponents = delumpContributions(sourcePseudoComponents, delumpResolution, gammaScope, conservation); } List boundaries = sharedImaginaryBoundaries - ? determineReferenceEqualMassBoundaries(referenceExtraction.pseudoComponents, delumpResolution) - : determineReferenceBoundaries(referenceExtraction.pseudoComponents); + ? determineReferenceEqualMassBoundaries(referenceExtraction.pseudoComponents, delumpResolution, basis, + gammaScope, conservation) + : determineReferenceBoundaries(referenceExtraction.pseudoComponents, basis, boundaryMode); List profiles = distributeToProfiles(sourcePseudoComponents, boundaries, - referenceExtraction.pseudoComponents.size()); + referenceExtraction.pseudoComponents.size(), basis); for (int i = 0; i < referenceExtraction.pseudoComponents.size(); i++) { PseudoComponentProfile profile = profiles.get(i); @@ -270,7 +290,8 @@ public static SystemInterface characterizeToReference(SystemInterface source, Sy SystemInterface characterized = characterizeToReferenceCore(source, reference, options.isInheritReferenceProperties(), options.isDelumpBeforeRecharacterization(), - options.getDelumpResolution(), options.isSharedImaginaryBoundaries()); + options.getDelumpResolution(), options.isSharedImaginaryBoundaries(), options.getDelumpBinningBasis(), + options.getDelumpGammaScope(), options.getDelumpConservation(), options.getReferenceBoundaryMode()); if (options.isTransferBinaryInteractionParameters()) { transferBinaryInteractionParameters(reference, characterized); @@ -956,76 +977,271 @@ private static List determineQuantileBoundaries(List - * This implements the delumping stage of the Pedersen et al. (Chapter 5) lumping/delumping scheme (Eqs. 5.35-5.37 - * describe the inverse mass-weighted lumping). For every parent lump the method: + * Identical to {@link #determineQuantileBoundaries(List, int)} but the cut points are placed on the requested + * {@link DelumpBinningBasis} key (molar mass or boiling point) instead of always on the boiling-point sorting key. + * The contributions are sorted by the chosen key before the cumulative-mass scan so the boundaries are returned in + * ascending key order. + * + * @param contributions the contributions to cut + * @param targetPseudoComponents the number of groups (one fewer boundary is returned) + * @param basis the binning key on which equal-mass cut points are placed + * @return the ascending equal-mass boundaries (size {@code targetPseudoComponents - 1}) + */ + private static List determineQuantileBoundaries(List contributions, + int targetPseudoComponents, DelumpBinningBasis basis) { + if (targetPseudoComponents <= 1 || contributions.isEmpty()) { + return Collections.emptyList(); + } + + List sorted = new ArrayList<>(contributions); + sorted.sort(Comparator.comparingDouble(c -> c.binningKey(basis))); + + double totalMass = 0.0; + for (PseudoComponentContribution contribution : sorted) { + totalMass += contribution.mass; + } + + if (!(totalMass > MASS_TOLERANCE)) { + return Collections.emptyList(); + } + + double[] targets = new double[targetPseudoComponents - 1]; + for (int i = 0; i < targets.length; i++) { + targets[i] = totalMass * (i + 1) / targetPseudoComponents; + } + + List boundaries = new ArrayList<>(targets.length); + double cumulative = 0.0; + int targetIndex = 0; + for (PseudoComponentContribution contribution : sorted) { + double nextCumulative = cumulative + contribution.mass; + while (targetIndex < targets.length && nextCumulative >= targets[targetIndex] - MASS_TOLERANCE) { + boundaries.add(contribution.binningKey(basis)); + targetIndex++; + } + cumulative = nextCumulative; + } + + while (boundaries.size() < targets.length) { + boundaries.add(Double.POSITIVE_INFINITY); + } + + return boundaries; + } + + /** Slope of the Pedersen molar-mass / carbon-number relation (Eq. 5.27): M[g/mol] = 14·C - 4. */ + private static final double CARBON_MASS_SLOPE = 14.0; + /** Offset of the Pedersen molar-mass / carbon-number relation (Eq. 5.27): M[g/mol] = 14·C - 4. */ + private static final double CARBON_MASS_OFFSET = 4.0; + /** Leading factor of the Katz-Firoozabadi boiling-point correlation (Eq. 5.28). */ + private static final double KATZ_FACTOR = 97.58; + /** Molar-mass exponent of the Katz-Firoozabadi boiling-point correlation (Eq. 5.28). */ + private static final double KATZ_MASS_EXPONENT = 0.3323; + /** Density exponent of the Katz-Firoozabadi boiling-point correlation (Eq. 5.28). */ + private static final double KATZ_DENSITY_EXPONENT = 0.04609; + /** + * Mild default molar-distribution slope B (Eq. 5.15, {@code ln z = A + B*C}) used only when a lump has no usable + * neighbour from which to estimate the local slope. A negative value gives the customary light-end bias. + */ + private static final double DEFAULT_MOLAR_SLOPE = -0.25; + /** Clamp on the estimated molar-distribution slope B to guard against degenerate neighbour ratios. */ + private static final double MAX_MOLAR_SLOPE = 2.0; + /** Fallback half-width (in carbon numbers) for a lump that has no neighbour on either side. */ + private static final double DEFAULT_CARBON_HALF_WIDTH = 0.5; + /** Fallback liquid density (g/cm3) for the Katz-Firoozabadi correlation when a lump density is unavailable. */ + private static final double DEFAULT_LIQUID_DENSITY = 0.8; + /** Lower clamp on the fitted Whitson gamma shape parameter α (Eq. 5.27 molar distribution). */ + private static final double MIN_GAMMA_SHAPE = 0.5; + /** Upper clamp on the fitted Whitson gamma shape parameter α to guard against a near-delta distribution. */ + private static final double MAX_GAMMA_SHAPE = 50.0; + /** Minimum cell probability mass below which a gamma slice is treated as empty. */ + private static final double GAMMA_WEIGHT_TOLERANCE = 1e-12; + + /** + * Split each coarse source pseudo-component into a grid of finer single-carbon-number (SCN) sub-fractions using a + * neighbour-bounded carbon-number interval and a Whitson gamma molar distribution, conserving the requested quantity + * (moles, mass, or both). + * + *

+ * This implements the delumping stage of the Pedersen et al. (Chapter 5) lumping/delumping scheme. Unlike a fixed + * ±40 % molar-mass window with a hard-coded exponential decay (which smears a lump over carbon numbers it + * does not contain and fabricates spurious light/heavy mass), every parent lump is delumped only inside its own + * carbon-number interval, inferred from its ordered neighbours in the sorted contribution list. For every parent + * lump i with carbon number {@code C_i = (M_i + 4)/14} (Eq. 5.27 inverted) the method: *

    - *
  1. builds a molar-mass grid that brackets the parent molar mass over a symmetric ±{@code window} - * window;
    2. - *
  2. assigns normalized, light-end-biased sub-fraction mole weights (their sum is one, so the parent moles are - * conserved);
    3. - *
  3. rescales the molar-mass grid by a single factor so that Σk nk Mk = - * nparent Mparent exactly, conserving the parent mass;
    4. - *
  4. spreads the normal boiling point monotonically with molar mass so the sub-fractions can cross reference cut - * boundaries (this is what removes the identity source-to-reference mapping), while holding density and the critical - * constants at the parent values.
    5. + *
  5. Bounds the sub-fraction range by the neighbours, not a fixed window. The lower edge is the carbon-number + * midpoint between lumps i-1 and i, the upper edge the midpoint between i and + * i+1. The first and last lumps use a one-sided rule (the available neighbour gap is mirrored to the open + * side). This guarantees sub-fractions never cross into a neighbour's range and never invent material outside the + * lump.
    6. + *
  6. Shapes the molar distribution with a Whitson gamma fitted on molar mass. When + * {@code gammaScope == GLOBAL} a single gamma is fitted by the method of moments to the whole C7+ lump set; when + * {@code gammaScope == NEIGHBOURS} (the default) a local gamma is fitted per lump from its immediate neighbours. The + * gamma is sliced between the SCN cell molar-mass edges: each cell's mole weight is the gamma probability mass and + * its molar mass is the gamma conditional mean over the cell (closed-form via the regularized lower incomplete gamma + * function). When a gamma cannot be fitted (degenerate variance, isolated lump) the method falls back to Pedersen's + * exponential molar distribution (Eq. 5.15, {@code ln z = A + B*C}) on the cell-centred grid. The mole weights are + * normalized, so the parent moles are conserved.
    7. + *
  7. Applies the conservation closure. With {@code DelumpConservation.BOTH} (the default) the sub-fraction + * molar masses are rescaled so that Σk nk Mk = nparent + * Mparent exactly (parent moles and mass both conserved, Eqs. 5.35-5.37). With + * {@code DelumpConservation.MOLES} the gamma molar masses are kept and the mass is allowed to float; with + * {@code DelumpConservation.MASS} the moles are rescaled so the mass is conserved exactly and the moles float.
    8. + *
  8. Assigns the normal boiling point via the non-linear Katz-Firoozabadi correlation (Eq. 5.28), + * {@code Tb = 97.58*M^0.3323*rho^0.04609}. When the parent has a boiling point its Katz-Firoozabadi shape is anchored + * to the parent value ({@code Tb_k = Tb_parent * katz(M_k)/katz(M_parent)}) so the sub-fractions stay on the + * parent/reference boiling-point scale while following the non-linear carbon-number dependence; otherwise the + * absolute Eq. 5.28 value is used. Density and the critical constants are inherited from the parent.
    9. *
* *

* The expanded list is re-sorted by sorting key (normal boiling point, falling back to molar mass) so it can be fed * directly to {@link #distributeToProfiles}. * - * @param sourcePseudoComponents the coarse source pseudo-components to delump + * @param sourcePseudoComponents the coarse source pseudo-components to delump, sorted ascending by sorting key * @param resolution number of sub-fractions per parent lump; values of 1 or less return the input unchanged + * @param gammaScope scope of the Whitson gamma molar-distribution fit (global, or local per lump) + * @param conservation quantity conserved exactly (parent moles, parent mass, or both) * @return the delumped, sorted list of sub-fractions */ private static List delumpContributions( - List sourcePseudoComponents, int resolution) { + List sourcePseudoComponents, int resolution, DelumpGammaScope gammaScope, + DelumpConservation conservation) { if (resolution <= 1) { return sourcePseudoComponents; } - final double window = 0.4; - final double decay = 1.5; + int lumpCount = sourcePseudoComponents.size(); + double[] carbonNumbers = new double[lumpCount]; + for (int i = 0; i < lumpCount; i++) { + carbonNumbers[i] = carbonNumberFromMolarMass(sourcePseudoComponents.get(i).molarMass); + } + + // Optional single global Whitson gamma fitted to the whole C7+ lump set (DelumpGammaScope.GLOBAL). The shift eta is + // the lower molar-mass edge of the lightest lump so every per-lump slice has a non-negative reduced variable. + double[] globalGamma = null; + if (gammaScope == DelumpGammaScope.GLOBAL) { + globalGamma = fitGlobalGamma(sourcePseudoComponents, carbonNumbers); + } - List expanded = new ArrayList<>(sourcePseudoComponents.size() * resolution); + List expanded = new ArrayList<>(lumpCount * resolution); - for (PseudoComponentContribution parent : sourcePseudoComponents) { + for (int i = 0; i < lumpCount; i++) { + PseudoComponentContribution parent = sourcePseudoComponents.get(i); double parentMoles = parent.moles; double parentMolarMass = parent.molarMass; - if (!(parentMoles > 0.0) || !(parentMolarMass > 0.0)) { + double parentCarbon = carbonNumbers[i]; + if (!(parentMoles > 0.0) || !(parentMolarMass > 0.0) || !Double.isFinite(parentCarbon)) { + expanded.add(parent); + continue; + } + + // 1. Neighbour-bounded carbon-number range (sub-fractions never cross into a neighbour's interval). + boolean hasLower = i > 0 && Double.isFinite(carbonNumbers[i - 1]) && carbonNumbers[i - 1] < parentCarbon; + boolean hasUpper = i < lumpCount - 1 && Double.isFinite(carbonNumbers[i + 1]) + && carbonNumbers[i + 1] > parentCarbon; + double upperHalf = hasUpper ? 0.5 * (carbonNumbers[i + 1] - parentCarbon) + : (hasLower ? 0.5 * (parentCarbon - carbonNumbers[i - 1]) : DEFAULT_CARBON_HALF_WIDTH); + double lowerHalf = hasLower ? 0.5 * (parentCarbon - carbonNumbers[i - 1]) + : (hasUpper ? 0.5 * (carbonNumbers[i + 1] - parentCarbon) : DEFAULT_CARBON_HALF_WIDTH); + double lowerCarbon = parentCarbon - lowerHalf; + double upperCarbon = parentCarbon + upperHalf; + if (!(upperCarbon > lowerCarbon)) { expanded.add(parent); continue; } + double width = upperCarbon - lowerCarbon; double[] weights = new double[resolution]; - double weightSum = 0.0; double[] molarMasses = new double[resolution]; - for (int k = 0; k < resolution; k++) { - double t = (double) k / (resolution - 1); - weights[k] = Math.exp(-decay * t); - weightSum += weights[k]; - molarMasses[k] = parentMolarMass * (1.0 - window + 2.0 * window * t); + double weightSum = 0.0; + + // 2. Pick the gamma shape: a single global fit, or a local fit over the {i-1, i, i+1} window. + double[] gamma = gammaScope == DelumpGammaScope.GLOBAL ? globalGamma + : fitLocalGamma(sourcePseudoComponents, carbonNumbers, i, lowerCarbon); + + boolean gammaUsable = false; + if (gamma != null) { + double alpha = gamma[0]; + double beta = gamma[1]; + double eta = gamma[2]; + // 3a. Slice the gamma between the SCN cell edges; cell weight = probability mass, MW = conditional mean. + for (int k = 0; k < resolution; k++) { + double edgeLoCarbon = lowerCarbon + width * k / resolution; + double edgeHiCarbon = lowerCarbon + width * (k + 1) / resolution; + double mwLoG = molarMassFromCarbonNumber(edgeLoCarbon) * 1000.0; + double mwHiG = molarMassFromCarbonNumber(edgeHiCarbon) * 1000.0; + double xLo = Math.max(0.0, (mwLoG - eta) / beta); + double xHi = Math.max(0.0, (mwHiG - eta) / beta); + double prob = Gamma.regularizedGammaP(alpha, xHi) - Gamma.regularizedGammaP(alpha, xLo); + double condMeanG; + if (prob > GAMMA_WEIGHT_TOLERANCE) { + double firstMoment = alpha * beta + * (Gamma.regularizedGammaP(alpha + 1.0, xHi) - Gamma.regularizedGammaP(alpha + 1.0, xLo)); + condMeanG = eta + firstMoment / prob; + } else { + prob = 0.0; + condMeanG = 0.5 * (mwLoG + mwHiG); + } + weights[k] = prob; + molarMasses[k] = condMeanG / 1000.0; + weightSum += prob; + } + gammaUsable = weightSum > GAMMA_WEIGHT_TOLERANCE; + } + + if (!gammaUsable) { + // 3b. Exponential fallback (Eq. 5.15, ln z = A + B*C) on the cell-centred carbon grid. + double slope = estimateMolarSlope(sourcePseudoComponents, carbonNumbers, i, hasLower, hasUpper); + weightSum = 0.0; + for (int k = 0; k < resolution; k++) { + double carbon = lowerCarbon + width * (k + 0.5) / resolution; + double weight = Math.exp(slope * (carbon - parentCarbon)); + weights[k] = weight; + molarMasses[k] = molarMassFromCarbonNumber(carbon); + weightSum += weight; + } + } + + if (!(weightSum > 0.0)) { + expanded.add(parent); + continue; } + // 4. Mole weights (normalized => parent moles conserved) and the requested conservation closure. + double parentMass = parentMoles * parentMolarMass; double rawMass = 0.0; for (int k = 0; k < resolution; k++) { - double moleFraction = weights[k] / weightSum; - rawMass += parentMoles * moleFraction * molarMasses[k]; + rawMass += parentMoles * (weights[k] / weightSum) * molarMasses[k]; + } + double massScale = 1.0; + double moleScale = 1.0; + if (conservation == DelumpConservation.BOTH) { + massScale = rawMass > MASS_TOLERANCE ? parentMass / rawMass : 1.0; + } else if (conservation == DelumpConservation.MASS) { + moleScale = rawMass > MASS_TOLERANCE ? parentMass / rawMass : 1.0; } - double parentMass = parentMoles * parentMolarMass; - double scale = rawMass > MASS_TOLERANCE ? parentMass / rawMass : 1.0; + // 5. Build the sub-fractions with Katz-Firoozabadi (Eq. 5.28) boiling points anchored to the parent. boolean parentHasTb = Double.isFinite(parent.normalBoilingPoint) && parent.normalBoilingPoint > 0.0; + double parentKatz = katzFiroozabadiBoilingPoint(parentMolarMass, parent.density); + boolean anchorTb = parentHasTb && Double.isFinite(parentKatz) && parentKatz > 0.0; for (int k = 0; k < resolution; k++) { double moleFraction = weights[k] / weightSum; - double subMoles = parentMoles * moleFraction; - double subMolarMass = molarMasses[k] * scale; - double ratio = subMolarMass / parentMolarMass; - double subTb = parentHasTb ? parent.normalBoilingPoint * ratio : parent.normalBoilingPoint; + double subMoles = parentMoles * moleFraction * moleScale; + double subMolarMass = molarMasses[k] * massScale; + double subKatz = katzFiroozabadiBoilingPoint(subMolarMass, parent.density); + double subTb; + if (anchorTb && Double.isFinite(subKatz)) { + subTb = parent.normalBoilingPoint * (subKatz / parentKatz); + } else if (Double.isFinite(subKatz) && subKatz > 0.0) { + subTb = subKatz; + } else { + subTb = parent.normalBoilingPoint; + } expanded.add(new PseudoComponentContribution(parent.name + "_d" + k, subMoles, subMolarMass, parent.density, subTb, parent.criticalTemperature, parent.criticalPressure, parent.acentricFactor, parent.criticalVolume, @@ -1039,21 +1255,375 @@ private static List delumpContributions( return expanded; } - private static List determineReferenceBoundaries(List referenceContributions) { + /** + * Package-private test seam that delumps a synthetic set of lumps and returns the resulting single-carbon-number + * sub-fractions as a plain numeric matrix, so unit tests can assert the neighbour-bounded ranges, the gamma molar + * decay, the Katz-Firoozabadi molar-mass / boiling-point relation, and the mole/mass conservation closures without + * needing access to the private {@link PseudoComponentContribution} type. + * + *

+ * Each input lump is described by its moles, molar mass (kg/mol), liquid density (g/cm3), and normal boiling point + * (K); the critical constants and other component properties are inherited from the parent and are immaterial to the + * delumping arithmetic, so they are filled with {@link Double#NaN}. The returned matrix has one row per sub-fraction + * with columns {@code {moles, molarMass[kg/mol], normalBoilingPoint[K], parentIndex}}, sorted ascending by sorting + * key. The {@code parentIndex} is the index of the originating lump in the (sorting-key sorted) input, so a test can + * group the sub-fractions back to their parent without re-deriving the molar-mass interval edges. + * + * @param moles per-lump moles + * @param molarMass per-lump molar mass in kg/mol + * @param density per-lump liquid density in g/cm3 + * @param boilingPoint per-lump normal boiling point in K (use {@link Double#NaN} for none) + * @param resolution number of sub-fractions per lump + * @param gammaScope scope of the Whitson gamma molar-distribution fit + * @param conservation quantity conserved exactly when a lump is delumped + * @return a matrix with one row per sub-fraction and columns {@code {moles, molarMass, normalBoilingPoint, + * parentIndex}} + */ + static double[][] delumpForTesting(double[] moles, double[] molarMass, double[] density, double[] boilingPoint, + int resolution, CharacterizationOptions.DelumpGammaScope gammaScope, + CharacterizationOptions.DelumpConservation conservation) { + List lumps = new ArrayList<>(moles.length); + for (int i = 0; i < moles.length; i++) { + lumps.add(new PseudoComponentContribution("L" + i, moles[i], molarMass[i], density[i], boilingPoint[i], + Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, + Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN)); + } + lumps.sort(Comparator.comparingDouble(PseudoComponentContribution::sortingKey)); + List result = delumpContributions(lumps, resolution, gammaScope, conservation); + double[][] out = new double[result.size()][4]; + for (int i = 0; i < result.size(); i++) { + PseudoComponentContribution sub = result.get(i); + out[i][0] = sub.moles; + out[i][1] = sub.molarMass; + out[i][2] = sub.normalBoilingPoint; + out[i][3] = parseParentIndex(sub.name); + } + return out; + } + + /** + * Parse the originating parent lump index from a {@link #delumpForTesting} sub-fraction name of the form + * {@code L_d}; returns {@code -1} when the name does not match. + * + * @param name the sub-fraction name + * @return the parent index, or {@code -1} when the name is not a delumping sub-fraction + */ + private static double parseParentIndex(String name) { + if (name == null || name.length() < 2 || name.charAt(0) != 'L') { + return -1.0; + } + int suffix = name.indexOf("_d"); + if (suffix <= 1) { + return -1.0; + } + try { + return Double.parseDouble(name.substring(1, suffix)); + } catch (NumberFormatException ex) { + return -1.0; + } + } + + /** + * Package-private test seam exposing the molar-mass / carbon-number conversion (Pedersen Eq. 5.27) so unit tests can + * compute the expected neighbour-bounded carbon-number midpoints. + * + * @param molarMass the molar mass in kg/mol + * @return the equivalent single-carbon number + */ + static double carbonNumberFromMolarMassForTesting(double molarMass) { + return carbonNumberFromMolarMass(molarMass); + } + + /** + * Package-private test seam exposing the carbon-number / molar-mass conversion (Pedersen Eq. 5.27) so unit tests can + * map carbon-number midpoints back onto molar-mass interval edges. + * + * @param carbonNumber the single-carbon number + * @return the equivalent molar mass in kg/mol + */ + static double molarMassFromCarbonNumberForTesting(double carbonNumber) { + return molarMassFromCarbonNumber(carbonNumber); + } + + /** + * Fit a single Whitson gamma molar distribution (Pedersen et al. Eq. 5.27 / Whitson 1983) to the whole set of valid + * C7+ lumps by the method of moments. The shift parameter η is anchored at the lower molar-mass edge of the + * lightest lump's neighbour-bounded interval so that every per-lump gamma slice has a non-negative reduced variable. + * + * @param lumps the parent lumps + * @param carbonNumbers the carbon number of each lump (parallel to {@code lumps}) + * @return {@code {alpha, beta, eta}} in g/mol units, or {@code null} when the distribution is degenerate + */ + private static double[] fitGlobalGamma(List lumps, double[] carbonNumbers) { + int n = lumps.size(); + double[] mwG = new double[n]; + double[] moles = new double[n]; + double minEdgeG = Double.POSITIVE_INFINITY; + int valid = 0; + for (int i = 0; i < n; i++) { + PseudoComponentContribution lump = lumps.get(i); + if (!(lump.moles > 0.0) || !(lump.molarMass > 0.0) || !Double.isFinite(carbonNumbers[i])) { + continue; + } + mwG[valid] = lump.molarMass * 1000.0; + moles[valid] = lump.moles; + double lowerCarbon = carbonNumbers[i] - DEFAULT_CARBON_HALF_WIDTH; + double edgeG = molarMassFromCarbonNumber(lowerCarbon) * 1000.0; + if (edgeG < minEdgeG) { + minEdgeG = edgeG; + } + valid++; + } + if (valid < 2 || !Double.isFinite(minEdgeG)) { + return null; + } + return fitGammaByMoments(mwG, moles, valid, minEdgeG); + } + + /** + * Fit a local Whitson gamma molar distribution to lump {@code i} from its immediate neighbours ({@code i-1, i, i+1}). + * The shift parameter η is anchored at the lump's lower molar-mass edge so the slice reduced variables are + * non-negative. + * + * @param lumps the parent lumps + * @param carbonNumbers the carbon number of each lump (parallel to {@code lumps}) + * @param i the index of the lump being delumped + * @param lowerCarbon the lower carbon-number edge of the lump's neighbour-bounded interval + * @return {@code {alpha, beta, eta}} in g/mol units, or {@code null} when the local distribution is degenerate + */ + private static double[] fitLocalGamma(List lumps, double[] carbonNumbers, int i, + double lowerCarbon) { + double[] mwG = new double[3]; + double[] moles = new double[3]; + int valid = 0; + for (int j = i - 1; j <= i + 1; j++) { + if (j < 0 || j >= lumps.size()) { + continue; + } + PseudoComponentContribution lump = lumps.get(j); + if (!(lump.moles > 0.0) || !(lump.molarMass > 0.0) || !Double.isFinite(carbonNumbers[j])) { + continue; + } + mwG[valid] = lump.molarMass * 1000.0; + moles[valid] = lump.moles; + valid++; + } + if (valid < 2) { + return null; + } + double etaG = molarMassFromCarbonNumber(lowerCarbon) * 1000.0; + return fitGammaByMoments(mwG, moles, valid, etaG); + } + + /** + * Method-of-moments fit of a shifted gamma distribution to a mole-weighted set of molar masses. + * + *

+ * Given the mole-weighted mean μ and variance σ² of the molar masses, the shifted gamma with origin + * η has shape α = (μ-η)²/σ² and scale β = σ²/(μ-η). The shape + * is clamped to {@code [MIN_GAMMA_SHAPE, MAX_GAMMA_SHAPE]} with the scale recomputed from the mean to keep μ = + * η + αβ. + * + * @param mwG the molar masses in g/mol (first {@code count} entries used) + * @param moles the mole weights (first {@code count} entries used) + * @param count the number of valid entries + * @param etaG the gamma origin (shift) in g/mol + * @return {@code {alpha, beta, eta}} in g/mol units, or {@code null} when the fit is degenerate + */ + private static double[] fitGammaByMoments(double[] mwG, double[] moles, int count, double etaG) { + double z = 0.0; + double mean = 0.0; + for (int k = 0; k < count; k++) { + z += moles[k]; + mean += moles[k] * mwG[k]; + } + if (!(z > 0.0)) { + return null; + } + mean /= z; + double variance = 0.0; + for (int k = 0; k < count; k++) { + double d = mwG[k] - mean; + variance += moles[k] * d * d; + } + variance /= z; + double shifted = mean - etaG; + if (!(variance > MASS_TOLERANCE) || !(shifted > 0.0)) { + return null; + } + double alpha = shifted * shifted / variance; + if (alpha < MIN_GAMMA_SHAPE) { + alpha = MIN_GAMMA_SHAPE; + } else if (alpha > MAX_GAMMA_SHAPE) { + alpha = MAX_GAMMA_SHAPE; + } + double beta = shifted / alpha; + if (!(beta > 0.0)) { + return null; + } + return new double[] { alpha, beta, etaG }; + } + + /** + * Convert a molar mass to an equivalent single-carbon number using the Pedersen molar-mass relation (Eq. 5.27) + * inverted: {@code C = (M[g/mol] + 4) / 14}. + * + * @param molarMassKgPerMol molar mass in kg/mol (NeqSim's internal unit) + * @return the equivalent carbon number, or {@link Double#NaN} when the molar mass is non-finite or non-positive + */ + private static double carbonNumberFromMolarMass(double molarMassKgPerMol) { + if (!Double.isFinite(molarMassKgPerMol) || molarMassKgPerMol <= 0.0) { + return Double.NaN; + } + double molarMassGramPerMol = molarMassKgPerMol * 1000.0; + return (molarMassGramPerMol + CARBON_MASS_OFFSET) / CARBON_MASS_SLOPE; + } + + /** + * Convert a carbon number to a molar mass using the Pedersen molar-mass relation (Eq. 5.27): + * {@code M[g/mol] = 14*C - 4}. + * + * @param carbonNumber the (possibly fractional) carbon number + * @return the molar mass in kg/mol (NeqSim's internal unit) + */ + private static double molarMassFromCarbonNumber(double carbonNumber) { + double molarMassGramPerMol = CARBON_MASS_SLOPE * carbonNumber - CARBON_MASS_OFFSET; + return molarMassGramPerMol / 1000.0; + } + + /** + * Evaluate the Katz-Firoozabadi boiling-point correlation (Pedersen et al. Eq. 5.28): + * {@code Tb = 97.58*M^0.3323*rho^0.04609}, with the molar mass in g/mol and the density in g/cm3. + * + * @param molarMassKgPerMol molar mass in kg/mol (NeqSim's internal unit) + * @param density liquid density in g/cm3; a default of 0.8 is used when it is non-finite or non-positive + * @return the correlated normal boiling point in kelvin, or {@link Double#NaN} when the molar mass is invalid + */ + private static double katzFiroozabadiBoilingPoint(double molarMassKgPerMol, double density) { + if (!Double.isFinite(molarMassKgPerMol) || molarMassKgPerMol <= 0.0) { + return Double.NaN; + } + double molarMassGramPerMol = molarMassKgPerMol * 1000.0; + double rho = Double.isFinite(density) && density > 0.0 ? density : DEFAULT_LIQUID_DENSITY; + return KATZ_FACTOR * Math.pow(molarMassGramPerMol, KATZ_MASS_EXPONENT) * Math.pow(rho, KATZ_DENSITY_EXPONENT); + } + + /** + * Estimate the Pedersen molar-distribution slope {@code B} of Eq. 5.15 ({@code ln z = A + B*C}) for lump {@code i} + * from the local lump-to-lump mole ratio over the carbon-number gap to its neighbours. Using both neighbours gives a + * centred estimate; a single neighbour gives a one-sided estimate; with no usable neighbour the mild + * {@link #DEFAULT_MOLAR_SLOPE} is returned. The result is clamped to ±{@link #MAX_MOLAR_SLOPE} to guard + * against degenerate ratios. + * + * @param lumps the sorted parent lumps + * @param carbonNumbers the carbon number of each lump (parallel to {@code lumps}) + * @param i the index of the lump whose slope is estimated + * @param hasLower whether lump {@code i-1} is a usable lighter neighbour + * @param hasUpper whether lump {@code i+1} is a usable heavier neighbour + * @return the estimated, clamped molar-distribution slope {@code B} + */ + private static double estimateMolarSlope(List lumps, double[] carbonNumbers, int i, + boolean hasLower, boolean hasUpper) { + double slope = DEFAULT_MOLAR_SLOPE; + if (hasLower && hasUpper) { + double molesLow = lumps.get(i - 1).moles; + double molesHigh = lumps.get(i + 1).moles; + double deltaCarbon = carbonNumbers[i + 1] - carbonNumbers[i - 1]; + if (molesLow > 0.0 && molesHigh > 0.0 && deltaCarbon > 0.0) { + slope = (Math.log(molesHigh) - Math.log(molesLow)) / deltaCarbon; + } + } else if (hasUpper) { + double molesMid = lumps.get(i).moles; + double molesHigh = lumps.get(i + 1).moles; + double deltaCarbon = carbonNumbers[i + 1] - carbonNumbers[i]; + if (molesMid > 0.0 && molesHigh > 0.0 && deltaCarbon > 0.0) { + slope = (Math.log(molesHigh) - Math.log(molesMid)) / deltaCarbon; + } + } else if (hasLower) { + double molesLow = lumps.get(i - 1).moles; + double molesMid = lumps.get(i).moles; + double deltaCarbon = carbonNumbers[i] - carbonNumbers[i - 1]; + if (molesLow > 0.0 && molesMid > 0.0 && deltaCarbon > 0.0) { + slope = (Math.log(molesMid) - Math.log(molesLow)) / deltaCarbon; + } + } + if (!Double.isFinite(slope)) { + slope = DEFAULT_MOLAR_SLOPE; + } + if (slope > MAX_MOLAR_SLOPE) { + slope = MAX_MOLAR_SLOPE; + } else if (slope < -MAX_MOLAR_SLOPE) { + slope = -MAX_MOLAR_SLOPE; + } + return slope; + } + + private static List determineReferenceBoundaries(List referenceContributions, + DelumpBinningBasis basis, ReferenceBoundaryMode boundaryMode) { if (referenceContributions.size() <= 1) { return Collections.emptyList(); } + if (boundaryMode == ReferenceBoundaryMode.CENTROID_SPAN && basis == DelumpBinningBasis.MOLAR_MASS) { + return determineReferenceCentroidSpanBoundaries(referenceContributions, basis); + } + List boundaries = new ArrayList<>(referenceContributions.size() - 1); for (int i = 0; i < referenceContributions.size() - 1; i++) { - double key1 = referenceContributions.get(i).sortingKey(); - double key2 = referenceContributions.get(i + 1).sortingKey(); + double key1 = referenceContributions.get(i).binningKey(basis); + double key2 = referenceContributions.get(i + 1).binningKey(basis); double boundary = Double.isFinite(key1) && Double.isFinite(key2) ? 0.5 * (key1 + key2) : Math.max(key1, key2); boundaries.add(boundary); } return boundaries; } + /** + * Place the reference cut boundaries so that each reference cut key is the centroid of its own span. + * + *

+ * The default {@link #determineReferenceBoundaries} places each edge at the arithmetic midpoint of two adjacent cut + * keys, which implicitly assumes every cut is equally wide. When the reference slate mixes very narrow and very wide + * cuts (e.g. a one-carbon cut next to a fifteen-carbon cut) the midpoint of the means is not the true span edge, so + * material is mis-binned into the narrower neighbour. Requiring each key to sit at the centroid of its span gives the + * recurrence {@code b_i = 2*key_i - b_(i-1)}, anchored at {@code b_0 = key_0 - 0.5*(key_1 - key_0)} and walked from + * the lightest to the heaviest cut. Each candidate is clamped to lie strictly between its two adjacent keys (and to + * stay monotone), falling back to the arithmetic midpoint otherwise, so the strict one-to-one inheritance ordering is + * preserved. The recurrence is linear in the cut key, so it is only valid on a basis that is itself linear in the + * binning quantity ({@link DelumpBinningBasis#MOLAR_MASS}). + * + * @param referenceContributions the reference pseudo-components, sorted ascending by binning key + * @param basis the (molar-mass) binning key on which the centroid recurrence and clamping bounds are evaluated + * @return the clamped centroid-span cut boundaries (size {@code referenceContributions.size() - 1}) + */ + private static List determineReferenceCentroidSpanBoundaries( + List referenceContributions, DelumpBinningBasis basis) { + int n = referenceContributions.size(); + double[] key = new double[n]; + for (int i = 0; i < n; i++) { + key[i] = referenceContributions.get(i).binningKey(basis); + } + List boundaries = new ArrayList<>(n - 1); + double prev = Double.isFinite(key[0]) && Double.isFinite(key[1]) ? key[0] - 0.5 * (key[1] - key[0]) : key[0]; + for (int i = 0; i < n - 1; i++) { + double midpoint = Double.isFinite(key[i]) && Double.isFinite(key[i + 1]) ? 0.5 * (key[i] + key[i + 1]) + : Math.max(key[i], key[i + 1]); + double candidate = 2.0 * key[i] - prev; + double boundary; + if (Double.isFinite(candidate) && candidate > key[i] && candidate < key[i + 1]) { + boundary = candidate; + } else { + boundary = midpoint; + } + if (!boundaries.isEmpty() && boundary <= boundaries.get(boundaries.size() - 1)) { + boundary = midpoint; + } + boundaries.add(boundary); + prev = boundary; + } + return boundaries; + } + /** * Place the reference cut boundaries as carbon-number-based equal-mass cut points on the reference fluid's imaginary * (fine-resolution) composition. @@ -1077,25 +1647,30 @@ private static List determineReferenceBoundaries(List determineReferenceEqualMassBoundaries( - List referenceContributions, int resolution) { + List referenceContributions, int resolution, DelumpBinningBasis basis, + DelumpGammaScope gammaScope, DelumpConservation conservation) { int cuts = referenceContributions.size() - 1; if (cuts < 1) { return Collections.emptyList(); } if (resolution <= 1) { - return determineReferenceBoundaries(referenceContributions); + return determineReferenceBoundaries(referenceContributions, basis, ReferenceBoundaryMode.MIDPOINT); } - List imaginary = delumpContributions(referenceContributions, resolution); - List equalMass = determineQuantileBoundaries(imaginary, referenceContributions.size()); + List imaginary = delumpContributions(referenceContributions, resolution, gammaScope, + conservation); + List equalMass = determineQuantileBoundaries(imaginary, referenceContributions.size(), basis); List boundaries = new ArrayList<>(cuts); for (int i = 0; i < cuts; i++) { - double lower = referenceContributions.get(i).sortingKey(); - double upper = referenceContributions.get(i + 1).sortingKey(); + double lower = referenceContributions.get(i).binningKey(basis); + double upper = referenceContributions.get(i + 1).binningKey(basis); double midpoint = Double.isFinite(lower) && Double.isFinite(upper) ? 0.5 * (lower + upper) : Math.max(lower, upper); @@ -1149,6 +1724,65 @@ private static List distributeToProfiles(List + * Identical to {@link #distributeToProfiles(List, List, int)} but the contributions are binned on the requested + * {@link DelumpBinningBasis} key. With {@link DelumpBinningBasis#BOILING_POINT} the legacy boiling-point sorting key + * and the original (assumed pre-sorted) behaviour are retained exactly; with {@link DelumpBinningBasis#MOLAR_MASS} a + * copy is first sorted ascending by molar mass so the monotonic boundary scan bins on molar mass. + * + * @param contributions the contributions to distribute + * @param boundaries the ascending cut boundaries on the same key as {@code basis} + * @param targetPseudoComponents the number of groups + * @param basis the binning key + * @return one profile per group (size {@code targetPseudoComponents}) + */ + private static List distributeToProfiles(List contributions, + List boundaries, int targetPseudoComponents, DelumpBinningBasis basis) { + if (targetPseudoComponents <= 0) { + return Collections.emptyList(); + } + + List builders = new ArrayList<>(targetPseudoComponents); + for (int i = 0; i < targetPseudoComponents; i++) { + builders.add(new PseudoComponentGroupBuilder(MASS_TOLERANCE)); + } + + if (contributions.isEmpty()) { + List emptyProfiles = new ArrayList<>(targetPseudoComponents); + for (int i = 0; i < targetPseudoComponents; i++) { + emptyProfiles.add(PseudoComponentProfile.empty()); + } + return emptyProfiles; + } + + List ordered = contributions; + if (basis != DelumpBinningBasis.BOILING_POINT) { + ordered = new ArrayList<>(contributions); + ordered.sort(Comparator.comparingDouble(c -> c.binningKey(basis))); + } + + int groupIndex = 0; + double currentBoundary = boundaries.isEmpty() ? Double.POSITIVE_INFINITY : boundaries.get(0); + + for (PseudoComponentContribution contribution : ordered) { + double key = contribution.binningKey(basis); + while (groupIndex < targetPseudoComponents - 1 && key > currentBoundary) { + groupIndex++; + currentBoundary = groupIndex < boundaries.size() ? boundaries.get(groupIndex) : Double.POSITIVE_INFINITY; + } + builders.get(groupIndex).addContribution(contribution, contribution.mass); + } + + List profiles = new ArrayList<>(targetPseudoComponents); + for (PseudoComponentGroupBuilder builder : builders) { + profiles.add(builder.buildWithOverrides(null, null, null, null)); + } + return profiles; + } + private static List combineProfiles(List> perFluidProfiles, double[] fluidMassTotals, double[] fluidMoleTotals) { if (perFluidProfiles.isEmpty()) { @@ -1318,6 +1952,21 @@ private double sortingKey() { return key; } + /** + * Return the key on which this contribution is binned for the requested basis. The molar-mass basis returns the + * molar mass directly (the conserved, monotonic quantity); the boiling-point basis returns the legacy + * {@link #sortingKey()} (normal boiling point, falling back to molar mass). + * + * @param basis the binning basis + * @return the binning key value + */ + private double binningKey(DelumpBinningBasis basis) { + if (basis == DelumpBinningBasis.MOLAR_MASS) { + return molarMass; + } + return sortingKey(); + } + /** * Returns a copy of this contribution with its mole amount (and hence mass) scaled by the given factor. All * intensive properties and the {@link #sortingKey()} are preserved, so the copy occupies the same position on the diff --git a/src/test/java/neqsim/thermo/characterization/CharacterizeToReferenceCentroidBoundaryTest.java b/src/test/java/neqsim/thermo/characterization/CharacterizeToReferenceCentroidBoundaryTest.java new file mode 100644 index 0000000000..c473e30d66 --- /dev/null +++ b/src/test/java/neqsim/thermo/characterization/CharacterizeToReferenceCentroidBoundaryTest.java @@ -0,0 +1,141 @@ +package neqsim.thermo.characterization; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertNotNull; +import static org.junit.jupiter.api.Assertions.assertNull; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import neqsim.thermo.characterization.CharacterizationOptions.ReferenceBoundaryMode; +import neqsim.thermo.component.ComponentInterface; +import neqsim.thermo.system.SystemInterface; +import neqsim.thermo.system.SystemPrEos; + +/** + * Tests the optional span-aware ({@link ReferenceBoundaryMode#CENTROID_SPAN}) reference cut-point placement of + * {@link PseudoComponentCombiner#characterizeToReference(SystemInterface, SystemInterface, CharacterizationOptions)} + * selected through {@link CharacterizationOptions#getReferenceBoundaryMode()}. + * + *

+ * The default {@link ReferenceBoundaryMode#MIDPOINT} places each cut edge at the arithmetic midpoint of the two + * adjacent cut keys, which implicitly assumes equal cut widths. When the reference slate mixes a narrow cut with a much + * wider one, the midpoint of the means is not the true span edge and material is mis-binned into the narrower + * neighbour. The {@link ReferenceBoundaryMode#CENTROID_SPAN} mode instead requires each cut key to be the centroid of + * its own span, which moves the contested edge and re-assigns the affected material. The span recurrence is linear in + * molar mass, so it is only applied on the {@link CharacterizationOptions.DelumpBinningBasis#MOLAR_MASS} basis (reached + * here by enabling delumping with resolution one, which leaves the source pseudo-components unspread). + */ +class CharacterizeToReferenceCentroidBoundaryTest { + + /** + * Reference slate with three cuts of very unequal width on the molar-mass axis: a narrow F1/F2 pair (100, 120 g/mol) + * next to a far heavier F3 (300 g/mol). Under MIDPOINT the F2/F3 edge lands at 210 g/mol; under CENTROID_SPAN the F2 + * span centred on 120 g/mol forces the edge down to 130 g/mol. + */ + private static SystemInterface contestedReference() { + SystemInterface fluid = new SystemPrEos(298.15, 60.0); + fluid.addComponent("methane", 0.6); + fluid.addTBPfraction("F1", 0.10, 0.100, 0.74); + fluid.addTBPfraction("F2", 0.10, 0.120, 0.76); + fluid.addTBPfraction("F3", 0.10, 0.300, 0.85); + return fluid; + } + + /** A single source pseudo-component at 150 g/mol, in the contested 130-210 g/mol band. */ + private static SystemInterface contestedSource() { + SystemInterface fluid = new SystemPrEos(298.15, 60.0); + fluid.addComponent("methane", 0.8); + fluid.addTBPfraction("S1", 0.20, 0.150, 0.80); + return fluid; + } + + /** Options that reach the MOLAR_MASS basis (delump, resolution one) with the requested boundary mode. */ + private static CharacterizationOptions options(ReferenceBoundaryMode mode) { + return CharacterizationOptions.builder().delumpBeforeRecharacterization(true).delumpResolution(1) + .inheritReferenceProperties(false).referenceBoundaryMode(mode).build(); + } + + /** Name of the single populated pseudo-component cut, or {@code null} if none received any mass. */ + private static String binnedCutName(SystemInterface fluid) { + String name = null; + for (int i = 0; i < fluid.getNumberOfComponents(); i++) { + ComponentInterface c = fluid.getComponent(i); + if ((c.isIsTBPfraction() || c.isIsPlusFraction()) && c.getNumberOfmoles() > 0.0) { + name = c.getName(); + } + } + return name; + } + + /** Total pseudo-fraction moles in a fluid. */ + private static double pseudoMoles(SystemInterface fluid) { + double moles = 0.0; + for (int i = 0; i < fluid.getNumberOfComponents(); i++) { + ComponentInterface c = fluid.getComponent(i); + if (c.isIsTBPfraction() || c.isIsPlusFraction()) { + moles += c.getNumberOfmoles(); + } + } + return moles; + } + + @Test + @DisplayName("referenceBoundaryMode defaults to MIDPOINT and reproduces the midpoint binning") + void testDefaultIsMidpoint() { + SystemInterface source = contestedSource(); + SystemInterface reference = contestedReference(); + + CharacterizationOptions defaults = CharacterizationOptions.builder().delumpBeforeRecharacterization(true) + .delumpResolution(1).inheritReferenceProperties(false).build(); + assertEquals(ReferenceBoundaryMode.MIDPOINT, defaults.getReferenceBoundaryMode(), + "referenceBoundaryMode should default to MIDPOINT"); + + SystemInterface defaultResult = PseudoComponentCombiner.characterizeToReference(source, reference, defaults); + SystemInterface midpointResult = PseudoComponentCombiner.characterizeToReference(source, reference, + options(ReferenceBoundaryMode.MIDPOINT)); + + assertEquals(binnedCutName(midpointResult), binnedCutName(defaultResult), + "the default mode must bin the source into the same cut as explicit MIDPOINT"); + assertEquals(pseudoMoles(midpointResult), pseudoMoles(defaultResult), 1e-9, + "the default mode must conserve the same pseudo moles as explicit MIDPOINT"); + } + + @Test + @DisplayName("CENTROID_SPAN moves the unequal-width edge and re-assigns the contested lump") + void testCentroidSpanReassignsContestedLump() { + SystemInterface source = contestedSource(); + SystemInterface reference = contestedReference(); + double sourceMoles = pseudoMoles(source); + + SystemInterface midpointResult = PseudoComponentCombiner.characterizeToReference(source, reference, + options(ReferenceBoundaryMode.MIDPOINT)); + SystemInterface centroidResult = PseudoComponentCombiner.characterizeToReference(source, reference, + options(ReferenceBoundaryMode.CENTROID_SPAN)); + + // MIDPOINT: F2/F3 edge at 210 g/mol -> the 150 g/mol source bins to F2. + assertEquals("F2_PC", binnedCutName(midpointResult), "MIDPOINT must bin the 150 g/mol lump into F2"); + assertNull(componentOrNull(midpointResult, "F3_PC"), "MIDPOINT must leave F3 empty"); + + // CENTROID_SPAN: F2/F3 edge at 130 g/mol -> the 150 g/mol source bins to F3. + assertEquals("F3_PC", binnedCutName(centroidResult), "CENTROID_SPAN must re-assign the 150 g/mol lump to F3"); + assertNull(componentOrNull(centroidResult, "F2_PC"), "CENTROID_SPAN must leave F2 empty"); + + // Re-assignment must conserve the source pseudo moles in both modes. + assertEquals(sourceMoles, pseudoMoles(midpointResult), sourceMoles * 1e-6, "MIDPOINT conserves pseudo moles"); + assertEquals(sourceMoles, pseudoMoles(centroidResult), sourceMoles * 1e-6, "CENTROID_SPAN conserves pseudo moles"); + + ComponentInterface f3 = componentOrNull(centroidResult, "F3_PC"); + assertNotNull(f3, "CENTROID_SPAN must populate F3"); + assertEquals(0.20, f3.getNumberOfmoles(), 1e-6, "the full 0.20 mol of the source lump must land in F3"); + } + + /** Returns the named component if present and carrying mass, otherwise {@code null}. */ + private static ComponentInterface componentOrNull(SystemInterface fluid, String name) { + for (int i = 0; i < fluid.getNumberOfComponents(); i++) { + ComponentInterface c = fluid.getComponent(i); + if (c.getName().equals(name) && c.getNumberOfmoles() > 0.0) { + return c; + } + } + return null; + } +} diff --git a/src/test/java/neqsim/thermo/characterization/PseudoComponentDelumpDistributionTest.java b/src/test/java/neqsim/thermo/characterization/PseudoComponentDelumpDistributionTest.java new file mode 100644 index 0000000000..a7f676bd12 --- /dev/null +++ b/src/test/java/neqsim/thermo/characterization/PseudoComponentDelumpDistributionTest.java @@ -0,0 +1,200 @@ +package neqsim.thermo.characterization; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpConservation; +import neqsim.thermo.characterization.CharacterizationOptions.DelumpGammaScope; + +/** + * Direct unit tests of the neighbour-aware, gamma-shaped delumping stage of {@link PseudoComponentCombiner}, exercised + * through the package-private {@link PseudoComponentCombiner#delumpForTesting} seam on a synthetic three-lump fluid. + * + *

+ * The four invariants required of a Pedersen-faithful (Chapter 5) delumping are checked independently: + *

    + *
  1. Neighbour-bounded ranges — the middle lump's sub-fraction molar masses stay inside the carbon-number + * midpoints to its lower and upper neighbours; no fabricated mass crosses into a neighbour's interval.
    2. + *
  2. Monotonic molar decay — with a light-biased gamma the sub-fraction mole weights decrease with + * increasing molar mass (Eq. 5.15 / Whitson gamma behaviour).
    3. + *
  3. Eq. 5.27 / 5.28 MW–Tb relation — the assigned normal boiling points follow the non-linear + * Katz-Firoozabadi molar-mass dependence ({@code Tb_i/Tb_j = (M_i/M_j)^0.3323} for equal-density sub-fractions of one + * parent), i.e. boiling point is not linear in molar mass.
    4. + *
  4. Mole and mass conservation — per parent lump the sub-fraction moles sum to the parent moles and the + * sub-fraction mass sums to the parent mass, to 1e-9, when {@link DelumpConservation#BOTH} is requested.
    5. + *
+ */ +class PseudoComponentDelumpDistributionTest { + + /** Katz-Firoozabadi molar-mass exponent (Eq. 5.28). */ + private static final double KATZ_MASS_EXPONENT = 0.3323; + + // Synthetic three-lump fluid: an evenly carbon-spaced light/medium/heavy slate (C ~ 7.4, 11.7, 17.4). + private static final double[] MOLES = { 0.50, 0.30, 0.20 }; + private static final double[] MOLAR_MASS = { 0.100, 0.160, 0.240 }; // kg/mol + private static final double[] DENSITY = { 0.73, 0.79, 0.83 }; // g/cm3 + private static final double[] BOILING_POINT = { 360.0, 430.0, 500.0 }; // K + private static final int RESOLUTION = 6; + + private static double[][] delump(DelumpGammaScope scope, DelumpConservation conservation) { + return PseudoComponentCombiner.delumpForTesting(MOLES, MOLAR_MASS, DENSITY, BOILING_POINT, RESOLUTION, scope, + conservation); + } + + /** + * Group the flat sub-fraction matrix back into the three parent lumps using the parent-index column returned by the + * test seam. Returns one {@code double[resolution][4]} block per parent in ascending parent-index order, preserving + * the ascending sorting-key order within each block. + * + * @param subFractions the flat sub-fraction matrix from {@link PseudoComponentCombiner#delumpForTesting} + * @return the sub-fractions grouped per parent lump + */ + private static double[][][] groupByParent(double[][] subFractions) { + double[][][] groups = new double[3][RESOLUTION][4]; + int[] counts = new int[3]; + for (double[] sub : subFractions) { + int parent = (int) Math.round(sub[3]); + assertTrue(parent >= 0 && parent < 3, "unexpected parent index " + parent); + groups[parent][counts[parent]++] = sub; + } + for (int p = 0; p < 3; p++) { + assertEquals(RESOLUTION, counts[p], "each parent lump must produce exactly " + RESOLUTION + " sub-fractions"); + } + return groups; + } + + @Test + @DisplayName("middle lump delumps only within its neighbour-bounded carbon range") + void testNeighbourBoundedRanges() { + // Use MOLES conservation so the gamma conditional-mean molar masses are kept unscaled and therefore stay strictly + // inside the neighbour-bounded carbon cells. (The BOTH closure rescales molar masses to conserve mass, which can + // nudge the top sub-fraction marginally past the upper midpoint; range-bounding is a property of the cell grid.) + double[][] subs = delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.MOLES); + + double cLow = PseudoComponentCombiner.carbonNumberFromMolarMassForTesting(MOLAR_MASS[0]); + double cMid = PseudoComponentCombiner.carbonNumberFromMolarMassForTesting(MOLAR_MASS[1]); + double cHigh = PseudoComponentCombiner.carbonNumberFromMolarMassForTesting(MOLAR_MASS[2]); + double lowerEdge = PseudoComponentCombiner.molarMassFromCarbonNumberForTesting(0.5 * (cLow + cMid)); + double upperEdge = PseudoComponentCombiner.molarMassFromCarbonNumberForTesting(0.5 * (cMid + cHigh)); + + double[][][] groups = groupByParent(subs); + double[][] middle = groups[1]; + for (double[] sub : middle) { + double mw = sub[1]; + assertTrue(mw >= lowerEdge - 1e-9, + "middle sub-fraction MW " + mw + " kg/mol fell below the lower-neighbour midpoint " + lowerEdge); + assertTrue(mw <= upperEdge + 1e-9, + "middle sub-fraction MW " + mw + " kg/mol exceeded the upper-neighbour midpoint " + upperEdge); + } + } + + @Test + @DisplayName("sub-fraction mole weights form a single-peaked (unimodal) gamma distribution") + void testUnimodalMolarDistribution() { + double[][] subs = delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.BOTH); + double[][][] groups = groupByParent(subs); + + for (int p = 0; p < 3; p++) { + double[][] lump = groups[p]; + // sub-fractions arrive sorted ascending by sorting key (boiling point), hence ascending molar mass. + for (int k = 1; k < lump.length; k++) { + assertTrue(lump[k][1] > lump[k - 1][1] - 1e-12, "sub-fraction molar masses must be ascending within a lump"); + } + // A fitted Whitson gamma (shape >= 0.5) is unimodal: the mole weights rise to a single peak then fall, with no + // second hump. Allow a tiny tolerance against floating-point wiggle around the peak. + int directionChanges = 0; + for (int k = 1; k < lump.length - 1; k++) { + boolean risingBefore = lump[k][0] - lump[k - 1][0] > 1e-12; + boolean fallingAfter = lump[k + 1][0] - lump[k][0] < -1e-12; + boolean fallingBefore = lump[k][0] - lump[k - 1][0] < -1e-12; + boolean risingAfter = lump[k + 1][0] - lump[k][0] > 1e-12; + if ((risingBefore && fallingAfter) || (fallingBefore && risingAfter)) { + directionChanges++; + } + } + assertTrue(directionChanges <= 1, + "parent " + p + ": mole-weight profile must be unimodal but had " + directionChanges + " turning points"); + } + } + + @Test + @DisplayName("boiling point follows the non-linear Katz-Firoozabadi MW dependence (Eq. 5.27/5.28)") + void testKatzFiroozabadiMwTbRelation() { + double[][] subs = delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.BOTH); + double[][][] groups = groupByParent(subs); + + for (int p = 0; p < 3; p++) { + double[][] lump = groups[p]; + double[] ref = lump[0]; + for (int k = 1; k < lump.length; k++) { + double mwRatio = lump[k][1] / ref[1]; + double tbRatio = lump[k][2] / ref[2]; + double expectedTbRatio = Math.pow(mwRatio, KATZ_MASS_EXPONENT); + assertEquals(expectedTbRatio, tbRatio, 1e-6, + "parent " + p + " slice " + k + ": Tb ratio must follow (M_i/M_j)^0.3323"); + // Confirm the relation is genuinely non-linear: a linear Tb-proportional-to-MW model would give + // tbRatio=mwRatio. + assertTrue(tbRatio < mwRatio - 1e-9 || Math.abs(mwRatio - 1.0) < 1e-9, + "Tb must grow slower than MW (non-linear), but tbRatio " + tbRatio + " >= mwRatio " + mwRatio); + } + } + } + + @Test + @DisplayName("BOTH conserves parent moles and mass per lump to 1e-9") + void testMoleAndMassConservationBoth() { + double[][] subs = delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.BOTH); + double[][][] groups = groupByParent(subs); + + for (int p = 0; p < 3; p++) { + double sumMoles = 0.0; + double sumMass = 0.0; + for (double[] sub : groups[p]) { + sumMoles += sub[0]; + sumMass += sub[0] * sub[1]; + } + assertEquals(MOLES[p], sumMoles, 1e-9, "parent " + p + ": moles not conserved"); + assertEquals(MOLES[p] * MOLAR_MASS[p], sumMass, 1e-9, "parent " + p + ": mass not conserved"); + } + } + + @Test + @DisplayName("MOLES conserves parent moles exactly; MASS conserves parent mass exactly") + void testSingleQuantityConservationModes() { + double[][][] molesGroups = groupByParent(delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.MOLES)); + for (int p = 0; p < 3; p++) { + double sumMoles = 0.0; + for (double[] sub : molesGroups[p]) { + sumMoles += sub[0]; + } + assertEquals(MOLES[p], sumMoles, 1e-9, "MOLES mode: parent " + p + " moles not conserved"); + } + + double[][][] massGroups = groupByParent(delump(DelumpGammaScope.NEIGHBOURS, DelumpConservation.MASS)); + for (int p = 0; p < 3; p++) { + double sumMass = 0.0; + for (double[] sub : massGroups[p]) { + sumMass += sub[0] * sub[1]; + } + assertEquals(MOLES[p] * MOLAR_MASS[p], sumMass, 1e-9, "MASS mode: parent " + p + " mass not conserved"); + } + } + + @Test + @DisplayName("global gamma scope also conserves moles and mass per lump") + void testGlobalGammaScopeConserves() { + double[][] subs = delump(DelumpGammaScope.GLOBAL, DelumpConservation.BOTH); + double[][][] groups = groupByParent(subs); + for (int p = 0; p < 3; p++) { + double sumMoles = 0.0; + double sumMass = 0.0; + for (double[] sub : groups[p]) { + sumMoles += sub[0]; + sumMass += sub[0] * sub[1]; + } + assertEquals(MOLES[p], sumMoles, 1e-9, "GLOBAL scope parent " + p + ": moles not conserved"); + assertEquals(MOLES[p] * MOLAR_MASS[p], sumMass, 1e-9, "GLOBAL scope parent " + p + ": mass not conserved"); + } + } +}