Unexpected kinetic reactions in PhreeqcRM

[hwjung]

Dear developers,

I have encountered an issue in my simulation where a kinetic reactant unexpectedly appears in the second grid cell, despite no kinetic reaction being assigned to this cell. This behavior does not occur in any other grid cells. Below, I have included the model outputs showing the mass distribution for the first, second, third, and fifth grid cells. The kinetic reactant should be placed on the fifth grid cells:

Cell #:                  1              2              3            5
pH:                   5.8            5.8            5.8            5.8
m_Ca+2 (mol/kgw):     4.00001e-08    4.00001e-08    4.00001e-08    4.00001e-08
m_Cl- (mol/kgw):      1.65223e-06    1.65223e-06    1.65223e-06    1.65223e-06
m_HCO3- (mol/kgw):    8.68809e-09    8.68809e-09    8.68809e-09    8.68809e-09
m_CO3-2 (mol/kgw):    2.58264e-13    2.58264e-13    2.58264e-13    2.58264e-13
si_Calcite:          -11.5112       -11.5112       -11.5112       -11.5112
k_Calcite:            0              3.38e-09       0              3.38e-09

I have attached the minimal C++ code and the input.in file that reproduce the issue on my Mac (M1 Max, macOS Sonoma 14.4.1, gcc-13, -std=c++11).
Any help would be greatly appreciated.

Best,
hwjung

Code: [Select]

#include "PhreeqcRM.h"
#include <fstream>

int main(int argc, char **argv) {
    int subNxyz = 150;
    const int nthreads = 1;

    std::unique_ptr<PhreeqcRM> phreeqc_rm(new PhreeqcRM(subNxyz, nthreads));
    int status = phreeqc_rm->SetComponentH2O("true");
    status = phreeqc_rm->SetFilePrefix("main");
    status = phreeqc_rm->OpenFiles();

    std::vector<double> P0(subNxyz,1.);
    std::vector<double> T0(subNxyz,25.);
    std::vector<double> rho0(subNxyz,1.);
    std::vector<double> por0(subNxyz,1.);
    std::vector<double> sat0(subNxyz,1.);
    std::vector<double> REV0(subNxyz,1.);

    status = phreeqc_rm->SetRepresentativeVolume(REV0);
    status = phreeqc_rm->SetPorosity(por0);
    status = phreeqc_rm->SetSaturation(sat0);
    status = phreeqc_rm->SetTime(0.);
    status = phreeqc_rm->SetDensity(rho0);
    status = phreeqc_rm->SetTemperature(T0);
    status = phreeqc_rm->SetPressure(P0);
    status = phreeqc_rm->SetSelectedOutputOn(true);
    status = phreeqc_rm->SetScreenOn(false);
    status = phreeqc_rm->SetPrintChemistryOn(false, false, false);
    status = phreeqc_rm->SetUnitsSolution(2);

    status = phreeqc_rm->LoadDatabase("phreeqc.dat");
    if (status != 0) { phreeqc_rm->ErrorHandler(status, "PhreeqcRM LoadDatabase failed"); }

    std::string PHQ_in;
    std::ifstream file("input.in");
    std::stringstream buffer;
    buffer << file.rdbuf();
    PHQ_in = buffer.str();
    status = phreeqc_rm->RunString(true, true, true, PHQ_in);
    if (status != 0) { phreeqc_rm->ErrorHandler(status, "PhreeqcRM RunString failed"); }

    int ncomps = phreeqc_rm->FindComponents();

    std::vector<int> ic1(subNxyz*7,-1);
    for (int iT=0; iT<subNxyz; ++iT) {
            ic1[iT] = 1; // Solution: 1
            if (iT == 4 ) {
                ic1[6*subNxyz+iT] = 1; // Kinetics: 1
            }
    }
    status = phreeqc_rm->InitialPhreeqc2Module(ic1);
    status = phreeqc_rm->RunCells();
    if (status != 0) { phreeqc_rm->ErrorHandler(status, "PhreeqcRM RunCells failed"); }

    std::vector<double> c0(subNxyz * ncomps, 0);
    status = phreeqc_rm->GetConcentrations(c0);
    std::vector<double> so2;
    phreeqc_rm->GetSelectedOutput(so2);

    std::vector<int> headloc={6}, solidloc;
    for (int i = 8; i < phreeqc_rm->GetSelectedOutputColumnCount()-1; ++i) { headloc.push_back(i); }
    int outputCount = headloc.size();
    if (outputCount==1) {outputCount = 0;}
    std::vector<std::string> headings;
    headings.resize(outputCount);
    for (int j = 0; j < outputCount; j++) {
        phreeqc_rm->GetSelectedOutputHeading(headloc[j], headings[j]);
        if (headings[j].substr(0, 2) == "k_") {
            solidloc.push_back(j);
        }
    }

    for (int iC=0; iC<outputCount; ++iC) {
        std::cout << headings[iC] << " " << so2[headloc[iC]*subNxyz+0] << " " << so2[headloc[iC]*subNxyz+1] <<  " " << so2[headloc[iC]*subNxyz+2] <<  " " << so2[headloc[iC]*subNxyz+4] << std::endl;
    }
    std::cout<< std::endl;

    return 0;
}


Code: [Select]

SOLUTION 1
    pH   5.8
    temp 25.
    unit mol/L
    Ca   4e-8
    Cl   6.42e-5 charge
    C    3.96e-8
end

KINETICS 1
Calcite
-tol   1e-15
-m0    3.38e-09
-m     3.38e-09
-parms 5.012e-05 3.311e-08 1.585e-10
-time  0.0897129

RATES
Calcite
   -start
1  rem PARM(1) - PARM(3) = rate constants in mol/cm2/s
10 si_cc = SI("Calcite")
20 rate = (PARM(1) * ACT("H+") + PARM(2) * ACT("HCO3-") + PARM(3))
30 rate = 1e6 * rate * (1 - 10^(si_cc))
40 moles = rate * TIME
50 SAVE moles

SELECTED_OUTPUT 1
    -file                 so.txt
    -pH                   true
    -kinetic_reactants    Calcite
    -saturation_indices   Calcite
    -molalities           Ca+2 Cl- HCO3- CO3-2


[dlparkhurst]

I'm not feeling energetic this morning, so I'm not planning to debug your code. I'll give you a couple possibilities.

I assume your indexing is correct, but in PhreeqcRM 3.8.0, there are new methods to set the initial reactants individually, instead of the 7 x nxyz array. You can use rm.InitialKinetics2Module(ik), where ik is an integer array of length nxyz. This approach would avoid any indexing error.

The other possibility is that you have run your input file through the workers, the initial phreeqc, and the utility instances, so KINETICS 1 is defined in all of the IPhreeqc workers. You should probably remove all solutions and reactants before you populate the grid with InitialPhreeqc2Module. It is a bit of a flaw in the initialization strategy that you may need to define SELECTED_OUTPUT or aqueous model changes in the input file, but the SOLUTION and reactant definitions may be hanging around in the worker cells. I recommend using the following before you populate the cells (it is in all of the examples).

Code: [Select]

std::string input = "DELETE; -all";
status = phreeqc_rm.RunString(true, false, true, input.c_str());


Code: [Select]

//
// Four ways to set initial conditions
//
// 1. Define mixing
//std::vector<int> ic1, ic2;
//ic1.resize(nxyz * 7, -1);
//ic2.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, ic2, f1);

// 2. No mixing is defined, so the following is equivalent
// status = phreeqc_rm.InitialPhreeqc2Module(ic1);

// 3. Simplest way to set these initial conditions
//
{
std::vector<int> solutions(nxyz, 1);
phreeqc_rm.InitialSolutions2Module(solutions);
std::vector<int> exchanges(nxyz, 1);
phreeqc_rm.InitialExchanges2Module(exchanges);
}


[hwjung]

Dr. Parkhurst,

Thanks for your recommendation regarding the "DELETE; -all" command. I wasn't entirely sure of its exact function before, but with your advice, my code is now working perfectly. Thank you so much for your assistance.