Closed defencedog closed 11 months ago
Furher findings, the error is loading .so
>>> from ctypes import *
>>> cdll.LoadLibrary('/data/data/com.termux/files/usr/local/lib/libiphreeqc-3.7.3.so')
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/data/data/com.termux/files/usr/lib/python3.9/ctypes/__init__.py", line 452, in LoadLibrary
return self._dlltype(name)
File "/data/data/com.termux/files/usr/lib/python3.9/ctypes/__init__.py", line 374, in __init__
self._handle = _dlopen(self._name, mode)
OSError: dlopen failed: cannot locate symbol "__trunctfdf2" referenced by "/data/data/com.termux/files/usr/local/lib/libiphreeqc-3.7.3.so"...
Solved loading issue, recompile with gcc-11 instead of clang
In termux pkg i gcc-default-11
above steps / log remains same...
./configure --prefix=$PREFIX/local && make -j4 && make install 2>&1 | tee iphreeqc-3.7.3-15968_aarch64.log
You can make clang default again by pkg remove gcc-default-11
Successfully run following .py (modified for nonCOM)
"""Compares gypsum solubility for WATEQ4F and Pitzer databases.
"""
# Import standard library modules first.
import os
# Import mode non-COM for LINUX / ANDROID
import phreeqpy.iphreeqc.phreeqc_dll as phreeqc_mod
phreeqc = phreeqc_mod.IPhreeqc('/data/data/com.termux/files/usr/local/lib/libiphreeqc-3.7.3.so')
# Then get third party modules.
import matplotlib.pyplot as plt
def selected_array(db_path, input_string):
"""Load database via non-COM and run input string.
"""
phreeqc.load_database(db_path)
phreeqc.run_string(input_string)
return phreeqc.get_selected_output_array()
def show_results(input_string):
"""Get results for different databases
"""
wateq4f_result = selected_array('wateq4f.dat', input_string)
pitzer_result = selected_array('pitzer.dat', input_string)
# Get data from the arrays.
nacl_conc = [entry[0] for entry in wateq4f_result][1:]
wateq4f_values = [entry[1] for entry in wateq4f_result][1:]
pitzer_values = [entry[1] for entry in pitzer_result][1:]
# Plot
plt.plot(nacl_conc, pitzer_values, 'k', nacl_conc, wateq4f_values,'k--')
plt.axis([0, 6, 0, .06])
plt.legend(('PITZER','WATEQ4F'), loc = (0.4, 0.4))
plt.ylabel('GYPSUM SOLUBILITY, MOLES PER KILOGRAM WATER')
plt.xlabel('NaCl, MOLES PER KILOGRAM WATER')
plt.savefig("Figure2.png")
plt.show()
if __name__ == '__main__':
# This will only run when called as script from the command line
# and not when imported from another script.
INPUT_STRING = """
SOLUTION 1
END
INCREMENTAL_REACTIONS
REACTION
NaCl 1.0
0 60*0.1 moles
EQUILIBRIUM_PHASES
Gypsum
USE solution 1
SELECTED_OUTPUT
-reset false
-total Na S(6)
END"""
show_results(INPUT_STRING)
And this multi-threading example again modified
# -*- coding: utf-8 -*-
"""A coupled advection-reaction model using the COM server.
The sketch below. shows the setup. The coupled model contains a advection and a
reaction model. The reaction model can have one or more Phreeqc calculators.
The calculators can work in parallel using the module `multiprocessing`.
+--------------------------------------------------------------------------+
| CoupledModel |
| +--------------------------------------------------------------------+ |
| | AdvectionModel | |
| | | |
| +--------------------------------------------------------------------+ |
| | ReactionModel | |
| | +--------------------+--------------------+--------------------+ | |
| | | PhreeqcCalculator1 | PhreeqcCalculator2 | PhreeqcCalculator3 | | |
+--------------------------------------------------------------------------+
Author: Mike Müller, mmueller@hydrocomputing.com
"""
import multiprocessing
import os
import time
import matplotlib.pyplot as plt
import phreeqpy.iphreeqc.phreeqc_dll as phreeqc_mod
phreeqc = phreeqc_mod.IPhreeqc('/data/data/com.termux/files/usr/local/lib/libiphreeqc-3.7.3.so')
class CoupledModel(object):
"""This is a coupled advection model.
Since it is just a simple example, we use a 1D model.
The same approach can be applied to a 2d or 3D model
as long as the advection model supports it.
The PHREEQC part is the same.
Furthermore, instead of a simple advection model,
we can have a more sophisticated transport model.
"""
def __init__(self, ncells, nshifts, initial_conditions, processes):
self.nshifts = nshifts
self.reaction_model = ReactionModel(ncells, initial_conditions,
processes)
self.reaction_model.make_initial_state()
init_conc = dict([(name, [value] * ncells) for name, value in
self.reaction_model.init_conc.items()])
self.advection_model = AdvectionModel(init_conc,
self.reaction_model.inflow_conc)
self.component_names = self.reaction_model.component_names
self.results = {}
for name in self.component_names:
self.results[name] = []
def run(self):
"""Go over all time steps (shifts).
"""
for shift in range(self.nshifts):
self.advection_model.advect()
self.advection_model.save_results(self.results)
self.reaction_model.modify(self.advection_model.conc)
self.advection_model.update_conc(self.reaction_model.conc)
self.reaction_model.finish()
class AdvectionModel(object):
"""Very simple 1D advection model.
This model can be replaced by a more sophisticated transport
model with two or three dimensions.
This could be an external code such as MT3D or anything
else that is moving concentrations through the subsurface.
"""
def __init__(self, init_conc, inflow_conc):
"""Set the initial and inflow concentrations.
Both concentrations are dictionaries with the specie names as keys.
Values are 1D arrays for `init_conc` and scalars for `inflow_conc`.
"""
self.conc = init_conc
self.inflow_conc = inflow_conc
self.outflow = {}
def update_conc(self, new_conc):
"""Update the concentrations after the reactions step.
This is very simple but could be more involved
if the transport model is more complex.
"""
self.conc = new_conc
def advect(self):
"""Shift one cell.
"""
for name in self.conc:
self.outflow[name] = self.conc[name][-1]
self.conc[name][1:] = self.conc[name][:-1]
self.conc[name][0] = self.inflow_conc[name]
def save_results(self, results):
"""Save the calculation results.
Typically, we would write our results into a file.
For simplicity we just add the current outflow that
we stored in `self.outflow` and add it to `results`,
which is a dictionary with all specie names as keys
and lists as values.
"""
for name in self.conc:
results[name].append(self.outflow[name])
class ReactionModel(object):
"""Calculate reactions using PHREEQC as computational engine.
We have no direct contact with IPhreeqc here.
We make one or more instances of `PhreeqcCalculator`
that are actually using IPhreeqc.
We can use more than one processor with `multiprocessing`.
"""
def __init__(self, ncells, initial_conditions, processes):
if processes > ncells:
raise ValueError('Number of processes needs to be less or equal '
'than number of cells. %d processes %d cells.'
% (processes, ncells))
if processes < 1:
raise ValueError('Need at least one process got %d' % processes)
self.parallel = False
if processes > 1:
self.parallel = True
self.ncells = ncells
self.initial_conditions = initial_conditions
self.processes = processes
self.inflow_conc = {}
self.init_conc = {}
self.conc = {}
self.component_names = []
self.calculators = []
self.cell_ranges = []
self._init_calculators()
self.make_initial_state()
def _init_calculators(self):
"""If we are going parallel we need several calculators.
"""
if self.parallel:
# Domain decomposition.
slave_ncells, reminder = divmod(self.ncells, self.processes)
root_ncells = slave_ncells + reminder
current_cell = root_ncells
root_calculator = PhreeqcCalculator(root_ncells,
self.initial_conditions)
self.calculators = [root_calculator]
self.cell_ranges = [(0, root_ncells)]
for process in range(self.processes - 1):
self.calculators.append(PhreeqcCalculatorProxy(slave_ncells,
self.initial_conditions))
self.cell_ranges.append((current_cell,
current_cell + slave_ncells))
current_cell += slave_ncells
assert current_cell == self.ncells
self.calculators.reverse()
self.cell_ranges.reverse()
else:
root_calculator = PhreeqcCalculator(self.ncells,
self.initial_conditions)
# Just one calculator and the entire range but still use a list
# to provide the same interface as the parallel case.
self.calculators = [root_calculator]
self.cell_ranges = [(0, self.ncells)]
def make_initial_state(self):
"""Get the initial values from the calculator(s).
"""
self.inflow_conc = self.calculators[0].inflow_conc
self.init_conc = self.calculators[0].init_conc
self.component_names = self.calculators[0].component_names
if self.parallel:
# Make sure all calculators are initialized the same.
for calculator in self.calculators[1:]:
assert self.inflow_conc == calculator.inflow_conc
assert self.init_conc == calculator.init_conc
assert self.component_names == calculator.component_names
def modify(self, new_conc):
"""Pass new conc after advection to the calculator.
"""
self.conc = {}
for name in self.component_names:
self.conc[name] = []
for cell_range, calculator in zip(self.cell_ranges, self.calculators):
current_conc = dict([(name, value[cell_range[0]:cell_range[1]]) for
name, value in new_conc.items()])
calculator.modify(current_conc)
for calculator in self.calculators:
conc = calculator.get_modified()
for name in self.component_names:
self.conc[name].extend(conc[name])
def finish(self):
"""This is necessary for multiprocessing.
Multiprocessing uses external processes. These need to be
explicitly closed to avoid hanging of the program at
the end.
"""
for calculator in self.calculators:
calculator.finish()
class PhreeqcCalculator(object):
"""All PHREEQC calculations happen here.
This is the only place where we interact wit IPhreeqc.
Each instance of this class might run in a different
process using `multiprocessing`.
"""
def __init__(self, ncells, initial_conditions):
"""
ncells - number of cells
initial_conditions - string containing PHREEQC input for
solution and exchange, see example below
"""
self.ncells = ncells
self.initial_conditions = initial_conditions
self.inflow_conc = {}
self.init_conc = {}
self.conc = {}
self.phreeqc = phreeqc_mod.IPhreeqc('/data/data/com.termux/files/usr/local/lib/libiphreeqc-3.7.3.so')
self.phreeqc.load_database(r"phreeqc.dat")
self.components = []
self.component_names = []
self._make_initial_state()
def _make_initial_state(self):
"""Copy solution to all cells and calculate initial conditions.
"""
self.phreeqc.run_string(self.initial_conditions)
self.components = self.phreeqc.get_component_list()
start = 1
end = self.ncells
code = ''
code += "COPY solution 1 %d-%d\n" % (start, end)
code += "COPY exchange 1 %d-%d\n" % (start, end)
code += "END\n"
code += "RUN_CELLS; -cells %d-%d\n" % (start, end)
code += self.make_selected_output(self.components)
self.phreeqc.run_string(code)
self.conc = self.get_selected_output()
all_names = self.conc.keys()
self.component_names = [name for name in all_names if name not in
('cb', 'H', 'O')]
code = ''
code += self.make_selected_output(self.components)
code += "RUN_CELLS; -cells 0-1\n"
self.phreeqc.run_string(code)
start_conc = self.get_selected_output()
for name in self.component_names:
self.inflow_conc[name] = start_conc[name][0]
self.init_conc[name] = start_conc[name][1]
def modify(self, new_conc):
"""Set new concentration after advection and re-calculate.
"""
conc = self.conc
end = self.ncells + 1
conc.update(new_conc)
modify = []
for index, cell in enumerate(range(1, end)):
modify.append("SOLUTION_MODIFY %d" % cell)
modify.append("\t-cb %e" % conc['cb'][index])
modify.append("\t-total_h %s" % conc['H'][index])
modify.append("\t-total_o %s" % conc['O'][index])
modify.append("\t-totals")
for name in self.component_names:
modify.append("\t\t%s\t%s" % (name, conc[name][index]))
modify.append("RUN_CELLS; -cells %d-%d\n" % (1, self.ncells))
code = '\n'.join(modify)
self.phreeqc.run_string(code)
self.conc = self.get_selected_output()
def get_modified(self):
"""Return calculated conc.
"""
return self.conc
@ staticmethod # this is just a function but belongs here
def make_selected_output(components):
"""
Build SELECTED_OUTPUT data block.
"""
headings = "-headings cb H O "
headings += '\t'.join(components)
selected_output = """
SELECTED_OUTPUT
-reset false
USER_PUNCH
"""
selected_output += headings + "\n"
# 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 get_selected_output(self):
"""Return calculation result as dict.
Header entries are the keys and the columns
are the values as lists of numbers.
"""
output = self.phreeqc.get_selected_output_array()
header = output[0]
conc = {}
for head in header:
conc[head] = []
for row in output[1:]:
for col, head in enumerate(header):
conc[head].append(row[col])
return conc
def finish(self):
"""Placeholder to give same interface as the multiprocessing version.
"""
pass
class PhreeqcCalculatorProxy(object):
"""Proxy that communicates with other processes.
We uses this proxy for parallel computations.
All code that is specific for parallel computing is located
in here.
"""
def __init__(self, ncells, initial_conditions):
"""Go parallel.
"""
self.in_queue = multiprocessing.JoinableQueue()
self.out_queue = multiprocessing.JoinableQueue()
self.process = multiprocessing.Process(
target=process_worker,
args=(ncells, initial_conditions, self.in_queue, self.out_queue))
self.process.start()
(self.inflow_conc,
self.init_conc,
self.component_names) = self.out_queue.get()
def modify(self, new_conc):
"""Run PHREEQC in another process.
"""
self.in_queue.put(new_conc)
def get_modified(self):
"""Return calculated conc.
"""
return self.out_queue.get()
def finish(self):
"""Terminate the process.
"""
self.in_queue.put(None)
self.process.join()
def process_worker(ncells, initial_conditions, in_queue, out_queue):
"""This runs in another process.
"""
print ('Started process with ID', os.getpid())
calculator = PhreeqcCalculator(ncells, initial_conditions)
out_queue.put((calculator.inflow_conc, calculator.init_conc,
calculator.component_names))
while True:
new_conc = in_queue.get()
# None is the sentinel. We are done
if new_conc is None:
break
calculator.modify(new_conc)
out_queue.put(calculator.conc)
def plot(ncells, outflow, specie_names):
"""Plot the results.
"""
colors = {'Ca': 'r', 'Cl': 'b', 'K': 'g', 'N': 'y', 'Na': 'm'}
x = [i / float(ncells) for i in
range(1, len(outflow[specie_names[0]]) + 1)]
args = []
for name in specie_names:
args.extend([x, outflow[name], colors[name]])
# pylint: disable-msg=W0142
# we do want *
plt.plot(*args)
plt.legend(specie_names, loc=(0.8, 0.5))
plt.ylabel('MILLIMOLES PER KILOGRAM WATER')
plt.xlabel('PORE VOLUME')
plt.show()
def measure_time(func, *args, **kwargs):
"""Convenience function to measure run times.
"""
import sys
if sys.platform == 'win32':
# time.clock is more accurate on Windows
# deprecation warning
timer_func = time.clock
else:
# but behaves differently on other platforms
timer_func = time.process_time
start = timer_func()
result = func(*args, **kwargs)
return result, time.process_time() - start
if __name__ == '__main__':
def main(ncells, nshifts, processes=2):
"""
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
"""
def run():
"""Do the work.
"""
model = CoupledModel(ncells, nshifts, initial_conditions,
processes)
model.run()
return model, model.results
(model, outflow), run_time = measure_time(run)
print ('Statistics')
print ('==========')
print ('number of cells: ', ncells)
print ('number of shifts: ', nshifts)
print ('number of processes:', processes)
print ('run_time: ', run_time)
plot(ncells, outflow, model.component_names)
main(ncells=400, nshifts=1200, processes=2)
results
~/.../com/python $ python parallel_advect_nonCOM.py
Started process with ID 11153
Statistics
==========
number of cells: 400
number of shifts: 1200
number of processes: 2
run_time: 34.015152873999995
Build log. Incorporated only first patch from https://github.com/ufz/iphreeqc/issues/2#issuecomment-867954005
Install log
While loading library it gives this error
Android .so check