# Examples for PhreeqPy¶

This is an example for the use PhreeqPy for one-dimensional advection. It is the example 11 from the PHREEQC users manual.

```"""Advection with DLL or COM server.

Using MODFIY we update the concentration on every
time step. We shift by one cell per step.
"""

from __future__ import print_function

import sys

# Simple Python 3 compatiblity adjustment.
if sys.version_info == 2:
range = xrange

import os
import timeit

MODE = 'dll'  # 'dll' or 'com'

if MODE == 'com':
import phreeqpy.iphreeqc.phreeqc_com as phreeqc_mod
elif MODE == 'dll':
import phreeqpy.iphreeqc.phreeqc_dll as phreeqc_mod
else:
raise Exception('Mode "%s" is not defined use "com" or "dll".' % MODE)

def make_initial_conditions():
"""
Specify initial conditions data blocks.

Uniform initial conditions are assumed.
"""
initial_conditions = """
TITLE Example 11.--Transport and ion 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  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
equilibrate 1
X                0.0011
END
"""
return initial_conditions

def make_selected_output(components):
"""
Build SELECTED_OUTPUT data block
"""
for i in range(len(components)):
selected_output = """
SELECTED_OUTPUT
-reset false
USER_PUNCH
"""
#
# charge balance, H, and O
#
code = '10 w = TOT("water")\n'
code += '20 PUNCH CHARGE_BALANCE, TOTMOLE("H"), TOTMOLE("O")\n'
#
# All other elements
#
lino = 30
for component in components:
code += '%d PUNCH w*TOT(\"%s\")\n' % (lino, component)
lino += 10
selected_output += code
return selected_output

def initialize(cells, first=False):
"""
Initialize IPhreeqc module
"""
phreeqc = phreeqc_mod.IPhreeqc()
initial_conditions = make_initial_conditions()
phreeqc.run_string(initial_conditions)
components = phreeqc.get_component_list()
selected_output = make_selected_output(components)
phreeqc.run_string(selected_output)
phc_string = "RUN_CELLS; -cells 0-1\n"
phreeqc.run_string(phc_string)
conc = get_selected_output(phreeqc)
inflow = {}
initial = {}
for name in conc:
if first:
inflow[name] = conc[name]
else:
inflow[name] = conc[name]
initial[name] = conc[name]
task += "COPY solution 1 %d-%d\n" % (cells, cells)
task += "COPY exchange 1 %d-%d\n" % (cells, cells)
task += "RUN_CELLS; -cells %d-%d\n" % (cells, cells)
conc = get_selected_output(phreeqc)
for name in conc:
value = [initial[name]] * len(conc[name])
conc[name] = value
return phreeqc, inflow, conc

"""Advect by shifting concentrations from previous time step.
"""
all_names = conc.keys()
names = [name for name in all_names if name not in ('cb', 'H', 'O')]
for name in conc:
# shift one cell
conc[name][1:] = conc[name][:-1]
conc[name] = inflow[name]
modify = []
for index, cell in enumerate(range(cells, cells + 1)):
modify.append("SOLUTION_MODIFY %d" % cell)
modify.append("\t-cb      %e" % conc['cb'][index])
modify.append("\t-total_h %f" % conc['H'][index])
modify.append("\t-total_o %f" % conc['O'][index])
modify.append("\t-totals")
for name in names:
modify.append("\t\t%s\t%f" % (name, conc[name][index]))
modify.append("RUN_CELLS; -cells %d-%d\n" % (cells, cells))
cmd = '\n'.join(modify)
phreeqc.run_string(cmd)
conc = get_selected_output(phreeqc)
return conc

def get_selected_output(phreeqc):
"""Return calculation result as dict.

Header entries are the keys and the columns
are the values as lists of numbers.
"""
output = phreeqc.get_selected_output_array()
conc = {}
for row in output[1:]:
return conc

def run(ncells, shifts, specie_names):
"""Do one run in one process.
"""
cells = (1, ncells)
phreeqc, inflow, conc = initialize(cells, first=True)
outflow = {}
for name in specie_names:
outflow[name] = []
for _counter in range(shifts):
conc = advect_step(phreeqc, inflow, conc, cells)
for name in specie_names:
outflow[name].append(conc[name][-1])
return outflow

def write_outflow(file_name, outflow):
"""Write the outflow values to a file.
"""
fobj = open(file_name, 'w')
fobj.write('\n')
fobj.write('\n')

def main(ncells, shifts):
"""Run different versions with and without multiprocessing
"""

def measure_time(func, *args, **kwargs):
"""Convinience function to measure run times.
"""
start = timeit.default_timer()
result = func(*args, **kwargs)
return result, timeit.default_timer() - start

print('Dimensions')
print('==========')
print('number of cells:   ', ncells)
print('number of shifts   ', shifts)
specie_names = ('Ca', 'Cl', 'K', 'N', 'Na')
outflow, run_time = measure_time(run, ncells, shifts, specie_names)
if not os.path.exists('data'):
os.mkdir('data')
write_outflow('data/out.txt', outflow)
print('run time:', run_time)
print("Finished simulation\n")

if __name__ == '__main__':
main(ncells=40, shifts=120)
```