PhreeqcRM: use RM_RunString with OpenMP

[Christinali91]

Hi David,

I try to use PhreeqcRM parallel capabilities in a two-phase reactive flow problem. I used EQUILIBRIUM_PHASES_MODIFY to update the gas phase(CO2). Everything works for the sequential run, but the results get weird when I use more than one thread.
As shown in the attached test case picture, I injected CO2 into pure water from the left corner. The color indicates pH value. In the sequential run, everything is perfectly propagating from the left corner to the upright corner, but in parallel, a weird value appears in the low right corner. I repeat the same simulation several times and found this type of weird value appears at random locations and random time steps.
I test the code without PhreeqcRM and without EQUILIBRIUM_PHASES_MODIFY. This type of issue doesn't occur in either of them so I have a suspicion on RM_RunString, but couldn't figure out how to fix it.
Below is the code I have:

Code: [Select]

        DO J = 1, NR_GC_NODES
                WRITE(DUMMY_STR,*) K-1
                input_IHC_temp = 'EQUILIBRIUM_PHASES_MODIFY  '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
                HC_GC_VALUES(K,1) = P(K)*1.0e-5
                !
                WRITE(DUMMY_STR,*)   log10(HC_GC_VALUES(K,1))       ! log10(pressure) 
                WRITE(DUMMY_STR_2,*)  HC_GC_VALUES(K,2)  !  TOT_CO2
                input_IHC_temp = input_IHC_temp //'-component  '//TRIM(ADJUSTL(NAME_HC_GC(1)))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-moles  ' //TRIM(ADJUSTL(DUMMY_STR_2))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-si  '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
                input_IHC_temp = input_IHC_temp // '-si_org '//TRIM(ADJUSTL(DUMMY_STR))//new_line(cc)
                input_IHC_temp = input_IHC_temp //'END' //new_line(cc)
                Ierr = RM_RunString(Id_phreeqc, 1, 0, 0, input_IHC_temp)
            ENDIF         
        ENDDO

       Ierr=RM_RunCells(Id_phreeqc)
I also tried put all the EQUILIBRIUM_PHASES_MODIFY # sections in one string array and do Ierr = RM_RunString(Id_phreeqc, 1, 0, 0, input_IHC) outside the loop, but it doesn't fix the issue. Could you give me some suggestions?

Thanks so much for your help and time!

Sincerely,
Christina

[dlparkhurst]

Try setting SetRebalanceFraction(0.0). This forces each cell to be processed by the same thread throughout the simulation. I don't think it should make a difference, but try it.

[Christinali91]

Hi David,

I tried SetRebalanceFraction(0.0). It fixed the problem! I can see it is probably not the most efficient way to parallel it, but I'll do the scaling analysis on my code and keep you updated.

Thanks a lot!

Best,
Christina

[dlparkhurst]

I'm not sure why disabling the rebalancing would fix the problem. The modify blocks should be run by each thread, and a block should be ignored if the block applies to a cell that is not in that thread.

Another approach is to modify only the equilibrium phases in each thread. Some pseudo code is given below that goes through each thread and modifies only the equilibrium phases in a thread. It would be most efficient to make one string for each thread that does not contain END statements. The loop over threads could be parallelized with OpenMP directives. Hopefully, rebalancing would also work.

Code: [Select]

  nthreads = RM_GetThreadCount(id)
  allocate(start_cell(n), end_cell(n))
  status = RM_GetStartCell(id, start_cell)
  status = RM_GetEndCell(id, start_cell)
  ! you can use OpenMP to parallize this loop
  do i = 1, nthreads
    thread_id = RM_GetIPhreeqcId(id, i-1)
do j = start_cell(i), end_cell(i)
  ! Accumulate equilibrium phases modify strings
enddo
RunString(thread_id, string)
  enddo