Skip to content

Replaced PhreeqcRM interface ctypes with phreeqcrm pip package #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
matteomasi wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Replaced PhreeqcRM interface ctypes with phreeqcrm pip package #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ba5362b
Added phreeqcrm to dependencies and ignored SwigPy deprecation warnin…
matteomasi Apr 10, 2026
152efac
Deleted PhreeqcRM compiled libraries
matteomasi Apr 10, 2026
df6c446
IRMStatus module no more necessary
matteomasi Apr 10, 2026
9982c01
Updated PhreecRM class and tests
matteomasi Apr 10, 2026
24a3b62
Updated examples with latest PhreeqcRM API changes
matteomasi Apr 10, 2026
f7d5aa7
Updated notebook examples
matteomasi Apr 10, 2026
94491ac
Updated documentation notebooks
matteomasi Apr 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 8 additions & 10 deletions docs/notebooks/ex1-titration.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -162,10 +162,8 @@
"\n",
"\n",
"# Initialize concentration vectors\n",
"cc = np.zeros(n_cells * n_comps, dtype=np.float64)\n",
"cs = np.zeros(n_cells * n_species, dtype=np.float64)\n",
"phr.rm_get_concentrations(cc)\n",
"phr.rm_get_species_concentrations(cs)\n",
"cc = phr.rm.GetConcentrations()\n",
"cs = phr.rm.GetSpeciesConcentrations()\n",
"\n",
"# HCl concentration vector\n",
"hcl = np.linspace(hcl_range[0], hcl_range[1], n_cells) # mol/L\n",
Expand All @@ -175,7 +173,7 @@
"indx_h = np.where(species == \"H+\")[0][0]\n",
"\n",
"# Add HCl\n",
"cs_r = cs.reshape(n_species, n_cells).T\n",
"cs_r = np.array(cs).reshape(n_species, n_cells).T\n",
"cs_r[:, indx_cl] += hcl # Cl-\n",
"cs_r[:, indx_h] += hcl # H+\n",
"cs1 = cs_r.T.reshape(n_cells * n_species) # Updated concentrations (after adding HCl)"
Expand Down Expand Up @@ -205,13 +203,13 @@
],
"source": [
"# Tranfer the modified species concentrations back to PhreeqcRM\n",
"phr.rm_species_concentrations2_module(cs1)\n",
"phr.rm_set_time(1.0)\n",
"phr.rm_set_time_step(1.0)\n",
"phr.rm.SpeciesConcentrations2Module(cs1)\n",
"phr.rm.SetTime(1.0)\n",
"phr.rm.SetTimeStep(1.0)\n",
"\n",
"# Run the simulation\n",
"start_time = time.time()\n",
"phr.rm_run_cells()\n",
"phr.rm.RunCells()\n",
"elapsed = time.time() - start_time\n",
"print(f\"Simulation completed in {elapsed:.2f} seconds\")"
]
Expand Down
28 changes: 13 additions & 15 deletions docs/notebooks/ex2-BTEX_dissolution.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

23 changes: 11 additions & 12 deletions docs/notebooks/ex3-BTEX_dissolution_transport_coupling.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

18 changes: 8 additions & 10 deletions examples/ex1-titration.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -162,10 +162,8 @@
"\n",
"\n",
"# Initialize concentration vectors\n",
"cc = np.zeros(n_cells * n_comps, dtype=np.float64)\n",
"cs = np.zeros(n_cells * n_species, dtype=np.float64)\n",
"phr.rm_get_concentrations(cc)\n",
"phr.rm_get_species_concentrations(cs)\n",
"cc = phr.rm.GetConcentrations()\n",
"cs = phr.rm.GetSpeciesConcentrations()\n",
"\n",
"# HCl concentration vector\n",
"hcl = np.linspace(hcl_range[0], hcl_range[1], n_cells) # mol/L\n",
Expand All @@ -175,7 +173,7 @@
"indx_h = np.where(species == \"H+\")[0][0]\n",
"\n",
"# Add HCl\n",
"cs_r = cs.reshape(n_species, n_cells).T\n",
"cs_r = np.array(cs).reshape(n_species, n_cells).T\n",
"cs_r[:, indx_cl] += hcl # Cl-\n",
"cs_r[:, indx_h] += hcl # H+\n",
"cs1 = cs_r.T.reshape(n_cells * n_species) # Updated concentrations (after adding HCl)"
Expand Down Expand Up @@ -205,13 +203,13 @@
],
"source": [
"# Tranfer the modified species concentrations back to PhreeqcRM\n",
"phr.rm_species_concentrations2_module(cs1)\n",
"phr.rm_set_time(1.0)\n",
"phr.rm_set_time_step(1.0)\n",
"phr.rm.SpeciesConcentrations2Module(cs1)\n",
"phr.rm.SetTime(1.0)\n",
"phr.rm.SetTimeStep(1.0)\n",
"\n",
"# Run the simulation\n",
"start_time = time.time()\n",
"phr.rm_run_cells()\n",
"phr.rm.RunCells()\n",
"elapsed = time.time() - start_time\n",
"print(f\"Simulation completed in {elapsed:.2f} seconds\")"
]
Expand Down
20 changes: 9 additions & 11 deletions examples/ex1-titration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,10 @@
MiBiReMo - Example - Calculate calcite titration curve.

This script models the reaction:
CaCO3(s) + 2HCl(aq) → CaCl2(aq) + CO2(g) + H2O(l)
CaCO3(s) + 2HCl(aq) -> CaCl2(aq) + CO2(g) + H2O(l)
assuming chemical reactions at equilibrium.

Last revision: 20/02/2026
Last revision: 10/04/2026
"""

import time
Expand Down Expand Up @@ -51,10 +51,8 @@
n_species = len(species) # Number of species

# Initialize concentration vectors
cc = np.zeros(n_cells * n_comps, dtype=np.float64)
cs = np.zeros(n_cells * n_species, dtype=np.float64)
phr.rm_get_concentrations(cc)
phr.rm_get_species_concentrations(cs)
cc = phr.rm.GetConcentrations()
cs = phr.rm.GetSpeciesConcentrations()

# Set HCl concentrations
hcl = np.linspace(hcl_range[0], hcl_range[1], n_cells) # mol/L
Expand All @@ -64,18 +62,18 @@
indx_h = np.where(species == "H+")[0][0]

# Update species concentrations with HCl
cs_r = cs.reshape(n_species, n_cells).T
cs_r = np.array(cs).reshape(n_species, n_cells).T
cs_r[:, indx_cl] += hcl # Cl-
cs_r[:, indx_h] += hcl # H+
cs1 = cs_r.T.reshape(n_cells * n_species)

# Run simulation with added HCl
phr.rm_species_concentrations2_module(cs1)
phr.rm_set_time(1.0)
phr.rm_set_time_step(1.0)
phr.rm.SpeciesConcentrations2Module(cs1)
phr.rm.SetTime(1.0)
phr.rm.SetTimeStep(1.0)

start_time = time.time()
phr.rm_run_cells()
phr.rm.RunCells()
elapsed = time.time() - start_time
print(f"Simulation completed in {elapsed:.2f} seconds")

Expand Down
28 changes: 13 additions & 15 deletions examples/ex2-BTEX_dissolution.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

16 changes: 7 additions & 9 deletions examples/ex2-BTEX_dissolution.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

Models the dissolution of benzene and ethylbenzene from NAPL phase into aqueous phase.

Last revision: 20/02/2026
Last revision: 10/04/2026
"""

import time
Expand Down Expand Up @@ -54,10 +54,8 @@
n_species = len(species)

# Get initial concentrations
cc = np.zeros(n_cells * n_comps, dtype=np.float64)
cs = np.zeros(n_cells * n_species, dtype=np.float64)
phr.rm_get_concentrations(cc)
phr.rm_get_species_concentrations(cs)
cc = np.array(phr.rm.GetConcentrations())
cs = np.array(phr.rm.GetSpeciesConcentrations())

# Time step setup
dt = sim_duration / n_steps * 24 * 3600.0 # Convert days to seconds
Expand All @@ -82,12 +80,12 @@
time_vector[step] = current_time

# Run simulation step
phr.rm_set_time(current_time)
phr.rm_set_time_step(dt)
status = phr.rm_run_cells()
phr.rm.SetTime(current_time)
phr.rm.SetTimeStep(dt)
phr.rm.RunCells()

# Store results
status = phr.rm_get_concentrations(cc)
cc = np.array(phr.rm.GetConcentrations())
concentration_results[step, :] = cc[component_map]

elapsed = time.time() - start_time
Expand Down
23 changes: 11 additions & 12 deletions examples/ex3-BTEX_dissolution_transport_coupling.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

21 changes: 10 additions & 11 deletions examples/ex3-BTEX_dissolution_transport_coupling.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
- Equilibrium dissolution (NAPL as equilibrium phases)
- Kinetic dissolution (NAPL kinetic dissolution)

Last revision: 20/02/2026
Last revision: 10/04/2026
"""

import time
Expand All @@ -32,9 +32,9 @@
dt = 0.1 # Coupling time step (days)
domain_length = 100.0 # Length of domain (m)
n_threads = 6 # Threads for calculation (-1 for all CPUs)
dispersivity = 0.05 # Dispersivity (m²)
dispersivity = 0.05 # Dispersivity (m2)
velocity = 1.0 # Groundwater velocity (m/d)
diffusion_coeff = 0.0 # Molecular diffusion (m²/s)
diffusion_coeff = 0.0 # Molecular diffusion (m2/s)
sim_duration = 100.0 # Simulation duration (days)

# Physical properties
Expand All @@ -47,7 +47,7 @@
probe_location = 49.75 # Probe location (m) for results extraction

# Derived parameters
dispersion_coeff = dispersivity * velocity # Dispersion coefficient (m²/d)
dispersion_coeff = dispersivity * velocity # Dispersion coefficient (m2/d)
n_steps = int(sim_duration / dt) # Number of time steps
contaminant_spot = int(0.5 * n_cells / domain_length) # Contaminant spot size

Expand Down Expand Up @@ -90,8 +90,7 @@ def run_simulation(pqi_file, kinetic=False):
n_comps = len(components)

# Get initial concentrations
cc = np.zeros(n_cells * n_comps, dtype=np.float64)
phr.rm_get_concentrations(cc)
cc = np.array(phr.rm.GetConcentrations())

# Prepare monitoring
monitored_species = ["Benz", "Ethyl"]
Expand Down Expand Up @@ -121,10 +120,10 @@ def run_simulation(pqi_file, kinetic=False):
time_vector[step] = current_time

# 1) Reactive step
phr.rm_set_time(current_time * 3600 * 24) # Convert to seconds
phr.rm_set_time_step(dt * 3600 * 24)
phr.rm_run_cells()
phr.rm_get_concentrations(cc)
phr.rm.SetTime(current_time * 3600 * 24) # Convert to seconds
phr.rm.SetTimeStep(dt * 3600 * 24)
phr.rm.RunCells()
cc = np.array(phr.rm.GetConcentrations())
conc_matrix = cc.reshape((n_comps, n_cells)).T

# 2) Transport step
Expand All @@ -144,7 +143,7 @@ def run_simulation(pqi_file, kinetic=False):

# Update concentrations in PhreeqcRM
cc = conc_matrix.T.flatten()
phr.rm_set_concentrations(cc)
phr.rm.SetConcentrations(cc)

elapsed = time.time() - start_time
print(f"Simulation completed in {elapsed:.2f} seconds")
Expand Down
138 changes: 0 additions & 138 deletions mibiremo/irmresult.py
View file
Open in desktop

This file was deleted.

Binary file removed mibiremo/lib/PhreeqcRM.dll
View file
Open in desktop
Binary file not shown.
Binary file removed mibiremo/lib/PhreeqcRM.so
View file
Open in desktop
Binary file not shown.
Loading
Loading