Implementation of Sequential iterative approach (SIA) using IPHREEQC.

[Ivan K.]

I'm implementing SIA in my C++ code using IPHREEQC. As far as I understand, in the beginning of each time step I have to make a copy of the initial chemical states in the cells, say a copy of my IPhreeqc instance. There is no copy constructor, and I read that the best way of making such a copy is to dump one instance and load this dump file or string into another. I read the description of the dump operation, but it still is not clear to me how it can be run in a simple way (I'm not a chemist, but a specialist in numerical modeling of flow and transport). 
Can anyone show me an example of the way this procedure of copying is done?

[dlparkhurst]

You do not necessarily need to make a copy of the chemical state. You can simply update the solution compositions, run reactions, and return the updated solution compositions. The updated reactants (equilibrium phases, kinetics, surfaces, etc) are retained in the module for the next step.

The reaction module PhreeqcRM is designed to be a reaction module for transport codes and has built-in parallelization.

If you need to make a copy, say in case you need to restart the time step, then you can dump the state of the module into a string. That string can be used in another module with RUN_STRING, and the state of the original module will be copied into another module (provided you have initialized the second module with the same database and possibly database modifications). I expect the dump and read process will be slow.

Alternatively, you could use a set of COPY commands in the original module to copy the solutions and reactants to unused cell numbers. That would probably be faster, but you would have to know which reactants you need to copy. You could then reverse the process to go back to the original state of the module.

Finally, there is the beginning of a copy method in PHREEQC that has never been completed. It may be sufficiently complete to serve your purposes.

[Ivan K.]

Well, I'd like to stay for a while with IPhreeqc, not switching to PhreeqcRM as far as it already works normally within the noniterative scheme (SNIA).
For the iterative scheme I have to store the initial solution compositions.
What I would like to see is an example of dump and read from dump operations.

[dlparkhurst]

I think something like this.

SetDumpStringOn(true);
RunString("DUMP; -all");
mydump = new std::string(GetMyDumpString());
...
RunString(mydump.c_str())

If you decide to use PhreeqcRM, I will write a couple of methods to save and retrieve state more efficiently.

[Ivan K.]

Thank you very much, this seems to be very useful for me at the moment.

[Ivan K.]

Dear David,
I shifted to PHREEQCRM in my tranport code (previously IPHREEQC was used). Are there any new functions already in PHREEQCRM which allow to store the initial chemical state and to copy it in order to perform SIA?

[dlparkhurst]

Depends on which versions we are comparing, but functions were added to get names of all of the various types of reactants--GetEquilibriumPhases, GetExchangeNames, GetExchangeSpecies, and others, which would allow writing SELECTED_OUTPUT/USER_PUNCH blocks to retrieve states (moles) of reactants. You could use _MODIFY datablocks to reset the states with these values. Doable, but probably not ideal.

There is a capability to copy a Phreeqc instance with the C++ method operator=. (There is not a capability to copy an IPhreeqc or PhreeqcRM instance.) The operator= has not been widely tested, so the first problem is whether it works correctly.

When using OpenMP for parallelization (or if you are using neither OpenMP nor MPI), you can get the worker IPhreeqcPhast instances with PhreeqcRM:GetWorkers(). There are 1 or more worker instances. You can get a pointer to a Phreeqc instance within the IPhreeqcPhast instance with IPhreeqcPhast:Get_PhreeqcPtr(). I think, you should be able to copy each worker Phreeqc instance with something like:

Phreeqc *w0 = phreeqc_rm->GetWorkers()[0]->Get_PhreeqcPtr();

and copy it back with

phreeqc_rm.GetWorkers()[0]->Get_PhreeqcPtr() = w0;

So, if you copy all of the worker Phreeqc instances, you should have copied the current state of PhreeqcRM. You probably need to exclude load balancing [SetRebalanceFraction(0.0)] to keep the set of cells in each worker instance the same.

I don't think it is possible to do the same with MPI without additional PhreeqcRM methods. It also would not be possible to do the same with C or Fortran without additional PhreeqcRM and wrapper methods.

[dlparkhurst]

Erratum. The last post had the wrong level of indirection.

Phreeqc w0;
w0 = *phreeqc_rm->GetWorkers()[0]->Get_PhreeqcPtr();

and copy it back with

*phreeqc_rm.GetWorkers()[0]->Get_PhreeqcPtr() = w0;

[dlparkhurst]

I'll add a capability to save the state of the chemistry in the next release with three methods:

StateSave(i);
StateApply(i);
StateDelete(i);