Using PhreeqcRM module with Fortran 90

[Christinali91]

Dear Dr. Parkhurst,

We are working on coupling PhreeqcRM to a reservior simulation code wrriten in Fortran 90. I've got most of the phreeqc functions working with the Fortran code but stuck on one problem related to keyword USER PUNCH. I've specify USER PUNCH in SELECTED OUTPUT section, but somehow the USER PUNCH components are not read by the program. When I print out the number of components using RM_GetSelectedOutputcolumncount(), this function only counts the number of molecular species, but doesn't count the number of user defined species using USER PUNCH keyword.

Below is the code I use to define USER PUNCH section.

Code: [Select]

    input = "SELECTED_OUTPUT 1" // new_line(c)
    input = input // "-RESET   FALSE" // new_line(c)
    input = input // "-high_precision" // new_line(c)
    input = input // "-user_punch" //new_line(c)
    input = input // 'USER_PUNCH  1'// new_line(c)

    DO I = 1, NELET
        input = input // '-Heading   '//STR_ELEMENT (I)// new_line(c)
    ENDDO

    DO I = 1, NPRINTG
        input = input // '-Heading   '//STR_PRINTG(I) // new_line(c)
    ENDDO

    DO  I=1, NO_HC_GC
        input = input //'-Heading   DELTA_'// TRIM(ADJUSTL(NAME_HC_GC(I))) // new_line(c)
        input = input //'-Heading    SI_'//TRIM(ADJUSTL(NAME_HC_GC(I))) // new_line(c)
    END DO
      !
      input = input // '-Heading   Charge' // new_line(c)
      input = input // '-Heading   Total_H' // new_line(c)
      input = input // '-Heading   Total_O' // new_line(c)

    DO  I=1,NELET
        CTR_INT = (I)*10
        WRITE(CTR_STR,*) CTR_INT
        input = input // CTR_STR//'  PUNCH TOTMOLE("'//TRIM(STR_ELEMENT (I))//'")' //new_line(c)
    END DO
    !! -------
    DO I=1, NPRINTG
        CTR_INT=(NELET+I)*10
        WRITE(CTR_STR,*) CTR_INT
        input = input // CTR_STR //'  PUNCH '//TRIM(STR_PRINTG(I)) //new_line(c)
    END DO
    !
    DO  I=1, NO_HC_GC
        CTR_INT=(NELET+NPRINTG+2*I-1)*10
        WRITE(CTR_STR,*) CTR_INT
        WRITE(DUMMY_STR,*) 1000.0d0
        input = input // CTR_STR //'  PUNCH EQUI_DELTA("'//TRIM(ADJUSTL(NAME_HC_GC(I)))//'")/'//TRIM(ADJUSTL(DUMMY_STR)) //new_line(c)
        CTR_INT=(NELET+NPRINTG+2*I)*10
        WRITE(CTR_STR,*) CTR_INT
        input = input // CTR_STR//'  PUNCH SI("'//TRIM(ADJUSTL(NAME_HC_GC(I)))//'")'//new_line(c)
    END DO
    !! ---- CHARGE BALANCE
    STARTING_EXTRA = 0
    ctr_int=(NELET+NPRINTG+2*NO_HC_GC+ STARTING_EXTRA+1)* 10
    WRITE(CTR_STR,*) CTR_INT
    input = input //CTR_STR//'  PUNCH CHARGE_BALANCE' //new_line(c)
    !
    ! C_total_H ----
    ctr_int=(NELET+NPRINTG+2*NO_HC_GC+STARTING_EXTRA+2)* 10
    WRITE(CTR_STR,*) CTR_INT
    input = input //CTR_STR//'  PUNCH TOTMOLE("H")' //new_line(c)  ! OR TOT
    !
    !C_total_O ----
    ctr_int=(NELET+NPRINTG+2*NO_HC_GC+ STARTING_EXTRA+3)* 10
    WRITE(CTR_STR,*) CTR_INT
    input = input //CTR_STR//'  PUNCH TOTMOLE("O")' //new_line(c)  ! OR TOT
    input = input // 'END' //new_line(c)
    input = input // 'KNOBS' //new_line(c)
    input = input // '-convergence_tolerance   1E-8' //new_line(c)


    Ierr = RM_RunString(ID_PHREEQC, 1, 1, 1, input)
    !write(*,'(a)') input

This is the output from the code:

Code: [Select]

SELECTED_OUTPUT 1
-RESET   FALSE
-high_precision
-user_punch
USER_PUNCH  1
-Heading   Ca
-Heading   Cl
-Heading   K
-Heading   N
-Heading   Na
-Heading   TOT("O")
-Heading   TOT("H")
-Heading   TOT("K")
-Heading   TOT("Na")
-Heading   -LA("H+")
-Heading   MOL("CaOH+")
-Heading   MOL("NaOH")
-Heading   MOL("KOH")
-Heading   MOL("H2O")
-Heading   Charge
-Heading   Total_H
-Heading   Total_O
          10                      PUNCH TOTMOLE("Ca")
          20                      PUNCH TOTMOLE("Cl")
          30                      PUNCH TOTMOLE("K")
          40                      PUNCH TOTMOLE("N")
          50                      PUNCH TOTMOLE("Na")
          60                      PUNCH TOT("O")
          70                      PUNCH TOT("H")
          80                      PUNCH TOT("K")
          90                      PUNCH TOT("Na")
         100                      PUNCH -LA("H+")
         110                      PUNCH MOL("CaOH+")
         120                      PUNCH MOL("NaOH")
         130                      PUNCH MOL("KOH")
         140                      PUNCH MOL("H2O")
         150                      PUNCH CHARGE_BALANCE
         160                      PUNCH TOTMOLE("H")
         170                      PUNCH TOTMOLE("O")
END
KNOBS
-convergence_tolerance   1E-8

When I do RM_GetSelectedOutputcolumncount() or RM_GetSelectedOutput(), I could only get 8 species instead of 17. Is there anything else I can modify to transport all 17 species in my code?

Thanks so much for your help!

[dlparkhurst]

The code that was generated looks to be okay.

You need to make sure the definition is run by the workers with a value of 1 for workers:

integer function RM_RunString   (   integer, intent(in)    id,
integer, intent(in)    initial_phreeqc,
integer, intent(in)    workers,
integer, intent(in)    utility,
character(len=*), intent(in)    input_string
)   

Then you need to do a little debugging. All of your values are defined in USER_PUNCH, so do you get different values if you PUNCH different values for the first 8 values? If you select some values with SELECTED_OUTPUT (like -mol H+), do they appear?

Other than a Fortran dimension error in your code, I don't have any other ideas.

[Christinali91]

Dear Dr. Parkhurst,

Thanks so much for the reply. I resolved the issue above by replacing RM_RunCells() with RM_RunString(). I am still a little bit confused about how to use these two functions? What's the main difference between RM_RunCells() and RM_RunString(). Is RM_RunString an equivalent version of AccumulateLines() in Iphreeqc?

My understanding is that we can define Selected_ouput using RM_RunString() function. I wrote a subroutine to generate selected_output block for phreeqc and then use RM_RunString to run it. However, I met an error message as below:

Code: [Select]

ERROR: Invalid argument.
ERROR: Selected output not found.
ERROR: Invalid argument.
ERROR: PhreeqcRM::GetSelectedOutputColumnCount

I printed and double-checked my selected_output block as below, but couldn't find where the invalid argument is. Could you give me some hints? I've already specified

Code: [Select]

Ierr = RM_SetSelectedOutputOn(Id_phreeqc_injection, 1) in the very beginning.

The printed out SELECTED_OUTPUT block(The code generated this block is omitted)

Code: [Select]

SELECTED_OUTPUT 1
-RESET FALSE
-high_precision
USER_PUNCH  1
-Heading   Ca
-Heading   Cl
-Heading   K
-Heading   N
-Heading   Na
-Heading   Charge
-Heading   Total_H
-Heading   Total_O
          10                     PUNCH TOTMOLE("Ca")
          20                     PUNCH TOTMOLE("Cl")
          30                     PUNCH TOTMOLE("K")
          40                     PUNCH TOTMOLE("N")
          50                     PUNCH TOTMOLE("Na")
          60                      PUNCH CHARGE_BALANCE
          70                      PUNCH TOTMOLE("H")
          80                      PUNCH TOTMOLE("O")
END
KNOBS
-convergence_tolerance   1E-8
RUN_CELLS; -CELLS            0                    ;END

The code where gives me the error:

Code: [Select]

Ierr=RM_SetCurrentSelectedOutputUserNumber(Id_phreeqc_injection,1)

WRITE(DUMMY_STR,*) K
input3 = input3 //'RUN_CELLS; -CELLS '//DUMMY_STR// ';END'//new_line(cc)

Ierr = RM_RunString(Id_phreeqc_injection,1, 1, 1, input3)

! Ierr = RM_RunCells(Id_phreeqc_injection)

N_COLUMN = RM_GetSelectedOutputColumnCount(Id_phreeqc_injection)
N_ROW = RM_GetSelectedOutputRowCount(Id_phreeqc_injection)


Thanks so much for your help and time!

[dlparkhurst]

The strategy for using PhreeqcRM is first to set up the initial conditions and selected output that you want, and then to sequentially transport the elements, run the chemistry, and process selected_output results.

So PhreeqcRM has three IPhreeqc instances, one to process initial conditions (initial phreeqc), one to store and run the contents of each cell (worker), and one we won't worry about. The intent is to use the initial phreeqc instance to read a file (RM_RunFile) or process a string (RM_RunString), usually with only the initial phreeqc instance, so that a number of solutions and chemical reactants (EQUILIBRIUM_PHASES, SURFACE, EXCHANGE, etc) are defined to the initial phreeqc instance. These solutions and reactants are identified by integer numbers SOLUTION 1, SOLUTION 10, etc. It is then possible to populate the the worker IPhreeqc instance(s) with the solutions and reactants you want in each cell of the transport model by using InitialPhreeqc2Module. Basically, you say I want solution 1, equilibrium_phases 2, and exchange 1 in cell number 1 of your reaction transport model, and so on for every cell. At this point, all the cells of the transport model have initial conditions, and usually, you are done using the initial IPhreeqc instance (except perhaps for boundary conditions). You would use RM_GetConcentrations to retrieve the concentrations of the cells (from the worker) to set up your first transport step.

Now, to be able to access selected output result for each cell, it is necessary to define SELECTED_OUTPUT and USER_PUNCH to the worker instances before you start running the chemistry (RM_RunCells). To do this, you need to use RM_RunString or RM_RunFile with one or more SELECTED_OUTPUT/USER_PUNCH definitions (you can define multiple selected output with different user numbers). You have to be sure that you set the option to run the string/file (that has the selected output definitions) in the worker instance.

At this point you should be able to transport, and then use RM_SetConcentrations (to put the tranported concentrations in the worker) and then RM_RunCells. The initial phreeqc instance is not used in RM_RunCells, only the worker, which contains all the information for each cell in the transport model. Essentially, RM_RunCells internally performs a PHREEQC RunCells for each cell in the transport model. When each cell is run, the current solution for a cell reacts with the other reactants in the cell, and the new results are written for selected output and the new solution and reactant compositions are saved for the cell (in the worker). After all the cells are run with RM_RunCells, you can use the various RM_ methods to retrieve the selected output results.

Please look at the Fortran advection example that is included with PhreeqcRM. It uses most of the PhreeqcRM methods, and has the logical sequence for their use.

[Christinali91]

Dear Dr. Parkhurst,

Thanks so much for the quick reply. Your explanation really gives me a better idea of the workflow.

We used to couple our transport simulator with IPhreeqc module. In the IPhreeqc module, we can use accumulatelines() function with SOLUTION_MODIFY block within the time iteration to modify the input file in the computer memory while maintaining all relevant information from an initialization step. How can I achieve this using the PhreeqcRM module? Do RM_SetConcentrations and RM_GetConcentrations take care of this?

Thanks a lot!

[dlparkhurst]

Yes, the intention is to use RM_GetConcentrations to retrieve total element concentrations (including H, O, and charge balance) from the module, transport the elements, send the transported concentrations to the module with RM_SetConcentrations, and then run reactions with RM_RunCells.

(No Dr. please.)

[Christinali91]

Hi David,

I've been studying the FORTRAN 90 example for PhreeqcRM. It mentioned that we could use the utility instance of PhreeqcRM for additional features. I am not sure if I understand this instance correctly. For example, I want to inject the CaCl solution from one cell into the domain. Should I use instance utility for it since the components, as well as SOLUTION, are different in this cell? Or should I just specify the SOLUTION for it using RM_InitialPhreeqc2Module with the same ID?

Thanks a lot!

Best,

[dlparkhurst]

No, I doubt you need to use the utility instance. The utility instance was added to do PHREEQC calculations not directly related to transport. Honestly, I have only used it once to be able to calculate the pH of water extracted from a well that was a mixture of water from multiple cells. I had the concentrations of the elements from the well, but I did not know the pH, so I used the utility instance to calculate the pH. (The pH is known for cells, so in most cases, this calculation is not necessary.)

It sounds like you are talking about boundary conditions. The intent was to define solutions for the boundary conditions in the initial phreeqc instance and then retrieve them with RM_InitialPhreeqc2Concentrations. This would give you sets of concentrations that you would use to define the transport system for each element/component. For instance, the concentration of SOLUTION 10 in the initial phreeqc instance could be retrieved and used to define the concentrations used in the transport equations for each element entering the domain through, for example, a flux boundary.

[Christinali91]

Hi David,

Thanks so much for helping me out! I am able to get the code running now. However, there is still one problem left I couldn't figure out.
I am able to get the correct concentration value using RM_GetConcentrations(), but the result from RM_GetSelectedOutput() is not correct as shown below:

Code: [Select]

GLOB_TIME           0
 Cell number            1
      Concentrations:
           1           H:   110.6821
           2           O:    55.3416
           3      Charge:    -0.0000
           4          Ca:     0.0006
           5          Cl:     0.0012
           6           K:     0.0000
           7           N:     0.0000
           8          Na:     0.0000
      Selected output:
            1                    Ca:      5.9821E-05
            2                    Cl:      1.1964E-04
            3                     K:      0.0000E+00
            4                     N:      0.0000E+00
            5                    Na:      0.0000E+00
            6                Charge:     -7.2835E-13
            7               Total_H:      1.1068E+01
            8               Total_O:      5.5342E+00

The correct Selected output should be:
 GLOB_TIME           0
 Cell number            1
      Selected output:
            1  :      6.0000E-04
            2  :      1.2000E-03
            3  :      0.0000E+00
            4  :      0.0000E+00
            5  :      0.0000E+00
            6  :      4.9738E-16
            7  :      1.1101E+02
            8  :      5.5507E+01

I've defined SELECTED_OUTPUT and USER_PUNCH to the worker instances before RM_RunCell() and specified the SELECTED_OUTPUT block I want to use. Is there anything else missing in my code?

Code: [Select]

    nlines = 50
    nlines = nlines + NELET
    nlines = nlines + NELET
    allocate (character(len=nlines * 100)::input3)

    input3 = 'SOLUTION_SPECIES;H2O + 0.01e- = H2O-0.01; log_k -9' // new_line(cc)
    input3 = input3 // "SELECTED_OUTPUT 1" // new_line(cc)
    input3 = input3 // "-RESET FALSE" // new_line(cc)
    input3 = input3 // "-high_precision" // new_line(cc)
    input3 = input3 // 'USER_PUNCH  1'// new_line(cc)

    DO I = 1, NELET
        input3 = input3 // '-Heading   '//STR_ELEMENT (I)// new_line(cc)
    ENDDO

      !
    input3 = input3 // '-Heading   Charge' // new_line(cc)
    input3 = input3 // '-Heading   Total_H' // new_line(cc)
    input3 = input3 // '-Heading   Total_O' // new_line(cc)

    starting_extra=0
    DO  I=1, NELET
        ctr_int=(starting_extra+I)* 10
        WRITE(CTR_STR,*) CTR_INT
        input3 = input3 //CTR_STR//' PUNCH TOTMOLE("'//TRIM(STR_ELEMENT (I))//'")'// new_line(cc)
    END DO

    ctr_int=(NELET+starting_extra+1)* 10
    WRITE(CTR_STR,*) CTR_INT
    input3 = input3 //CTR_STR //'  PUNCH CHARGE_BALANCE' // new_line(cc)

    !C_total_H ----
    ctr_int=(NELET+starting_extra+2)* 10
    WRITE(CTR_STR,*) CTR_INT
    input3 = input3 //CTR_STR //'  PUNCH TOTMOLE("H")'// new_line(cc)  ! OR TOT

    !C_total_O ----
    ctr_int=(NELET+starting_extra+3)* 10
    WRITE(CTR_STR,*) CTR_INT
    input3 = input3 //CTR_STR //'  PUNCH TOTMOLE("O")'// new_line(cc)  ! OR TOT
    input3 = input3 // 'END' // new_line(cc)
    input3 = input3 // 'KNOBS' // new_line(cc)
    input3 = input3 // '-convergence_tolerance   1E-8'// new_line(cc)  !! SET TO DEFAULT

    Ierr = RM_RunString(Id_phreeqc_injection, 1, 1, 0, input3)


Code: [Select]

       Ierr=RM_SetCurrentSelectedOutputUserNumber(Id_phreeqc_injection, 1)
       Ierr = RM_RunCells(Id_phreeqc_injection)

       Ierr = RM_GetConcentrations(Id_phreeqc_injection, concent_injet)

       N_COLUMN = RM_GetSelectedOutputColumnCount(Id_phreeqc_injection)
       N_ROW = RM_GetSelectedOutputRowCount(Id_phreeqc_injection)

       ALLOCATE(VECTOR_IPHREEQC_TEMP(1:N_ROW,1:N_COLUMN))
       VECTOR_IPHREEQC_TEMP = 0.0D0
       Ierr = RM_GetSelectedOutput(Id_phreeqc_injection, VECTOR_IPHREEQC_TEMP)

Thanks again for your time and help!

Best,

[dlparkhurst]

I think you are looking at different units.

The units of the values of RM_GetConcentrations are determined by RM_SetUnitsSolution, and can be mg/L, mol/L, or mass fraction.

The units you are printing in SELECTED_OUTPUT are moles. Now, the representative volume for the cell is 1 L by default, and it appears you are using a porosity of 0.1, so you have approximately 0.1 liter of solution in a cell. The number of moles of O in 0.1 L of H2O is about 5.5, and the number of moles of all elements will be the number of moles in 0.1 L of solution.

Note that the conversion from moles per kg water in PHREEQC to mg/L or mol/L involves the solution volume (and solution density), whereas mass fraction is directly calculable from the number of moles of elements. If you want to compare values of mg/L, for Cl, you need to print TOTMOL("Cl")*GFW("Cl")*1000/SOLN_VOL in USER_PUNCH.

[Christinali91]

Hi David,

Thanks a lot for the explanation. My code is running correctly now.
During the debugging, I found that there might be a typo in the PhreeqcRM manual for Fortran.
In the function definition of RM_RunString(), the first instance indicates initial_phreeqc and the second one indicates workers as below:

Code: [Select]

integer function RM_RunString ( integer, intent(in) id,
integer, intent(in) initial_phreeqc,
integer, intent(in) workers,
integer, intent(in) utility,
character(len=*), intent(in) input_string
)
While in the Fortran example for phreeqcRM, the sequence is inverse. The first instance is noted as workers as below:

Code: [Select]

Fortran Example:
 string = "DELETE; -all"
 status = RM_RunString(id, 1, 0, 1, string)  ! workers, initial_phreeqc, utility
I thought this may cause some confusion for others as well.

I also have a question about openMP. When I use multithread in phreeqcRM, I get a lot of print out messages every timestep as below:

Code: [Select]

          Estimated efficiency of chemistry without communication: 99.1209
          Cells shifted between threads     0
          Time rebalancing load             2.7895e-05
          Estimated efficiency of chemistry without communication: 99.2936
          Cells shifted between threads     0
          Time rebalancing load             3.8147e-05
          Estimated efficiency of chemistry without communication: 99.0095
          Cells shifted between threads     0
          Time rebalancing load             2.81334e-05
          Estimated efficiency of chemistry without communication: 98.7377
          Cells shifted between threads     0
          Time rebalancing load             3.00407e-05
          Estimated efficiency of chemistry without communication: 98.781
          Cells shifted between threads     0
          Time rebalancing load             2.69413e-05
          Estimated efficiency of chemistry without communication: 98.7068
          Cells shifted between threads     0
          Time rebalancing load             3.8147e-05
          Estimated efficiency of chemistry without communication: 99.0646
          Cells shifted between threads     0
          Time rebalancing load             4.60148e-05
Is there a way to turn it off? These print-outs really slow down my code in longer runs.

Thanks again for your help!         

[dlparkhurst]

RM_SetScreenOn.

Yikes, the order of the flags is actually a bug in the code. Regardless of documentation or comments, the order for RM_RunString is workers, initial_phreeqc, utility.

[Christinali91]

Hi David,

Is there any setting we need to do to enable OpenMP? In my current code, the efficiency of OpenMP is much lower than the results shown in your paper.  For 10 processors the goal is a 10-fold speed-up, or 1000%, but the biggest speed-up I got now is 40% for 10 processors, or 4% for each processor.

I doubt the OpenMP is actually enabled, but I check the number of threads it uses with 'top' command. It shows that this job is indeed using 10 threads.

After assigning nthreads = 10, Id_phreeqc = RM_Create(NR_GC_NODES, nthreads). I also add the openmp flag in the make file as below:

Code: [Select]

OPENMPPHREEQC = -qopenmp -I$(MKLROOT)/include
$(FC) $(FAST) $(OPENMPPHREEQC) $? -c  -module $(MOD) -o $@

Is there anything else I need to do to enable OpenMP run? Do I need something like !$OMP in the code? 

Thanks so much for your help and time!

[dlparkhurst]

Assuming you used optimization and defined the compiler directive USE_OPENMP when you compiled the library, seems like you should be using OpenMP correctly. From some of your output, I think this must be true. It would not try to rebalance the load unless OpenMP (or MPI) were invoked.

All I can think is that the chemistry calculation is not a large fraction of the calculation load. Either the tranpsort far outweighs the chemistry, or there are not enough cells or a difficult enough calculation to get meaningful timings.

You may need to add some timing in your code to evaluate where time is spent. The function omp_get_wtime() returns a double, and the difference between two calls is a number of seconds. You could also time an unrelated loop to test the effectiveness of OpenMP with different numbers of processors, something like

  !$OMP PARALLEL DO
  do i = 1, n
     dummy(i) = sqrt(dble(i))
  end do
  !$OMP END PARALLEL DO

where n is a large number.

[Christinali91]

Hi David,

I tested other OpenMP code on our server. It turns out that the efficiency problem may come from the setting of the server, so I plan to test the code on Mac as well and compare the  performance.

During the compiletion on Mac, I met some compiling errors. Since the compiletion topic is not belong to this section, I'll open another discussion under installation.

Thanks so much for your help and timely reply!

Best,

[dlparkhurst]

I would be surprised if the server did not make use of OpenMP. My guess is that you have not turned on some compiler flag. I attach a program that works on my computer. Using this computer, the time to run the test with PhreeqcRM decreases from about 4 seconds to 1 second as processors increase from 1 to 8. Similarly run times decrease from 4 to 1 second for OpenMP running a timing loop without PhreeqcRM.

For PhreeqcRM you need make sure your Makefile has RM_interface.F90 added and the OpenMP version of the PhreeqcRM library is linked. You also need phreeqc.dat, currently defined as ../phreeqc.dat in the code, but you can change that to wherever phreeqc.dat is located. For PhreeqcRM run times increase with the number of cells (nxyz) and with the index in the loop in USER_PUNCH (emulating more complicated chemistry).

Can't help you with Mac compilation.

[Christinali91]

Dear David,

Thanks so much for the test code. It helps me a lot while digging into the scaling issue of our code. We finally figured out the reason of the unlinear scaling after coupling PhreeqcRM to the reservior simulation code. I ignored the change of CPU frequency while doing the OpenMP test, which had a hugh impact on the running speed. After fixed that, PhreeqcRM performs almost perfect in the test run.

Thanks again for all your time and help!

Best,