源程序
```python
"""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 compatibility adjustment.
if sys.version_info[0] == 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
"""
headings = "-headings cb H O "
for i in range(len(components)):
headings += components[i] + "\t"
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 initialize(cells, first=False):
"""
Initialize IPhreeqc module
"""
phreeqc = phreeqc_mod.IPhreeqc(r'D:\Program Files\USGS\IPhreeqcCOM 3.7.3-15968')
phreeqc.load_database(r"phreeqc.dat")
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][0]
else:
inflow[name] = conc[name][1]
initial[name] = conc[name][1]
task = initial_conditions + "\n"
task += "COPY solution 1 %d-%d\n" % (cells[0], cells[1])
task += "COPY exchange 1 %d-%d\n" % (cells[0], cells[1])
task += "END\n"
task += "RUN_CELLS; -cells %d-%d\n" % (cells[0], cells[1])
task += selected_output
phreeqc.run_string(task)
conc = get_selected_output(phreeqc)
for name in conc:
value = [initial[name]] * len(conc[name])
conc[name] = value
return phreeqc, inflow, conc
def advect_step(phreeqc, inflow, conc, cells):
"""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][0] = inflow[name]
modify = []
for index, cell in enumerate(range(cells[0], cells[1] + 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[0], cells[1]))
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()
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 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):
# advect
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')
header = outflow.keys()
for head in header:
fobj.write('%20s' % head)
fobj.write('\n')
for lineno in range(len(outflow[head])):
for head in header:
fobj.write('%20.17f' % outflow[head][lineno])
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)
出现报错
```python
D:\Users\99225\anaconda3\envs\py3.8\python.exe E:/python/test/advect.py
Traceback (most recent call last):
File "E:/python/test/advect.py", line 221, in <module>
main(ncells=40, shifts=120)
File "E:/python/test/advect.py", line 213, in main
outflow, run_time = measure_time(run, ncells, shifts, specie_names)
File "E:/python/test/advect.py", line 205, in measure_time
result = func(*args, **kwargs)
File "E:/python/test/advect.py", line 171, in run
phreeqc, inflow, conc = initialize(cells, first=True)
File "E:/python/test/advect.py", line 93, in initialize
phreeqc = phreeqc_mod.IPhreeqc(r'D:\Program Files\USGS\IPhreeqcCOM 3.7.3-15968')
File "D:\Users\99225\anaconda3\envs\py3.8\lib\site-packages\phreeqpy\iphreeqc\phreeqc_dll.py", line 47, in __init__
phreeqc = ctypes.cdll.LoadLibrary(dll_path)
File "D:\Users\99225\anaconda3\envs\py3.8\lib\ctypes\__init__.py", line 451, in LoadLibrary
return self._dlltype(name)
File "D:\Users\99225\anaconda3\envs\py3.8\lib\ctypes\__init__.py", line 373, in __init__
self._handle = _dlopen(self._name, mode)
PermissionError: [WinError 5] 拒绝访问。
Dimensions
==========
number of cells: 40
number of shifts 120
Process finished with exit code 1
1、先确认这个目录是否存在
r'D:\Program Files\USGS\IPhreeqcCOM 3.7.3-15968'
2、如果目录存在，查看一下目录属性，是否有只读权限？
3、试试切换到管理员模式的命令行窗口
运行
D:\Users\99225\anaconda3\envs\py3.8\python.exe E:/python/test/advect.py