paleolimbot / tidyphreeqc

Tidy geochemical modeling using PHREEQC
22 stars 3 forks source link

Fehler in phr_input_section #18

Closed loewche closed 1 year ago

loewche commented 1 year ago

Dear sir or madame,

I am trying to use your package tidyphreeqc. However, I run in the following issue:

This example works! There is no error message. This example runs smoothly:

phr_selected_output( activities = c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "Al+3"), saturation_indices = c("Al(OH)3(a)", "Gibbsite"), molalities = c("Al(OH)3(a)", "Gibbsite"), pe = TRUE, alkalinity = TRUE)

This is another example I want to wish examined which won`t work:

phr_selected_output( activities = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-"), saturation_indices = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-"), molalities = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-") )

Here I get the following error message:

Error in phr_input_section("SELECTED_OUTPUT", number = .number, name = "", : 'number' must be a vector of length 1

What does the error mean and how I can fix this? I already note the lines here:

https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-56.html

To resolve this error I tried several solutions, nothing won't work.

A clue would be nice what a number of components you mean?

phr_selected_output(.number = 1, ...)

What's the difference between them both? Why the first lines are running and the lines from the second example won't?

What number of component? The numbers in the selected output, the number of species, the number of the components in the input? Which number is required to fulfill the expectations?

Sincerely yours

Ignimbrit commented 1 year ago

Hi loewche,

I was trying to understand your problem but am having trouble replicating your error. Both selected output definitions work just fine for me. Here is what I tried:

library(tidyphreeqc)

# So we define a solution, the hydrochemical state of which we want to model
pqi_solution <- phr_solution(
  pH = 7, temp = 25, pe = 4.5, redox = "O(-2)/O(0)", units = "ppm", 
  density = 1.02, Ca = 80, 'S(6)' = "96.     as SO4", 
  'S(-2)' = "1.      as S", `N(5) N(3)` = "14.     as N",
  `O(0)` = 8.0, 
  C = "61.0    as HCO3      CO2(g)     -3.5",
  Fe = "55.     ug/kgs as Fe S(6)/S(-2) Pyrite",
  '-water' = "0.5     # kg"
)

# Now we need to instruct PHREEQC, what information from the model run
# it should actually record for us in the output.
# To mimic the problem, we try two different kinds of output recordings

# this one allegedly works
pqi_selout_a <- phr_selected_output(
  activities = c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "Al+3"),
  saturation_indices = c("Al(OH)3(a)", "Gibbsite"),
  molalities = c("Al(OH)3(a)", "Gibbsite"),
  pe = TRUE,
  alkalinity = TRUE)

# while this one for reasons yet unclear apparently does not
pqi_selout_b <- phr_selected_output(
  activities = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-"),
  saturation_indices = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-"),
  molalities = c("Ca+2", "CaSO4", "HCO3-", "CO3-2", "OH-", "HSO4-")
)

# Now we have three different input objects. We need to group them
# and send them to phreeqc.
# There can always ever be one selected output 
# (that's something of an oversimplification but please bear with me).
# So we need two model runs for the different SELECTED_OUTPUT blocks we
# created. Both are tested against the same SOLUTION specification

# Model input Run A
pqi_input_a = phr_input(pqi_solution, pqi_selout_a)

# Model input Run B
pqi_input_b = phr_input(pqi_solution, pqi_selout_b)

# Let's run the models and capture the output

result_a = phr_run(pqi_input_a)

result_b = phr_run(pqi_input_b)

# No error for me here. Lets check it result_b will print
print(result_b)

result_b prints:

<phr_run_output>
PHREEQC run with 1 selected output(s)
Raw output at '/tmp/Rtmp9tJVEd/filef211b992135'
as_tibble():
# A tibble: 1 × 27
  selected_output   sim state   soln dist_x  time  step    pH    pe `m_Ca+2(mol/kgw)` `m_CaSO4(mol/kgw)` `m_HCO3-(mol/kgw)`
  <chr>           <int> <chr>  <int>  <dbl> <dbl> <int> <dbl> <dbl>             <dbl>              <dbl>              <dbl>
1 n1                  1 i_soln     1     NA    NA    NA     7   4.5           0.00185           0.000147          0.0000518
# ℹ 15 more variables: `m_CO3-2(mol/kgw)` <dbl>, `m_OH-(mol/kgw)` <dbl>, `m_HSO4-(mol/kgw)` <dbl>, `la_Ca+2` <dbl>, la_CaSO4 <dbl>,
#   `la_HCO3-` <dbl>, `la_CO3-2` <dbl>, `la_OH-` <dbl>, `la_HSO4-` <dbl>, `si_Ca+2` <dbl>, si_CaSO4 <dbl>, `si_HCO3-` <dbl>,
#   `si_CO3-2` <dbl>, `si_OH-` <dbl>, `si_HSO4-` <dbl>

all perfectly valid.

Could you give me a little more information about what you actually fed to phreeqc? Maybe provide us with a minimum (not) working example of what you would expect to work / would like to work, so that we can walk through it step by step and see where it actually breaks.

loewche commented 1 year ago

Hej,

first, I am glad that you got some time to look for my problem. And second I have to apologize, that I didn't make my problem more clear.

So second try:

Here is the first code: These lines are running perfectly:

tidyphreeqc::phr_run( tidyphreeqc::phr_solution_list(Al = 0.1, pe = 4, pH = seq(1, 13, by = 0.1)), phr_selected_output( activities = c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "Al+3"), saturation_indices = c("Al(OH)3(a)", "Gibbsite"), molalities = c("Al(OH)3(a)", "Gibbsite"), pe = TRUE, alkalinity = TRUE))

PHREEQC run with 1 selected output(s) Raw output at 'C:\Users\XXXX\AppData\Local\Temp\Rtmp8e6MZd\file197c30d2e3f' as_tibble(): # A tibble: 121 × 18 selected_output sim state soln dist_x time step pH pe `Alk(eq/kgw)` `m_Al(OH)3(a)(mol/kgw)` 1 n1 1 i_soln 1 NA NA NA 1 4 -0.118 0 2 n1 1 i_soln 2 NA NA NA 1.1 4 -0.0927 0 3 n1 1 i_soln 3 NA NA NA 1.2 4 -0.0729 0 4 n1 1 i_soln 4 NA NA NA 1.3 4 -0.0573 0 5 n1 1 i_soln 5 NA NA NA 1.4 4 -0.0450 0 6 n1 1 i_soln 6 NA NA NA 1.5 4 -0.0354 0 7 n1 1 i_soln 7 NA NA NA 1.6 4 -0.0279 0 8 n1 1 i_soln 8 NA NA NA 1.7 4 -0.0219 0 9 n1 1 i_soln 9 NA NA NA 1.8 4 -0.0173 0 10 n1 1 i_soln 10 NA NA NA 1.9 4 -0.0136 0 # ℹ 111 more rows # ℹ 7 more variables: `m_Gibbsite(mol/kgw)` , `la_Al(OH)4-` , `la_Al(OH)3` , # `la_Al(OH)2+` , `la_Al+3` , `si_Al(OH)3(a)` , si_Gibbsite # ℹ Use `print(n = ...)` to see more rows In another try, I thought to run the following lines: And I get this? tidyphreeqc::phr_run( tidyphreeqc::phr_solution(pH = 7, temp = 20), phr_equilibrium_phases(Gypsum = c(0)))%>% phr_print_output() phr_selected_output(.number = "Gypsum", molalities = "OH-", "H+", "Ca2+", "CasO4", "S(-2)") SELECTED_OUTPUT Gypsum - c(\"SOLUTION 1\", \" pH 7\", \" temp 20\", \"EQUILIBRIUM_PHASES 1\", \" Gypsum 0\") C:\Users\peter\AppData\Local\Temp\Rtmp8e6MZd\file197c285e143f list() NULL -molalities OH- - H+ - Ca2+ - CasO4 - S(-2) Why the output works fine in the first example, and why not in the second one? What is their difference? Admittedly, these are very simple examples. But I am trying to get back into phreeqC. The end goal is to use your package on one hand to create a web interface in Shiny for normal users. The second goal is to simulate the dissolution of Feldspar on a drill (5000 m) core to check the silicon content from the XRD or XRF analysis. Actually I started to work in phreeqC directly and pass the text files into the editor, after that I import the outcome back in R, which is a bit too much. I tried several packages, it seems that your package is more useful for me than the other ones. That's the reason I won't give up at this stage. Thanks again for going over!
Ignimbrit commented 1 year ago

Hi loewche,

it's not my library but maybe I can help anyway. There seems to be some structural problem with your mental model of either phreeqc or tidyphreeqc. That's perfectly understandable. I too needed quite some time to wrap my head around. Here is what I did to make your second, apparently dysfunctional code work. Please have a look if this is what you want (I still can't quite tell for sure but maybe we are getting closer :-))

# Let's do this very slowly, step by step, so that we won't miss anything

# First we have a solution
le_sol = phr_solution(pH = 7, temp = 20)

# Second, we add a solid, an equilibrium phase
le_mineral = phr_equilibrium_phases(Gypsum = c(0))
# This strictly means gypsum is saturated but no solid gypsum is present
# I suppose it's valid code but kinda can't see where this is going.

# Third, we want phreeqc to write down the output
le_selout = phr_selected_output(
  equilibrium_phases = "Gypsum",
  molalities = "OH-", "H+", "Ca2+", "CasO4", "S(-2)"
  )

# to wrap it up, we complie the phreeqc-program as one monolithic string.
pqi = phr_input(le_sol, le_mineral, le_selout)

# make phreeqc run it
pqo = phr_run(pqi)
print(tibble::as_tibble(pqo)|> _[,8:15])

To me this prints:

# A tibble: 2 × 8
     pH    pe `m_OH-(mol/kgw)` `m_Ca2+(mol/kgw)` `m_CasO4(mol/kgw)` `m_S(-2)(mol/kgw)` Gypsum d_Gypsum
  <dbl> <dbl>            <dbl>             <dbl>              <dbl>              <dbl>  <dbl>    <dbl>
1  7      4       0.0000000685                 0                  0                  0   0      0     
2  7.08  11.2     0.0000000996                 0                  0                  0   9.99  -0.0148```
loewche commented 1 year ago

Hej,

sorry, for misinterpreting the link to report a bug in the package. I simply followed in R the link BugReports: https://github.com/paleolimbot/tidyphreeqc/issues. And you answered, hmm. Anyway, I am happy that someone is trying to help and understand what I am ran into,

Foremost, I am a bit confused! Why you split the functions into three parts? Because the package offers a very smart possibility to do this in one (tidy) way.

tidyphreeqc::phr_run( tidyphreeqc::phr_solution(pH = 7, temp = 20), phr_equilibrium_phases(Gypsum = c(0)))%>%

phr_print_output()

phr_selected_output(.number = "Gypsum", molalities = "OH-", "H+", "Ca2+", "CasO4", "S(-2)")

And yet, there is an error. The output won't run. In opposite to your solution, the output runs. I don't really understand the difference here. In essence, they both do the same thing. But actually they do not.

THX,

again, to work with me to find a solution for this!