Coupling IPHREEQC/PHREEQCRM with a multi-phase flow solver

[ahmadreza_shojaee]

Hi

I am currently working on a research project where I have coupled IPHREEQC with a multiphase flow solver in MATLAB to simulate reactive transport in porous media. In some cases, the number of cells becomes quite large, and performing batch calculations for each cell using IPHREEQC significantly increases the computational cost.

I would like to ask whether PHREEQCRM would be a more suitable option in this case to reduce computational time, given its support for parallel reaction calculations. Additionally, is it possible to interface PHREEQCRM with MATLAB, either directly or through a wrapper, to maintain integration with my current simulation framework?

Any guidance or recommendations would be greatly appreciated.

Kind regards,
Ahmad

[dlparkhurst]

Yes, PhreeqcRM would be a better choice.

There is a C interface to PhreeqcRM, which should make it suitable for use with Matlab. I think it has already been done here: https://github.com/simulkade/PhreeqcMatlab.

[ahmadreza_shojaee]

Dear David,

Thank you for your reply.

As a follow-up to my initial message, I wanted to clarify that in my current setup, I am using IPhreeqc to perform independent single-cell calculations across the domain, without any transport (advection or diffusion) modeled directly. My main objective now is to improve computational efficiency.

Given this context, I?m wondering whether PHREEQCRM can be configured purely for parallel single-cell chemistry calculations, without invoking transport processes. Would setting the transport coefficients to zero be the correct approach?

If possible, I would also appreciate any examples or references demonstrating such a use case in PHREEQCRM.

Kind regards,
Ahmad

[dlparkhurst]

I don't know quite enough about what you are doing. How do the calculations change from one calculation for a "cell" to the next?

PhreeqcRM is set up to change the solution concentrations from one time step to the next with the SetConcentrations method. New solution concentrations after the chemistry calculation (RunCells) are retrieved with GetConcentrations, and other results related to the calculations can be retrieved through GetValue (BMIPhreeqcRM) or selected output. If your setup conforms to this sequence, then  PhreeqcRM would work well. It is possible to change conditions other than solution (and a

If more than solution concentrations need to be changed, then PhreeqcRM may not be suitable, or you may need to use it in a different way. PhreePlot may be similar to your needs. It calculates a grid of conditions, say pH vs pO2 for example. In this case, a string was needed to define each calculation, and it was more efficient to use the IPhreeqc instances (OMP version) individually. The method GetIPhreeqcId lets you obtain each IPhreeqc instance, and then the calculation was parallelized in an OMP do loop using IPhreeqc methods RunAccumulated and GetSelectedOutput. It was still convenient to use PhreeqcRM for initialization and control of multiple IPhreeqc instances rather than simply parallelizing PhreePlot itself.

If you explain what needs to change for each calculation, I might be able to give a little more insight.

[ahmadreza_shojaee]

Dear David,

Thank you very much for your detailed explanation and your willingness to assist.

To give you a clearer picture of my workflow: I am modelling multiphase flow in porous media using a reservoir simulator. For example, consider a case with 100 cells. After each time step of the flow simulation, I obtain updated values of pressure, temperature, mass of water, mineral amounts, and species concentrations for all cells.

Using this data, I construct 100 individual PHREEQC input strings (one for each cell) and run them sequentially in MATLAB using IPhreeqc. Each calculation is independent, and the purpose is to evaluate geochemical reactions within each cell based on its local conditions. As the number of cells increases, this sequential batch processing becomes a major computational bottleneck.

I am therefore exploring whether PhreeqcRM would allow me to run all these cell-based calculations in parallel, using the provided data, without needing to manually process each input string one at a time. My goal is to speed up the simulation by leveraging parallel execution while still treating the chemistry as independent across cells (i.e., no transport between them).

If this type of workflow is compatible with PhreeqcRM, I would greatly appreciate any guidance on how best to implement it.

Kind regards,
Ahmad

[dlparkhurst]

You should look at PhreeqcRM. It is designed to accept solution concentrations, temperature, pressure, and density to convert data to PHREEQC concentration units and run geochemical calculations.

Normally, the moles of minerals is calculated and maintained by PhreeqcRM.

[ahmadreza_shojaee]

Dear David,

It seems that setting the dispersivity and diffusion coefficients to zero, along with using closed-system boundary conditions, works for running independent cell calculations in PhreeqcRM.

My model includes RATES and KINETICS definitions, and it runs successfully. However, I?m currently unsure how to output specific variables?such as calcium (Ca) concentration?for each cell over all time steps. Could you please advise on how this can be achieved in PhreeqcRM?

Any guidance would be greatly appreciated.

Kind regards,
Ahmad

[dlparkhurst]

Documentation for PhreeqcRM is in the folder Doxygen/HTML in the directory in which you extracted PhreeqcRM.

You can use the method GetOutputVars to determine those variables that you can get with the method GetValue when using BMIPhreeqcRM.

A BMIPhreeqcRM instance also contains all of the methods of PhreeqcRM. The PhreeqcRM methods to extract data use SelectedOutput methods. You can define a SELECTED_OUTPUT/USER_PUNCH data block with the input file (run by the workers) and you can then extract the data with GetSelectedOutput and other associated methods. You must have used similar methods when you were using IPhreeqc.

There are quite a few posts on PhreeqcRM in this forum that might be useful.

[dlparkhurst]

BTW, PhreeqcRM does not know anything about transport, so it does not have methods to set dispersivity and it does not apply boundary conditions. It does no transport, it is only intended to run chemistry for each cell.

I have no idea what you are doing. Are you using TRANSPORT in Phreeqc? That is not the intended use of PhreeqcRM.

PhreeqcRM has been applied to more than a dozen transport models, so I think it should work for you in the way that it is intended.

[ahmadreza_shojaee]

Thank you, David for the comprehensive response.

To be more specific, I do not need any transport in my calculations. I am looking for a way to run chemistry for each single cell in parallel. Currently, for large models (say with 10000 cells), I have run time problem as my model is running each cell separately. I am looking for a way to run the chemistry of those cells without any transport. what if I use the following way:

Code: [Select]

Rate #where I define the rate models
#Input data for cell 1
Solution 1
Equilibrium_Phases 1
Gas_Phase 1
Kinetics 1
#Input data for cell 2
Solution 2
Equilibrium_Phases 2
Gas_Phase 2
Kinetics 2
....................
....................
....................
#Input data for cell 100
Solution 100
Equilibrium_Phases 100
Gas_Phase 100
Kinetics 100

Does this approach work to run cells at the same time?

If yes, how can I access the data for each cell?

Kind regards,
Ahmad

[dlparkhurst]

No, that does not run all of the cells in parallel. It is input for PHREEQC, it is not using PhreeqcRM, which has the parallel capabilities.

I don't think you have understood anything that I have said.

Here is a simple example of the use of PhreeqcRM, which is in the Tests directory of the PhreeqcRM installation. The cells are initialized, and run in parallel, the concentrations are extracted with GetConcentrations, where advection is calculated you would do your transport calculation, transported concentrations are sent to PhreeqcRM with SetConcentrations, and chemistry is run again (in parallel) with RunCells. After the last time step, concentrations are extracted from PhreeqcRM.

Please study this and other PhreeqcRM examples, and study the documentation found in Doxygen/HTML.

Code: [Select]

#if defined(USE_MPI)
#include <mpi.h>
#endif
#include <stdlib.h>
#include <iostream>
#include <string>
#include <vector>
#include "PhreeqcRM.h"
#include "IPhreeqc.hpp"
#include "IPhreeqcPhast.h"

void simpleadvection_cpp(std::vector<double>& c, std::vector<double> bc_conc, int ncomps, int nxyz, int dim);

int SimpleAdvect_cpp()
{
// Based on PHREEQC Example 11
try
{
// --------------------------------------------------------------------------
// Create PhreeqcRM
// --------------------------------------------------------------------------
int nxyz = 20;
#ifdef USE_MPI
// MPI
PhreeqcRM phreeqc_rm(nxyz, MPI_COMM_WORLD);
MP_TYPE comm = MPI_COMM_WORLD;
int mpi_myself;
if (MPI_Comm_rank(MPI_COMM_WORLD, &mpi_myself) != MPI_SUCCESS)
{
exit(4);
}
if (mpi_myself > 0)
{
phreeqc_rm.MpiWorker();
return EXIT_SUCCESS;
}
#else
// OpenMP
int nthreads = 3;
PhreeqcRM phreeqc_rm(nxyz, nthreads);
#endif
// Set properties
IRM_RESULT status = phreeqc_rm.SetComponentH2O(false);
phreeqc_rm.UseSolutionDensityVolume(false);
// Open files
status = phreeqc_rm.SetFilePrefix("SimpleAdvect_cpp");
phreeqc_rm.OpenFiles();
// Set concentration units
status = phreeqc_rm.SetUnitsSolution(2);           // 1, mg/L; 2, mol/L; 3, kg/kgs
status = phreeqc_rm.SetUnitsExchange(1);           // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// Set conversion from seconds to user units (days)
double time_conversion = 1.0 / 86400;
status = phreeqc_rm.SetTimeConversion(time_conversion);
// Set initial porosity
std::vector<double> por;
por.resize(nxyz, 0.2);
status = phreeqc_rm.SetPorosity(por);
// Set cells to print chemistry when print chemistry is turned on
std::vector<int> print_chemistry_mask;
print_chemistry_mask.resize(nxyz, 1);
status = phreeqc_rm.SetPrintChemistryMask(print_chemistry_mask);
int nchem = phreeqc_rm.GetChemistryCellCount();
// --------------------------------------------------------------------------
// Set initial conditions
// --------------------------------------------------------------------------
// Set printing of chemistry file
status = phreeqc_rm.SetPrintChemistryOn(false, true, false); // workers, initial_phreeqc, utility
// Load database
status = phreeqc_rm.LoadDatabase("phreeqc.dat");
// Run file to define solutions and reactants for initial conditions, selected output
status = phreeqc_rm.RunFile(true, true, true, "advect.pqi");
// Clear contents of workers and utility
std::string input = "DELETE; -all";
status = phreeqc_rm.RunString(true, false, true, input.c_str());
// Determine number of components to transport
int ncomps = phreeqc_rm.FindComponents();
// Get component information
const std::vector<std::string>& components = phreeqc_rm.GetComponents();
for (int i = 0; i < ncomps; i++)
{
std::ostringstream strm;
strm.width(10);
strm << components[i] << "\n";
phreeqc_rm.OutputMessage(strm.str());
}
phreeqc_rm.OutputMessage("\n");
// Set array of initial conditions
std::vector<int> ic1;
ic1.resize(nxyz * 7, -1);
std::vector<double> f1;
f1.resize(nxyz * 7, 1.0);
for (int i = 0; i < nxyz; i++)
{
ic1[i] = 1;              // Solution 1
ic1[nxyz + i] = -1;      // Equilibrium phases none
ic1[2 * nxyz + i] = 1;     // Exchange 1
ic1[3 * nxyz + i] = -1;    // Surface none
ic1[4 * nxyz + i] = -1;    // Gas phase none
ic1[5 * nxyz + i] = -1;    // Solid solutions none
ic1[6 * nxyz + i] = -1;    // Kinetics none
}
status = phreeqc_rm.InitialPhreeqc2Module(ic1);
// Initial equilibration of cells
double time = 0.0;
double time_step = 0.0;
std::vector<double> c;
c.resize(nxyz * components.size());
status = phreeqc_rm.SetTime(time);
status = phreeqc_rm.SetTimeStep(time_step);
status = phreeqc_rm.RunCells();
status = phreeqc_rm.GetConcentrations(c);
// --------------------------------------------------------------------------
// Set boundary condition
// --------------------------------------------------------------------------
std::vector<double> bc_conc;
std::vector<int> bc1;
int nbound = 1;
bc1.resize(nbound, 0);                      // solution 0 from Initial IPhreeqc instance
status = phreeqc_rm.InitialPhreeqc2Concentrations(bc_conc, bc1);
// --------------------------------------------------------------------------
// Transient loop
// --------------------------------------------------------------------------
int nsteps = 10;
std::vector<double> temperature, pressure;
temperature.resize(nxyz, 20.0);
pressure.resize(nxyz, 2.0);
phreeqc_rm.SetTemperature(temperature);
phreeqc_rm.SetPressure(pressure);
time_step = 86400.;
status = phreeqc_rm.SetTimeStep(time_step);
for (int steps = 0; steps < nsteps; steps++)
{
// Transport calculation here
{
std::ostringstream strm;
strm << "Beginning transport calculation             " << phreeqc_rm.GetTime() * phreeqc_rm.GetTimeConversion() << " days\n";
strm << "          Time step                         " << phreeqc_rm.GetTimeStep() * phreeqc_rm.GetTimeConversion() << " days\n";
phreeqc_rm.LogMessage(strm.str());
phreeqc_rm.SetScreenOn(true);
phreeqc_rm.ScreenMessage(strm.str());
}
simpleadvection_cpp(c, bc_conc, ncomps, nxyz, nbound);
// Transfer data to PhreeqcRM for reactions
bool print_selected_output_on = (steps == nsteps - 1) ? true : false;
bool print_chemistry_on = (steps == nsteps - 1) ? true : false;
status = phreeqc_rm.SetSelectedOutputOn(print_selected_output_on);
status = phreeqc_rm.SetPrintChemistryOn(print_chemistry_on, false, false); // workers, initial_phreeqc, utility
status = phreeqc_rm.SetConcentrations(c);         // Transported concentrations
time += time_step;
status = phreeqc_rm.SetTime(time);
// Run cells with transported conditions
{
std::ostringstream strm;
strm << "Beginning reaction calculation              " << time * phreeqc_rm.GetTimeConversion() << " days\n";
phreeqc_rm.LogMessage(strm.str());
phreeqc_rm.ScreenMessage(strm.str());
}
status = phreeqc_rm.RunCells();
// Transfer data from PhreeqcRM for transport
status = phreeqc_rm.GetConcentrations(c);
}
// Clean up
status = phreeqc_rm.CloseFiles();
status = phreeqc_rm.MpiWorkerBreak();
}
catch (PhreeqcRMStop)
{
std::string e_string = "SimpleAdvect_cpp failed with an error in PhreeqcRM.";
std::cerr << e_string << std::endl;
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD, 1);
#endif
return IRM_FAIL;
}
#if 0
catch (...)
{
std::string e_string = "SimpleAdvect_cpp failed with an unhandled exception.";
std::cerr << e_string << std::endl;
#ifdef USE_MPI
MPI_Abort(MPI_COMM_WORLD, 1);
#endif
return IRM_FAIL;
}
#endif
return EXIT_SUCCESS;
}
void
simpleadvection_cpp(std::vector<double>& c, std::vector<double> bc_conc, int ncomps, int nxyz, int dim)
{
for (int i = nxyz - 1; i > 0; i--)
{
for (int j = 0; j < ncomps; j++)
{
c[j * nxyz + i] = c[j * nxyz + i - 1];                 // component j
}
}
// Cell zero gets boundary condition
for (int j = 0; j < ncomps; j++)
{
c[j * nxyz] = bc_conc[j * dim];                                // component j
}
}




[ahmadreza_shojaee]

Dear David,

Thank you for the clarification, and apologies for the confusion.

I now see where I misunderstood. I?ll review the example in the Tests directory and dive into the PhreeqcRM documentation to better understand the proper workflow.

Thanks again for your guidance.

Kind regards,
Ahmad