Conceptual Models > Incorporation PHREEQC in programming languages

SetConcentrations correctly in 2D model


Hello David,

I am trying to model reactions in a 2D domain using PhreeqcRM. My 1D transport model works fine with PhreeqcRM, while my 2D model works without PhreeqcRM. However, I get an “unknown exception” error (I am using PhreeqcMatlab – with the advection_cpp file as my guide) when I couple the 2D model with PhreeqcRM. I know there may likely be a syntax problem going on, but I can’t figure out how to fix it. That is, I am not sure in what form I need to set the concentration before Runcells().
•   For instance, in 1D, the initial concentration at time_step 0 is an n x ncomps size array, where n = number of cells, and ncomps = number of components. Should this remain the same in 2D, with n = total number of 2D cells? Should this be the same size for the rest of the time_steps?
•   What about the mapping? I know grid2chem is 0:nxyz-1 in 1D. Not sure how it’ll look like in 2D - or if it really does matter?
•   What about ic1, ic2, f1, bc1, bc2, bc_f1, bc_conc? Are their array sizes expected to change from what they are in 1D?

In a nut shell, what is the key to setting the 2D concentration (or how should the set concentration look like) before employing Runcells()? This is where my code fails.

Thanks for any leads.

For starters, when you create an instance of PhreeqcRM, you must define nxyz, which is the total number of cells in your transport model. If it is a rectangular grid, nxyz = nrows * ncolumns.

You then must decide on a numbering system, 1 to nxyz, for the transport grid. For Fortran, the numbering is usually 1,1 (1); 2,1 (2); 3,1 (3), ... . For C/C++ it is usually 1,1 (1); 1,2 (2); 1,3 (3) ... . I don't know the order of 2D arrays in Matlab. You will have to make sure the component concentrations  end up in the right order to be able to transport each component. Whatever you pick, by default, there are nxyz chemistry cells with the same numbering as the transport cells. If there are inactive cells, or symmetry, you can use a mapping to limit the number of cells for which chemistry is run (nchem <= nxyz).

The arguments to SetConcentrations and GetConcentrations are always arrays that are nxyz x ncomps, regardless of whether or not there are inactive cells. The concentration data is translated from the nxyz transport cells to the nchem chemistry cells (SetConcentrations) and from the chemistry cells back to the transport cells (GetConcetrations).

Similarly, ic1, ic2, ... arrays are nxyz x 7. 


[0] Message Index

Go to full version