# PhreeqcUsers Discussion Forum

## Processes => Reactive transport modelling => Topic started by: moein on May 07, 2019, 02:40:19 PM

Title: Calling PhreeqcRM shared library (.so) from Python
Post by: moein on May 07, 2019, 02:40:19 PM
I have been trying to call PhreeqcRM functions from Python and I have been successful for most of the functions except one, here you can see how I do it (and Im just putting the part of the code related to that function):

from ctypes import *
import numpy as np
libc = CDLL('libphreeqcrm-3.5.0.so')

ic1 = np.zeros(nxyz*7,dtype=np.int32).reshape(nxyz, 7)
ic2 = np.zeros(nxyz*7,dtype=np.int32).reshape(nxyz, 7)
f1 = np.zeros(nxyz*7,dtype=np.float64).reshape(nxyz, 7)

ic1[:,0] = 1       # Solution 1
ic1[:,1] = -1      # Equilibrium phases none
ic1[:,2] = 1       # Exchange 1
ic1[:,3]= -1      # Surface none
ic1[:,4]= -1      # Gas phase none
ic1[:,5]= -1      # Solid solutions none
ic1[:,6]= -1      # Kinetics none
ic2[:,0]= -1      # Solution none
ic2[:,1]= -1      # Equilibrium phases none
ic2[:,2]= -1      # Exchange none
ic2[:,3]= -1      # Surface none
ic2[:,4]= -1      # Gas phase none
ic2[:,5]= -1      # Solid solutions none
ic2[:,6]= -1      # Kinetics none
f1[:,0]= 1.0      # Mixing fraction ic1 Solution
f1[:,1]= 1.0      # Mixing fraction ic1 Equilibrium phases
f1[:,2]= 1.0      # Mixing fraction ic1 Exchange 1
f1[:,3]= 1.0      # Mixing fraction ic1 Surface
f1[:,4]= 1.0      # Mixing fraction ic1 Gas phase
f1[:,5]= 1.0      # Mixing fraction ic1 Solid solutions
f1[:,6]= 1.0      # Mixing fraction ic1 Kinetics

status = libc.RM_InitialPhreeqc2Module(id,ic1.flatten('F').ctypes,ic2.flatten('F').ctypes,f1.flatten('F').ctypes)

here is what I get:

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

ERROR: PhreeqcRM failed.
ERROR: Processing initial conditions.
ERROR: InitialPhreeqc2Module: Failure in PhreeqcRM

ERROR: PhreeqcRM failed.
ERROR: PhreeqcRM::InitialPhreeqc2Module
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

it seems that the ic1 array, when its transported to the function, instead of -1, becomes 1, I was wondering if anyone could help me to resolve the issue. in the following you can find the phreeqc input file:

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
TITLE Example 11.--Transport and cation exchange.
SOLUTION 0  CaCl2
units            mmol/kgw
temp             25.0
pH               7.0     charge
pe               12.5    O2(g)   -0.68
Ca               0.6
Cl               1.2
SOLUTION 1-40  Initial solution for column
units            mmol/kgw
temp             25.0
pH               7.0     charge
pe               12.5    O2(g)   -0.68
Na               1.0
K                0.2
N(5)             1.2
END
EXCHANGE 1-40
-equilibrate 1
X                0.0011
END
SELECTED_OUTPUT
-reset           false
#        -step
-totals          Na Cl K Ca
USER_PUNCH