v3.6.0

• Correctness release. A full literature-verified review of every calculation module found and fixed a set of long-standing correlation bugs. Several outputs change materially - notably Motiee hydrate temperatures, Gray and Beggs-Brill VLP pressures, linear gas rates, stock-tank GOR and undersaturated oil viscosity. Each fix below was confirmed against the original paper (or two independent secondary sources) before the change was made. All frozen regression baselines and documentation examples were re-pinned to the corrected values.

• Gas Z-factor fixes:

  • Hall-Yarborough now actually converges. The Newton loop exited after exactly one step (a loop-variable bug present since the original implementation), so HY returned "initial guess + one step", in error by up to 14% near the critical region (Tpr 1.05, Ppr 2). The solver now iterates to a 5e-4 relative step with a warning on non-convergence. HY Z values shift accordingly (mid-range shifts are small, e.g. 0.86774 to 0.86775 at 2000 psia / sg 0.75 / 200 degF).
  • DAK Newton derivative corrected (the R5 term derivative used rhor**3 where the algebra gives rhor**2 - (A11*rhor^2)^2 terms). Converged Z values shift by < 1e-6; convergence is restored to quadratic (about 5 iterations instead of 18 near-critical). A non-convergence warning was added.
  • WYW Z-factor method removed from the public API. The Wang, Ye & Wu (2021) explicit fit behaves incorrectly at low Ppr (Z tending to 0.944 as Ppr tends to 0) and existed mainly as the HY initial guess, which it remains internally. zmethod='WYW' now raises ValueError.
  • Single-sided tc/pc overrides honoured on Rust paths. Effective tc/pc are now resolved once via gas_tc_pc and forwarded to the Rust kernels, so supplying only one override gives identical answers with and without the accelerator (previously the Rust path silently ignored a lone override). The DAK and HY array paths also now route through the Rust batch kernels for all cmethods, roughly 2x faster than the previous per-element dispatch and bit-identical.

• gas_rate_linear over-predicted by exactly 2 pi. The linear-flow branch reused the radial constant 1422, which embeds the 2 pi of radial geometry. A separate linear-flow constant (2 pi x 1422) is now used; linear gas rates drop by a factor of ~6.28 (doc example: 8.20 to 1.31 Mscf/d). Radial rates are unchanged.

• Motiee hydrate temperature corrected to the published psia/degF basis. The correlation was evaluated in kPa/degC with a transcribed intercept, reading roughly 30-36 degF hot (1000 psia, sg 0.65: 96.2 degF, vs Towler 61.1 and the Katz chart 62-66). The verified Motiee (1991) form now gives 60.2 degF at the same point. Hydrate formation pressures and inhibitor dose estimates from hydmethod='MOTIEE' shift substantially; the default TOWLER method was always correct.

• Oil correlation fixes:

  • oil_rs_st returned ln(Rst). Valko-McCain (2003) Eq 3-2 defines ln(Rst) as the polynomial; the missing exponential is now applied (doc example: 4.18 to 65.13 scf/stb, consistent with the 0.1618 x Rsp heuristic). Confirmed against the original paper, including the psia separator-pressure basis (nomenclature, p168). Metric output is now also converted correctly.
  • Petrosky-Farshad undersaturated viscosity used natural log where McCain Eq 3.24b specifies log10, producing unphysical viscosity growth above Pb (6x over 6000 psi in the worst verified case). PVTO undersaturated viscosities and viscosity-dependent doc examples shift.
  • Vogel IPR edge cases. The Pb clamp only applied to array pwf in oil_rate_radial and never in oil_rate_linear, so scalar undersaturated calls under-predicted by ~22%; saturated reservoirs (pr <= pb) could return negative rates at pwf = pr. Both functions now share an elementwise helper: composite Darcy+Vogel above Pb, pure Vogel driven by pr when saturated (rate is zero at pwf = pr).
  • oil_rs_bub VALMC inversion is now analytic (quadratic then cubic root selection) with machine-precision round-trip, replacing a secant loop that diverged silently above the correlation's representable Pb maximum (returning e.g. rsb of 476,000). Requests beyond the maximum now raise ValueError. The docstring's claimed STAN fallback (which never existed) is removed.
  • oil_rs with rsmethod='STAN' now honours user rsb/pb (scaled like VALMC), removing an Rs discontinuity at Pb. oil_bo with bomethod='STAN' now applies the McCain cofb undersaturated correction above Pb instead of staying constant.
  • Oil density above Pb now evaluates the bubble-point anchor with rsb (matching the Rust path and McCain's formulation), keeping density continuous at Pb; oil_co uses a one-sided difference near Pb so compressibility stays on a single smooth branch. oil_deno clamps the temperature term so inputs below 60 degF return real values. Minor breaking change: the unused pi parameter was removed from oil_co.
  • import pyrestoolbox.oil is now lazy (pandas/tabulate/gas/brine load only when needed), about 4x faster.

• Brine fixes:

  • Salinity-method routing in SoreideWhitson. The documented 'sechenov' and 'auto' options silently applied no salting-out at all (freshwater answers at any salinity); they now normalise to 'gamma_phi'. framework='proposed' + 'embedded' claimed a gamma_phi fallback that never happened - it now happens. framework='dropin' + 'embedded' now wires the delta-kij tables into the flash. framework='sw_original' with the default 'gamma_phi' double-counted salting-out (CH4 solubility ~3x understated at 20 wt%); the Sechenov gamma is now skipped for gases whose kij already embeds salinity. The engine raises ValueError on unrecognised salinity methods. Regression tests pin x_CH4 at 20 wt% NaCl for every framework/method combination.
  • Spycher-Pruess 99-109 degC blend corrected (verified against the 2004 and 2010 papers): the blending weights were inverted (returning the high-temperature model at 99 degC and vice versa), and the Python fugacity blend was a no-op through list aliasing, which was also a real Python/Rust divergence inside the band. Both fixed in both languages; xCO2 is continuous across the band boundaries.
  • Saturated brine compressibility derivative (Spivey Eq 4.33 chain) corrected: the Sechenov pressure-derivative was lumped into the denominator instead of subtracted, a 3.5-6% error in saturated cw for saline brines.
  • Rs mole-basis corrections: CO2_Brine_Mixture converted Spycher's ionised-basis xCO2 with paired-NaCl brine moles (Rs ~ -6.7% at 20 wt%); SoreideWhitson converted the salt-free flash x with brine moles (Rs ~ +7.7% at 20 wt%). Both reduce to the previous behaviour at zero salinity.
  • Constant consistency: 273 to 273.15 K in the brine_props chain, the 0.1015 MPa typos to 0.1013, standard-conditions CH4 molar volume at 288.706 K / 0.101325 MPa, and unified MW and Tb constants that had drifted between modules. Freshwater viscosity shifts ~0.14%; the two brine models now agree on identical inputs.
  • No more NaN cascade below the water vapour pressure: brine_props at low pressure / high temperature now returns finite values with zero dissolved methane. Pressures above the IAPWS-IF97 Region 1 limit (100 MPa = 14,503 psia) raise a clear error naming the limit.
  • Removed ~850 lines of dead code from the VLE engine and salting library (unrouted kij variants, superseded methods, demo blocks).

• Nodal VLP fixes (each applied identically in Python and the Rust accelerator):

  • Gray (1974) holdup rewritten to the published API 14B form. The previous implementation built its velocity group without surface tension and applied the R/(R+1) transform twice, over-predicting holdup ~30x in mist flow and putting Gray BHPs far out of family (condensate-well example: 1940 to 1067 psia, now beside HB/WG/BB at 1077-1209). The effective (wet film) roughness is also now the published piecewise form.
  • Beggs-Brill liquid velocity number used sigma in lbf/ft where the 1.938 constant requires dyne/cm (~11x inflation corrupting the inclination correction; BHP previously fell as water rate rose). The flow-pattern map now matches the published revised map (two regions were misclassified as intermittent instead of distributed), and the Payne et al (1979) 0.924 correction now applies to all uphill patterns including distributed. BB BHPs shift ~1% in typical cases (doc example: 1213 to 1224 psia).
  • Hagedorn-Brown two-phase Reynolds constant corrected from 96778 to 1900.8 (published 2.2e-2 with mass rate in lbm/day). Mostly masked in fully rough turbulence; HB BHPs shift on smooth pipe (doc example: 954.8 to 962.1 psia).
  • Divergence guard: an infeasible injection case previously returned a large negative BHP silently (e.g. -239,000 psia). The march now raises RuntimeError when pressure falls below atmospheric, identically from the Python and Rust paths.
  • The low-rate static-column fallback now warns when liquid rates are being ignored; fbhp/outflow_curve/operating_point take gsg from a supplied gas_pvt; fthp's thp_min metric sentinel is fixed (default is now None); negative rates raise ValueError.
  • Internal: the Rust VLP port now mirrors the Python _segment_march scaffold (one shared march, single calc_segments/condensate_dropout, ~75 named constants with citations), removing ~1,200 duplicated lines - the structure that bred the 3.4.0/3.5.0 constant-drift bugs. All eight entry points reproduce pure-Python to better than 1e-9.

• lorenz_2_layers built its saturation grid with np.arange, which appended a spurious extra layer for 28 of nlayers in 2-200 (e.g. 7 layers averaging k=86 when asked for 6 averaging 100). Now np.linspace; a regression sweep pins nlayers 2-50.

• oil_matbal metric Gi was understated by exactly 5.6146x (sm3 multiplying an internal rb/scf Bg), inflating OOIP ~34% on gas-injection cases in metric units. Field-unit results were always correct; a field-vs-metric equivalence test now pins this.

• simtools fixes: rel_perm_table now honours sgcr for SGOF and returns exactly the documented row count in all configurations (two off-by-one cases fixed); make_bot_og in metric mode now actually converts the Bg column it relabels to rm3/sm3 (~178x), no longer overwrites the harmonised sg_sp, and reports the harmonised pb/rsb in its results dict; zip_check_sim_deck no-files path raises FileNotFoundError instead of TypeError; zip_check_sim_deck and ix_extract_problem_cells default to non_interactive=True (behaviour change - pass False for the old prompts). Internal: simtools.py is split into domain sub-files (_decks, _relperm, _aquifer, _vfp, _pvt_tables) with the public API unchanged.

• DCA: forecast() cumulative/EUR now uses the analytic cumulative functions (uptime-scaled) instead of a right-rectangle rate sum that biased EUR ~ -1.7% low at dt=1; fit_decline(method='best') compares R2 on the common t > 0 subset so the comparison is apples-to-apples when t contains 0.

• Package infrastructure: pyrestoolbox.__version__ now exists (single-sourced from pyproject.toml) and the Rust extension reports its own version via get_status(); exceptions raised inside Rust-accelerated functions now propagate instead of silently falling back to Python; validate_pe_inputs rejects NaN and gains a scalar fast path ~40x faster; passing an Enum of the wrong class raises a descriptive ValueError instead of KeyError; bisect_solve raises when the root is not bracketed; a new GitHub Actions workflow runs the full suite (pure-Python and Rust-accelerated, including parity tests) on every push and pull request; run_all_tests.py now delegates to pytest and runs everything.

• 870 validation tests (up from 838 in 3.5.0), including 113 Python/Rust parity tests.