Using PhreeqcRM in Fortran

[Zhaoyang]

Dear all,


In order to couple with a Fortran-based transport code, I have compiled PhreeqcRM.

In the PhreeqcRM.sln, there are 7 projects: ALL-BUILD, FortranAdvect, INSTALL, PhreeqcRM, RUN-TESTS, TestRM, ZERO-CHECK. When I build the RUN-TESTS, all modules based on different languages (C, C++, Fortran) will run. Can I separate them due to focusing on Fortran? In addition, I am wondering whether I need to do in next is to replace advection_f90.F90 with my transport code.

Thanks for your help in advance.

Bests,
Zhaoyang

[dlparkhurst]

I think the only option you have is whether to test Fortran, which you want to do. With the Visual Studio solution that is provided, you simply want to make the PhreeqcRM library and test that it works correctly.

To link with your transport code, you need to create a new Visual Studio project. You would include your code (with calls to the PhreeqcRM subroutines), the Fortran file RM_interface.F90 , and you will need to link to the PhreeqcRM library that you have compiled.

[Zhaoyang]

Dear Dr. Parkhurst,

Thanks for your fast reply. Due to running the advection test successfully, I am sure that the PhreeqcRM library and test work correctly .

By the advection test, I found that the input file (advect.pqi) for PhreeqcRM is the same as that in Phreeqc. Does it mean that the input file for PhreeqcRM in coupling model is still the same as that in Phreeqc?

In addition, I want to obtain the following information from the advection test:
-How to call PhreeqcRM in the transport code, especially two functions (i.e., RM_SetConcentrations and RM_GetConcentrations) which are used to transfer concentrations between PhreeqcRM and transport code?
-How to link the transport code to the PhreeqcRM library that I have compiled?
Can I obtaion these two information by the advectionn test?

Thanks for your help in advance.

Bests,
Zhaoyang

[dlparkhurst]

My brother and sister are Dr. Parkhurst, but I am not. No Dr. please.

Yes, the strategy with PhreeqcRM is to use a PHREEQC input file to define solutions and reactions for initial and boundary conditions.

Advection_f90.F90 shows how to use many of the methods of PhreeqcRM, including RM_GetConcentrations and RM_SetConcentrations. You must transport total H, total O, (H2O, excess H, and excess O), plus every other element in solution. A Fortran array is used to transfer the concentrations from PhreeqcRM (RM_GetConcentrations) to your Fortran and to PhreeqcRM (RM_SetConcentrations) from your Fortran. All methods are documented in the phreeqcrm-3.6.2-15100/Doxygen/html directory of the distribution or online at https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm/index.htm

You will have to develop a Visual Studio solution for your simulator. You will need to link to PhreeqcRM.lib that you compiled to your code. The INSTALL option of PhreeqcRM Visual Studio will place the library in a location of your choice (although you must specify that choice when you run CMake, see doc/README.txt for details). In the VS solution for your simulator, you must specify the directory where the library is located and the library itself in the link options for your project.

[Zhaoyang]

Dear Parkhurst,

Thanks for your kindly help.

Now, I am understanding the code of advection_f90.F90. They are many return values assigned to “status”, e.g., status = RM_SetUnitsSolution(id, 2), status = RM_SetUnitsPPassemblage(id, 1), status = RM_SetUnitsExchange(id, 1) and so on. In this way, the value of “status” varies with each assignment. Could you explian this for me? Thanks a lot.

Best regards,
Zhaoyang

[dlparkhurst]

The PhreeqcRM calls are functions that return a value, so you must use a variable like status to receive the return value from each function.

The return values are usually error codes. If the value is zero, then the function executed successfully. Negative values indicate various error conditions. For simplicity, most return values were not checked, but to be careful, it would be good to check each return value for successful completion, and write an error message if there is a failure.

The following decodes an error if the return value is not zero:

Code: [Select]

if (status < 0) status = RM_DecodeError(id, status)

Something like the following would extract an error message from PhreeqcRM:

Code: [Select]

  ! Demonstration of error handling if ErrorHandlerMode is 0
  if (status .ne. 0) then
     l = RM_GetErrorStringLength(id)
#ifdef FORTRAN_2003
     allocate (character(len=l) :: errstr)
#endif
     write(*,*) "Start of error string: "
     status = RM_GetErrorString(id, errstr)
     write(*,"(A)") trim(errstr)
     write(*,*) "End of error string."
#ifdef FORTRAN_2003
     deallocate(errstr)
#endif
     status = RM_Destroy(id)
     stop
  endif

Ideally, you would write a subroutine that could be used to check each return code, and print the error message on failure. Even without checking every return code, the value of the return variable can help when you are debugging.

[Zhaoyang]

I really appreciate for your help, and I am getting more and more familiar with PhreeqcRM.
 
I don't understand a remark in the species_f90.F90, namely "Based on PHREEQC Example 11, transporting species rather than components". What does the "transporting species rather than components" mean? What is the difference between the advection_f90.F90 and species_f90.F90? Thanks a lot.

Cheers,
Zhaoyang

[dlparkhurst]

For diffusion calculations, species rather than components need to be transported, because each species may diffuse at a different rate. There are fewer components than species, possibly by a factor of 10, so it is a faster calculation if transporting each component is sufficient for your needs (likely in an advective system) than transporting each species (diffusion or electro-diffusion system).

As an example, for a pure water system, advection_f90 would transport H, O, and charge, whereas species_f90 would transport OH-, H+, H2O, H2(aq), and O2(aq).

Code: [Select]

----------------------------Distribution of species----------------------------

                                               Log       Log       Log    mole V
   Species          Molality    Activity  Molality  Activity     Gamma   cm³/mol

   OH-             1.013e-07   1.012e-07    -6.995    -6.995    -0.000     -4.14
   H+              1.001e-07   1.000e-07    -7.000    -7.000    -0.000      0.00
   H2O             5.551e+01   1.000e+00     1.744     0.000     0.000     18.07
H(0)          1.416e-25
   H2              7.079e-26   7.079e-26   -25.150   -25.150     0.000     28.61
O(0)          0.000e+00
   O2              0.000e+00   0.000e+00   -42.080   -42.080     0.000     30.40


[Zhaoyang]

In the set initial conditions of the advection_f90.F90, I am not clear for the array c in RM_GetConcentrations. Which function will return the array c? Does this array return from RM_RunFile? I am also not clear for the arry bc_conc in RM_InitialPhreeqc2Concentrations for seting the boundary condition. Which function will return the array bc_conc?

Code: [Select]

allocate(c(nxyz, ncomps)) !size of c: nxyz*ncomps
  status = RM_SetTime(id, time)
  status = RM_SetTimeStep(id, time_step)
  status = RM_RunCells(id) !Runs a reaction step for all of the cells in the reaction module.
  status = RM_GetConcentrations(id, c)

  nbound = 1  ! The first cell
  allocate(bc1(nbound), bc2(nbound), bc_f1(nbound))
  allocate(bc_conc(nbound, ncomps)) 
  bc1 = 0           ! solution 0 from Initial IPhreeqc instance
  bc2 = -1          ! no bc2 solution for mixing
  bc_f1 = 1.0       ! mixing fraction for bc1
  !Fills an array (bc_conc) with concentrations from solutions in the InitialPhreeqc instance.
  status = RM_InitialPhreeqc2Concentrations(id, bc_conc, nbound, bc1, bc2, bc_f1)

In addition, the initial and boundary conditions are obtained by implementing a few functions like RM_RunFile in the advection example. As mentioned before, I want to copule the 3D transport model with PhreeqcRM. For this coupling model, I will set the initial and boundary conditions for the 3D transport model. In this way, do I need to set the initial and boundary conditions for PhreeqcRM?

Thanks very much for your help in advance.

[dlparkhurst]

RM_GetConcentrations returns c.

Code: [Select]

status = RM_GetConcentrations(id, c)

The strategy is to use RM_RunFile to define a set of solutions and reactants in the "initial" iphreeqc instance that will be the initial and boundary conditions for the model. However, you need to populate the cells (in the "worker" instance(s)) by transferring the definitions from the initial to the worker instance with the method RM_InitialPhreeqc2Module or RM_InitialPhreeqcCell2Module. Each of these methods lets you put a specific combination of solution and reactants in each cell of the model. For example, ignoring other reactants for now, if you want a uniform solution composition throughout your 3D model, you could define one solution (say SOLUTION 1, in a file named initial.pqi) with RM_RunFile(id, 0, 1, 0, "initial.pqi"), and then copy it to all of the cells with RM_InitialPhreeqc2Module.

Concentrations to be used for boundary conditions (say you define another solution, SOLUTION 2, in initial.pqi) can be retrieved from the initial iphreeqc instance with RM_InitialPhreeqc2Concentrations.

[Zhaoyang]

However, I just checked the RM_GetConcentrations in https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqcrm/classphreeqcrm.html#a2880083aff4c3709dccdbfacd887cd0a. This told me that RM_GetConcentrations returns status (0 is success, negative is failure), not c.
Unfortunately, I didn't understand what you said about the initial and boundary conditions setup in my coupling model. Does it mean that I need to set the initial and boundary conditions for both transport and reactive code? Thanks tons.

[dlparkhurst]

The concentrations are returned through the second argument of RM_GetConcentrations(id, conc), in this example, the array conc.

Sorry, I can't explain it any differently. You set the initial concentrations for the module with RM_InitialPhreeqc2Module, and you get the concentrations for your transport model with RM_GetConcentrations. If you still don't understand, use the debugger to step through advection_f90.F90.

[dlparkhurst]

Fortran allows INTENT of arguments to be declared IN, OUT, or INOUT, meaning the argument is input to the function, output from the function or either input or output. For IN, the argument cannot be changed in the function, for OUT, the values entering the function are assumed undefined, and can only be set. For INOUT, a variable can be input, output, or both. For PhreeqcRM the return value of a function is generally used for error codes (negative values), although positive values may have other meanings (a positive return value for RM_GetChemistryCellCount is the number of chemistry cells).

Here is an extract of the declaration of RM_GetConcentrations found in RM_interface.F90 showing that the second argument is defined as INTENT(out), while the first argument is INTENT(in):

Code: [Select]

INTEGER FUNCTION RM_GetConcentrations(id, c)
...
    INTEGER, INTENT(in) :: id
    DOUBLE PRECISION, INTENT(out), DIMENSION(:,:) :: c


[Zhaoyang]

Yes, I have seen that in RM_interface.F90 and hence deleted the last post. Thanks very much for your help.

[Zhaoyang]

In the set initial conditions of the advection_f90.F90, the number of components for transport are printed in string1, while the components(i) and gfw(i) are printed in string. I have tried to change string to string1, but the results are incorrect. Why should the components(i) and gfw(i) be printed in string, instead of string1?

In addition, since I focus on material concentrations like NO3- and NH4+, should I use RM_GetSpeciesCount and RM_GetSpeciesName, instead of RM_FindComponent? In other words, should I use species rather than components?

Thanks tons for your explainations.

Code: [Select]

write(string1, "(A,I10)") "Number of components for transport:               ", RM_GetComponentCount(id)
  status = RM_OutputMessage(id, trim(string1))
  ! Get component information
  allocate(components(ncomps))
  allocate(gfw(ncomps))
  status = RM_GetGfw(id, gfw) ! Returns the gram formula weights (g/mol) for the components
  do i = 1, ncomps  ! ncomps, the number of components
     status = RM_GetComponent(id, i, components(i)) ! Retrieve an item from the reaction-module component
     write(string,"(A10, F15.4)") components(i), gfw(i)
     status = RM_OutputMessage(id, string)
  enddo

[dlparkhurst]

I'm not going to debug your changes. The way it is written in advection_f90.F90 works correctly. I think you should be able to replace string with string1 (both instances) in the do loop, but it is up to you to debug if a problem arises.

If you are not transporting with multi-component diffusion, I do not think you need to transport species. Transporting components (which assumes an advection/dispersion-dominated system) is sufficient to uniquely define the species concentrations in each cell.