Skip to content

ROADMAP.md

Latest commit

 

History

History
267 lines (230 loc) · 18.4 KB

ROADMAP.md

File metadata and controls

267 lines (230 loc) · 18.4 KB
Raw

PhreeqcMatlab — Refactoring & Continuation Roadmap

Status of the codebase at time of writing: development peaked in 2021 and tapered to a last commit in April 2024. Layer 1 (the FFI wrappers) is mature; Layer 2 (orchestration) is partial; Layer 3 (the object model) stalled roughly one-third complete and internally inconsistent. This document is the plan to stabilize, refactor, and continue the work.

Decisions locked in

  • Native libraries: bump to PhreeqcRM/IPhreeqc 3.8.6-17100. Compiled binaries are installed at /usr/local/lib (libphreeqcrm-3.8.6.so, libiphreeqc.so) with headers in /usr/local/include, and source tarballs are in ~/download/. The simulkade/PhreeqcRM release repo will be updated to 3.8.6 later so startup.m can download it.
  • Keep the value-class idiom (phrm = phrm.RM_...(), reassign-or-lose). No migration to handle classes. This must be documented loudly and guarded, since it is the #1 footgun.
  • First focus: stabilize the foundation (Milestone 1) before extending features.

Version-bump compatibility (already verified)

Symbol-level diff of the current wrapper against the installed 3.8.6 .so:

  • All 119 currently-wrapped RM_ functions still exist in 3.8.6 → no symbol-level breakage.
  • 26 new non-BMI functions available to wrap, notably the per-reactant initializers RM_InitialSolutions2Module, RM_InitialEquilibriumPhases2Module, RM_InitialExchanges2Module, RM_InitialSurfaces2Module, RM_InitialGasPhases2Module, RM_InitialSolidSolutions2Module, RM_InitialKinetics2Module, and getters RM_GetTemperature/GetPressure/GetPorosity/GetViscosity/GetDensityCalculated/GetSaturationCalculated.
  • 53 BMI functions (BMIPhreeqcRM) — a standardized Basic Model Interface, available as a modern alternative binding path (Milestone 4, optional).

Signature changes are still possible even where names match; they are only provable by running, which is exactly what the Milestone 1 test harness exists to catch.

Milestone 1 — Stabilize the foundation

Goal: a version-bumped, test-covered, CI-checked, packageable baseline with the known-broken code either fixed or safely guarded. Nothing here adds user-facing features.

1a. Library bump to 3.8.6 — DONE

  • Headers: verified the current curated libs/RM_interface_C.h/IPhreeqc.h are a valid subset of the 3.8.6 headers — all 119 wrapped functions are byte-for-byte signature-compatible (only whitespace differs on 4 density/saturation protos, and their old symbols persist in the 3.8.6 .so). No wrapper changes needed for the bump. New-function prototypes get added to the header in Milestone 4 when those functions are wrapped.
  • Rewrote startup.m: collapsed the four duplicated download blocks into one ensure_library helper; added a local-lib fallback (PHREEQCMATLAB_LIB_PATH → /usr/local/lib) that copies the installed .so into libs/; added a libs/.phreeqc_version stamp so a version bump auto-refreshes; version pinned via *_VERSION constants.
  • Discovered + solved the libstdc++ / GLIBCXX_3.4.32 issue (MATLAB's bundled libstdc++ tops out at 3.4.30; the modern-GCC 3.8.6 build needs 3.4.32). Added run_matlab.sh launcher that sets LD_PRELOAD to the system libstdc++, plus a preflight warning in startup.m.
  • Smoke-tested in MATLAB R2026a: startup → loadlibrary → RM_LoadDatabase → RM_RunString → RM_FindComponents returns correct components end to end.
  • TODO (deferred): checksum verification of downloaded binaries — deferred until the simulkade/PhreeqcRM release is updated to 3.8.6 (the download path isn't exercised yet).

1b. Fix the known-broken code (audit findings) — DONE

  • @SingleCellResult/SingleCellResult.m — renamed classdef SingleCellResults → matching SingleCellResult; rewrote as a minimal loadable placeholder (full fields in M3).
  • @Phase/Phase.m:49 — phase string → phase_string; also corrected the keyword to the plural PHREEQC EQUILIBRIUM_PHASES.
  • @Phase/Phase.m:114-118 — selected_output_object() now writes content(1..4) so all four lines survive.
  • @Phase/Phase.m:274-277 — read_json() now populates phase_names (and matching moles) instead of the nonexistent components.
  • @Gas/Gas.m:41 — wrapped partial_pressure(i) in num2str; guarded the empty-pressure case and put each GAS_PHASE identifier on its own line.
  • PhreeqcSingleCell.m and @SingleCell/SingleCell.m — removed the redundant RM_Create().
  • Broad catch in run()/run_in_phreeqc()/equilibrate_in_phreeqc() (Solution.m, Phase.m, Surface.m) — now surfaces the real ME.message via a PhreeqcMatlab:runFailed warning instead of a generic line. (Full type-stable return-contract redesign deferred to M2.)
  • Removed the template-boilerplate method1/Property1 from @Exchange, @Kinetics, @PhaseResult, @SingleCellResult so the classes load cleanly.
  • Verified in MATLAB R2026a: all touched classes construct; Phase/Gas phreeqc_string, Phase.read_json, and the result classes work; all 4 demo tests (SimpleAdvect, Advect, Species, Gas_m) pass against 3.8.6.
  • Deferred to M3: @SingleCell.run() (undefined iph_string, ignores varargin) — left non-functional for now; it is completed as part of the M3 object-model work.

1c. Test harness — DONE

  • Added tests/PhreeqcMatlabTest.m (matlab.unittest) with physically-grounded golden values: pure-water pH=7 & μ≈1e-7 (IPhreeqc); gypsum SI=0 / anhydrite SI≈−0.30 at 25 °C (IPhreeqc); gypsum solubility Ca=S≈15 mmol/L + component set (PhreeqcSingleCell); and a deterministic Solution.phreeqc_string() well-formedness check. Fixture: tests/fixtures/gypsum_anhydrite.pqc.
  • Single entry point tests/run_all_tests.m (runtests + errors on any failure, for CI). All 4 tests pass against 3.8.6.
  • The demo scripts (SimpleAdvect, Advect, Species, Gas_m) remain as examples; runtests ignores them (not test classes).
  • Findings for M2 (robustness): GetConcentrations segfaults if called before the module is initialized (RM_InitialPhreeqc2Module + RM_RunCells) — the wrapper should guard the precondition rather than pass an unbacked buffer to the native side. Also GetComponents returns a column cell array (orientation matters for callers).

1d. Packaging — DONE

  • package_toolbox.m builds PhreeqcMatlab.mltbx via matlab.addons.toolbox.ToolboxOptions (no hand-written .prj needed). Validated locally: 3.8 MB, ships source + databases + C headers + tests + examples; native binaries are fetched by startup.m, not bundled.
  • CHANGELOG.md added (Keep-a-Changelog style).
  • CI (GitHub Actions) intentionally dropped per maintainer preference; tests are run locally via ./run_matlab.sh -batch "addpath('tests'); run_all_tests".

Milestone 2 — Refactor the object model core (Layer 3)

Goal: make Layer 3 consistent and robust so the remaining classes can be completed cheaply. This is the highest-leverage refactor; do it before finishing individual classes.

  • Abstract Reactant base class — DONE. src/@Reactant/Reactant.m holds the shared name/number identity, an abstract phreeqc_string(), and a concrete input_string() that returns one assembled Phreeqc string. Solution, Phase, Surface, Gas, Exchange, Kinetics now subclass it (name/number removed from each). Surface overrides input_string() to concatenate its three coupled blocks in order. Exchange/Kinetics carry a loud not-implemented phreeqc_string (real bodies land in M3) so they stay instantiable. Covered by reactantPolymorphism (7/7 tests pass).
  • Unify the run/equilibrate verb. Today Solution uses run() while Surface/Phase/Gas use equilibrate_with(). Deferred into M3, where Exchange/Kinetics/Gas/SingleCell are implemented anyway — a shared equilibrate_with(solution) template on Reactant built on input_string() will land with them.
  • Shared string-builder utility — DONE. Added src/Tools/PhreeqcBlock.m, a fluent value builder (kv/kvopt/flag/line) with empty-field suppression, consistent scalar/vector numeric formatting, and deterministic spacing. Refactored Solution, Gas, Phase, Surface (all 3 sub-blocks), and SelectedOutput phreeqc_string() onto it; removed the pe /density 0 malformed lines and the fragile num2str(vector) calls. Every block round-trips through IPhreeqc with no parse error; covered by stringBuilder and surfaceStringRoundTrip tests.
  • Centralize the initial-condition vector — DONE (approach a). Added src/Tools/InitialConditions.m with named slot constants (SOLUTION…KINETICS), a detect(C) input scan, and vectors(present, nxyz) that builds ic1/ic2/f1 (single-cell 1×7 or multi-cell nxyz×7). PhreeqcSingleCell, InitializePhreeqcAdvection, InitializePhreeqcFVTool now share it (was three copies of the keyword scan); Solution.run/Surface.equilibrate_with use the slot constants instead of magic indices. Behavior-preserving, covered by the initialConditionsHelper test. (Approach b — wiring the new per-reactant RM_Initial*2Module functions — is left for M4 when those get wrapped.)
  • Robust result parsing — DONE (Solution; Surface partially). Added src/Tools/map_value.m (safe containers.Map lookup with a fallback). Solution.results_from_phreeqcrm now reads every SELECTED_OUTPUT column via map_value, so a renamed/absent column yields NaN for that field instead of throwing and discarding the whole result. Surface's EDL charge/potential lookups likewise robustified; its positional keys()/values() slicing is explicitly flagged FRAGILE to revisit in M3 with a CD-MUSIC equilibrate reference test (changing it blind is unsafe). While here, fixed two real bugs the now-tested Solution.run path exposed: a missing RM_FindComponents before RunCells (segfault) and a block-concatenation regression (END merged with the next keyword) — combine_phreeqc_strings made newline-robust. Solution.run now returns a populated SolutionResult; covered by solutionRunResults/mapValueSafeLookup. (Swapping scraped columns for the new RM_GetTemperature/... getters is left for M4.)
  • Enum classes: use or lose — DONE (deleted). @solution_units, @phase_units, @exchange_units, @kinetics_units, @sites_units, @edl_layer were unused dead code and wiring them would have broken the strcmpi string comparisons / read_json string assignments. Removed (recoverable via git); src/classes and its addpath are gone. The unit-number conventions they documented remain in inline comments at the RM_SetUnits* calls.
  • Unify JSON — DONE. Added src/Tools/assign_json_fields.m (shared JSON-field→property copier); refactored Solution, Phase, Surface read_json onto it (only the Composition/MasterSpecies/Reactions expansions stay bespoke). Added Solution.to_struct + write_json (round-trips; Composition via containers.Map so element names survive) and a Solution.from_json(name[,file]) factory. Deleted the dead Tools/read_json_ex.m copy and the empty Dan/HDan/Kraka stubs in solutions.json. Covered by jsonRoundTrip.

Milestone 2 complete. 11/11 tests pass.

Milestone 3 — Complete the stubbed classes ✅ COMPLETE

Goal: bring the half-built classes up to the Solution/Surface standard, on the M2 base class.

  • @Reactant equilibration template — added ic_slot() (each reactant's RM_InitialPhreeqc2Module slot) and a shared equilibrate_with(solution) / protected run_with_solution() that centralize the PhreeqcRM boilerplate previously duplicated in Solution.run/Surface.equilibrate_with.

  • @Exchange — real properties (sites/moles + optional master species, exchange reactions, log_k, dh), phreeqc_string(), three-block input_string(), read_json()/from_json(), inherited equilibrate_with(). Added database/exchange.json.

  • @Kinetics — reaction/rate fields (-m0, -m, -parms, -tol, -steps, RATES block), phreeqc_string(), input_string(), read_json()/from_json(), and equilibrate_in_phreeqc() (IPhreeqc, integrates -steps). Added database/kinetics.json with the manual calcite rate.

  • @Gas — implemented equilibrate_in_phreeqc(), inherited equilibrate_with(), fixed the broken selected_output_string(), added read_json()/from_json() and moved damp_CO2()/flue_gas() into database/gases.json.

  • @Phase — equilibrate_with() returns a PhaseResult (moles, moles transferred, SI) plus a SolutionResult, via one combined SELECTED_OUTPUT. (combine_selected_output() remains a % TBD helper — not needed by the object model.)

  • @SingleCell.run() — name-value constructor; builds the combined phreeqc string from all contained reactants, registers initial conditions per slot, RM_RunCells, returns a populated SingleCellResult.

  • Result classes @PhaseResult, @SingleCellResult — real parsed-output structures alongside @SolutionResult/@SurfaceResult.

  • Surface.equilibrate_with refactored (was flagged FRAGILE). Root cause found: it never worked, because combine_surface_solution_string and selected_output_string joined keyword blocks with strjoin's default space delimiter (and a bare concatenation), so PHREEQC misparsed the SURFACE_SPECIES / USER_PUNCH blocks — the positional keys()/values() slicing was moot. Both now join with newlines; the result parsing reads the ordered GetSelectedOutputHeadings/GetSelectedOutput and selects the m_/la_/element column groups by header prefix instead of slicing the alphabetically-sorted containers.Map. A CD-MUSIC calcite surface now equilibrates cleanly with seawater (10 surface species, mole fractions summing to 1, EDL charges/potentials populated); pinned by surfaceEquilibrateCdMusic. Reference models live in examples/phreeqc/chalk_cd_music/ (Wolthers 2008 + Heberling 2011).

Tests: 17/17 pass (6 new — phase/exchange/kinetics/gas equilibration, SingleCell.run, JSON factories).

Milestone 4 — Continue & extend ✅ (BMI deferred)

  • Ship the 3.8.6 C header. The committed libs/RM_interface_C.h was still 3.7.x, so loadlibrary only exposed the old API. Replaced with the 3.8.6 header (214 prototypes) + its irm_dll_export.h (shipped with IRM_DLL_EXPORT empty so MATLAB's thunk compiler can parse it). loadlibrary now binds 193 functions (was ~119).
  • Wrap the new non-BMI 3.8.6 functions in @PhreeqcRM: scalar-field getters (GetTemperature/GetPressure/GetPorosity/GetViscosity/GetDensityCalculated/ GetSaturationCalculated), RM_GetCurrentSelectedOutputUserNumber/RM_SetNthSelectedOutput, RM_Get/SetIthConcentration, RM_Get/SetIthSpeciesConcentration, RM_SetDensityUser/RM_SetSaturationUser, and the seven per-reactant RM_Initial*2Module initializers.
  • .pqm custom-input parser — ParsePqmConfig reads both the 1D (cells/shifts) and multi-D (Nx/Ny/Lx/Ly) forms into a config struct; ApplyRmSettings pushes the PhreeqcRM settings. ReadAdvectionFile now delegates to them (removing its duplicated sscanf ladder and a redundant RM_Create).
  • Multi-D reactive transport. InitializePhreeqcFVTool cleaned up (removed the stale "NOT DONE YET" banner and a stray end); new PhreeqcFVToolTransport couples FVTool advection/diffusion to RM_RunCells by operator splitting. FVTool is auto-provisioned by startup.m (cloned into external/FVTool, gitignored) so it works out of the box; fvtool_available still guards it and the driver errors helpfully if provisioning failed. 2D example: examples/transport/reactive_transport_2d.{m,pqm,pqc}. Verified end-to-end with FVTool: the 2D CaCl₂-flush / cation-exchange run reproduces the expected behaviour (Na displaced out, Ca breaks through, Cl → inflow) — regression test reactiveTransport2D.
  • (Deferred) BMI binding path. The 3.8.6 header exposes 64 RM_Bmi* functions; evaluate wrapping BMIPhreeqcRM as a modern alternative to the RM_ C interface. Left for later — the RM_ interface fully covers current needs.

Tests: 21/21 pass (3 new — .pqm parser 1D/2D, FVTool guard; plus newApi386Getters).

Milestone 5 — Documentation ✅ COMPLETE

Woven throughout, consolidated here:

  • Loudly document the value-class reassign-or-lose idiom and the explicit-destroy requirement — a prominent "two things that bite everyone" callout in the README, plus the class-header docs added across M2–M4 and a dedicated section in docs/architecture.md.
  • API reference + contributor guide: docs/object-model.md (the Layer-3 Reactant contract, every definition class, the result classes and JSON templates), docs/architecture.md (the three-layer map and the RM_ naming rule), and CONTRIBUTING.md (setup, conventions, and how to add a wrapper / definition class / test).
  • README.md refreshed: accurate 3.8.6 install/launch (incl. run_matlab.sh), a high-level object-model quickstart, a "running the tests" section, and an up-to-date roadmap pointer (the stale to-do list is gone).
  • CLAUDE.md kept current (base class, JSON path, 3.8.6 header, .pqm/FVTool helpers).

Remaining (owner task, not code): migrate/expand the GitHub wiki that the README links to — the repo docs above are the source material for it.

Suggested sequencing

  1. M1a + M1b first (bump + fixes) — small, high-confidence, unblocks a trustworthy baseline.
  2. M1c (tests) immediately after — locks the baseline and validates the bump.
  3. M1d (CI/packaging) — automates the gate.
  4. M2 (object-model core refactor) — the pivot that makes M3 cheap.
  5. M3 (complete classes), then M4 (new features / transport), with M5 continuous.