diff --git a/.github/workflows/pipeline.yml b/.github/workflows/pipeline.yml index f1adfcc..190a984 100644 --- a/.github/workflows/pipeline.yml +++ b/.github/workflows/pipeline.yml @@ -10,23 +10,50 @@ on: jobs: test-job: runs-on: ubuntu-latest + + defaults: + run: + shell: bash -l {0} + strategy: matrix: - python-version: ["3.8", "3.9", "3.10", "3.11"] + python-version: ["3.9", "3.10", "3.11"] steps: - uses: actions/checkout@v3 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + + - name: Get Date + id: get-date + run: echo "today=$(/bin/date -u '+%Y%m%d')" >> $GITHUB_OUTPUT + shell: bash + + - name: Setup Conda Environment + uses: conda-incubator/setup-miniconda@v3 with: - python-version: ${{ matrix.python-version }} + miniforge-variant: Mambaforge + use-mamba: true + activate-environment: cadet-python + channels: conda-forge, + + - name: Cache conda + uses: actions/cache@v3 + env: + # Increase this value to reset cache if environment.yml has not changed + CACHE_NUMBER: 0 + with: + path: ${{ env.CONDA }}/envs + key: python_${{ matrix.python-version }}-${{ steps.get-date.outputs.today }}-${{ env.CACHE_NUMBER }} + - name: Test Wheel install and import run: | + mamba install python==${{ matrix.python-version }} pip install wheel python setup.py bdist_wheel cd dist pip install CADET_Python*.whl python -c "import cadet; print(cadet.__version__)" + cd .. - name: Test with pytest run: | - pip install pytest - pytest tests --rootdir=tests + pip install .[testing] + mamba install cadet -c conda-forge + pytest tests --rootdir=tests -m "not slow" diff --git a/.gitignore b/.gitignore index baf0fe2..c5426f5 100644 --- a/.gitignore +++ b/.gitignore @@ -1,23 +1,12 @@ -################################################################################ -# This .gitignore file was automatically created by Microsoft(R) Visual Studio. -################################################################################ - + /.vs *.h5 +*.csv /__pycache__ -/fix_h5.py -/HICWang.py -/Pilot_Bypass_noBT_noC.py -/Test.py -/test_lwe.py -/comp2_1.csv -/comp2_2.csv -/comp2_comb.csv -/comp2_fraction.csv -/comp_1.csv -/ron.py /CADET.egg-info /dist /build/lib/cadet /cadet/__pycache__ /CADET_Python.egg-info +.idea +*tmp* \ No newline at end of file diff --git a/cadet/cadet.py b/cadet/cadet.py index 1a3530f..a4257fe 100644 --- a/cadet/cadet.py +++ b/cadet/cadet.py @@ -1,22 +1,25 @@ -from addict import Dict - +import copy +import json +import os +import platform +import pprint +import shutil +import subprocess import warnings +from pathlib import Path + with warnings.catch_warnings(): - warnings.filterwarnings("ignore",category=FutureWarning) + warnings.filterwarnings("ignore", category=FutureWarning) import h5py import numpy -import subprocess -import pprint -import copy -import json +from addict import Dict import filelock import contextlib -from pathlib import Path - from cadet.cadet_dll import CadetDLL + class H5(): pp = pprint.PrettyPrinter(indent=4) @@ -117,11 +120,17 @@ def __setitem__(self, key, value): obj = obj[i] obj[parts[-1]] = value + def is_dll(value): suffix = Path(value).suffix return suffix in {'.so', '.dll'} + class CadetMeta(type): + """ + A meta class for the CADET interface. This allows calls to Cadet.cadet_path = "..." to set + the cadet_path for all subsequent Cadet() instances. + """ _cadet_runner_class = None _is_file_class = None @@ -136,15 +145,11 @@ def cadet_path(cls): @cadet_path.setter def cadet_path(cls, value): - def __init__(cls): - cls._cadet_runner_class = None - cls._is_file_class = True - if cls._cadet_runner_class is not None and cls._cadet_runner_class.cadet_path != value: del cls._cadet_runner_class if is_dll(value): - cls._cadet_runner_class = CadetDLL(value) + cls._cadet_runner_class = CadetDLL(value) cls._is_file_class = False else: cls._cadet_runner_class = CadetFile(value) @@ -154,13 +159,23 @@ def __init__(cls): def cadet_path(cls): del cls._cadet_runner_class + class Cadet(H5, metaclass=CadetMeta): - #cadet_path must be set in order for simulations to run + # cadet_path must be set in order for simulations to run def __init__(self, *data): super().__init__(*data) self._cadet_runner = None + + self._is_file = None # Is CLI or DLL + # self.cadet_path # from Bill, declared in meta class -> path to CLI-file or DLL-file + + self.install_path = None # from Jo -> root of the CADET installation. + + self.cadet_cli_path = None + self.cadet_dll_path = None + self.cadet_create_lwe_path = None + self.return_information = None - self._is_file = None @property def is_file(self): @@ -173,8 +188,10 @@ def is_file(self): def cadet_runner(self): if self._cadet_runner is not None: return self._cadet_runner - if self._cadet_runner_class is not None: + if hasattr(self, "_cadet_runner_class") and self._cadet_runner_class is not None: return self._cadet_runner_class + self.autodetect_cadet() + return self._cadet_runner @property def cadet_path(self): @@ -188,7 +205,7 @@ def cadet_path(self, value): del self._cadet_runner if is_dll(value): - self._cadet_runner = CadetDLL(value) + self._cadet_runner = CadetDLL(value) self._is_file = False else: self._cadet_runner = CadetFile(value) @@ -198,6 +215,137 @@ def cadet_path(self, value): def cadet_path(self): del self._cadet_runner + def autodetect_cadet(self): + """ + Autodetect installation CADET based on operating system and API usage. + + Returns + ------- + cadet_root : Path + Installation path of the CADET program. + """ + executable = 'cadet-cli' + if platform.system() == 'Windows': + executable += '.exe' + + # Searching for the executable in system path + path = shutil.which(executable) + + if path is None: + raise FileNotFoundError( + "Could not autodetect CADET installation. Please provide path." + ) + + cli_path = Path(path) + + cadet_root = None + if cli_path is not None: + cadet_root = cli_path.parent.parent + self.install_path = cadet_root + + return cadet_root + + @property + def install_path(self): + """str: Path to the installation of CADET. + + This can either be the root directory of the installation or the path to the + executable file 'cadet-cli'. If a file path is provided, the root directory will + be inferred. + + Raises + ------ + FileNotFoundError + If CADET cannot be found at the specified path. + + Warnings + -------- + If the specified install_path is not the root of the CADET installation, it will + be inferred from the file path. + + See Also + -------- + check_cadet + """ + return self._install_path + + @install_path.setter + def install_path(self, install_path): + """ + Set the installation path of CADET. + + Parameters + ---------- + install_path : str or Path + Path to the root of the CADET installation. + It should either be the root directory of the installation or the path + to the executable file 'cadet-cli'. + If a file path is provided, the root directory will be inferred. + """ + if install_path is None: + self._install_path = None + self.cadet_cli_path = None + self.cadet_dll_path = None + self.cadet_create_lwe_path = None + return + + install_path = Path(install_path) + + if install_path.is_file(): + cadet_root = install_path.parent.parent + warnings.warn( + "The specified install_path is not the root of the CADET installation. " + "It has been inferred from the file path." + ) + else: + cadet_root = install_path + + self._install_path = cadet_root + + cli_executable = 'cadet-cli' + lwe_executable = 'createLWE' + + if platform.system() == 'Windows': + cli_executable += '.exe' + lwe_executable += '.exe' + + cadet_cli_path = cadet_root / 'bin' / cli_executable + if cadet_cli_path.is_file(): + self.cadet_cli_path = cadet_cli_path + self.cadet_path = cadet_cli_path + else: + raise FileNotFoundError( + "CADET could not be found. Please check the path" + ) + + cadet_create_lwe_path = cadet_root / 'bin' / lwe_executable + if cadet_create_lwe_path.is_file(): + self.cadet_create_lwe_path = cadet_create_lwe_path.as_posix() + + if platform.system() == 'Windows': + dll_path = cadet_root / 'bin' / 'cadet.dll' + dll_debug_path = cadet_root / 'bin' / 'cadet_d.dll' + else: + dll_path = cadet_root / 'lib' / 'lib_cadet.so' + dll_debug_path = cadet_root / 'lib' / 'lib_cadet_d.so' + + # Look for debug dll if dll is not found. + if not dll_path.is_file() and dll_debug_path.is_file(): + dll_path = dll_debug_path + + # Look for debug dll if dll is not found. + if dll_path.is_file(): + self.cadet_dll_path = dll_path.as_posix() + + if platform.system() != 'Windows': + try: + cadet_lib_path = cadet_root / 'lib' + if cadet_lib_path.as_posix() not in os.environ['LD_LIBRARY_PATH']: + os.environ['LD_LIBRARY_PATH'] += \ + os.pathsep + cadet_lib_path.as_posix() + except KeyError: + os.environ['LD_LIBRARY_PATH'] = cadet_lib_path.as_posix() + def transform(self, x): return str.upper(x) @@ -209,14 +357,14 @@ def load_results(self): if runner is not None: runner.load_results(self) - def run(self, timeout = None, check=None): + def run(self, timeout=None, check=None): data = self.cadet_runner.run(simulation=self.root.input, filename=self.filename, timeout=timeout, check=check) - #self.return_information = data + # self.return_information = data return data - def run_load(self, timeout = None, check=None, clear=True): + def run_load(self, timeout=None, check=None, clear=True): data = self.cadet_runner.run(simulation=self.root.input, filename=self.filename, timeout=timeout, check=check) - #self.return_information = data + # self.return_information = data self.load_results() if clear: self.clear() @@ -227,14 +375,15 @@ def clear(self): if runner is not None: runner.clear() + class CadetFile: def __init__(self, cadet_path): self.cadet_path = cadet_path - def run(self, filename = None, simulation=None, timeout = None, check=None): + def run(self, filename=None, simulation=None, timeout=None, check=None): if filename is not None: - data = subprocess.run([self.cadet_path, filename], timeout = timeout, check=check, capture_output=True) + data = subprocess.run([self.cadet_path, filename], timeout=timeout, check=check, capture_output=True) return data else: print("Filename must be set before run can be used") @@ -245,9 +394,10 @@ def clear(self): def load_results(self, sim): sim.load(paths=["/meta", "/output"], update=True) + def convert_from_numpy(data, func): ans = Dict() - for key_original,item in data.items(): + for key_original, item in data.items(): key = func(key_original) if isinstance(item, numpy.ndarray): item = item.tolist() @@ -264,9 +414,10 @@ def convert_from_numpy(data, func): ans[key] = item return ans -def recursively_load_dict( data, func): + +def recursively_load_dict(data, func): ans = Dict() - for key_original,item in data.items(): + for key_original, item in data.items(): key = func(key_original) if isinstance(item, dict): ans[key] = recursively_load_dict(item, func) @@ -274,6 +425,7 @@ def recursively_load_dict( data, func): ans[key] = item return ans + def set_path(obj, path, value): "paths need to be broken up so that subobjects are correctly made" path = path.split('/') @@ -285,7 +437,8 @@ def set_path(obj, path, value): temp[path[-1]] = value -def recursively_load( h5file, path, func, paths): + +def recursively_load(h5file, path, func, paths): ans = Dict() if paths is not None: for path in paths: @@ -306,8 +459,8 @@ def recursively_load( h5file, path, func, paths): ans[key] = recursively_load(h5file, local_path + '/', func, None) return ans -def recursively_save(h5file, path, dic, func): +def recursively_save(h5file, path, dic, func): if not isinstance(path, str): raise ValueError("path must be a string") if not isinstance(h5file, h5py._hl.files.File): diff --git a/examples/CSTR.py b/examples/CSTR.py deleted file mode 100644 index 9b35d96..0000000 --- a/examples/CSTR.py +++ /dev/null @@ -1,483 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -#use to render results -import matplotlib.pyplot as plt - -import types - -#location of the cadet binary -cadet_location = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(): - obj = types.SimpleNamespace() - obj.filename = "CSTR_V1_C1.h5" - obj.volume = 1 - obj.conc = 1 - obj.linear = 0 - obj.flowIn = 1 - obj.flowOut = 1 - obj.porosity = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V1_C2.h5" - obj.volume = 1 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V1_C2_F2_F1.h5" - obj.volume = 1 - obj.conc = 2 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V10_C1_F2_F1.h5" - obj.volume = 10 - obj.conc = 1 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V10_C2.h5" - obj.volume = 10 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C1_F2_F2.h5" - obj.volume = 100 - obj.conc = 1 - obj.flowIn = 2 - obj.flowOut = 2 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C2.h5" - obj.volume = 100 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C2_F2_F1.h5" - obj.volume = 100 - obj.conc = 2 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_sens.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_con_down.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - 1e-6 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_con_up.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 + 1e-6 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_lin_down.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_lin_up.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 + 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_por_up.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 + 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_por_down.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_filter_up.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 + 1e-6 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_filter_down.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - 1e-6 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - -def createSimulation(obj): - with h5py.File(obj.filename, 'w') as f: - createInput(f, obj) - -def createInput(f, obj): - input = f.create_group("input") - createModel(input, obj) - createReturn(input) - createSolver(input) - createSensitivity(input) - -def createModel(input, obj): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - # Part of CADET 3.1 - solver = model.create_group('solver') - set_int(solver, 'GS_TYPE', 1) - set_int(solver, 'MAX_KRYLOV', 0) - set_int(solver, 'MAX_RESTARTS', 0) - set_double(solver, 'SCHUR_SAFETY', 1.0e-8) - - createConnections(model, obj) - createInlet(model, obj) - createColumn(model, obj) - createOutlet(model) - -def createConnections(model, obj): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections, obj) - -def createSwitch(connections, obj): - """Create a switch in the system. In 3.0 these are very limited but in 3.1 they allow arbitrary connections between unit operations. - The format is a set of 4 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID. If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, obj.flowIn, - 1, 2, -1, -1, obj.flowOut]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model, obj): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 1) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [obj.conc,]) - set_double(sec, 'LIN_COEFF', [obj.linear,]) - set_double(sec, 'QUAD_COEFF', [0.0,]) - set_double(sec, 'CUBE_COEFF', [0.0,]) - -def createColumn(model, obj): - column = model.create_group('unit_001') - - set_double(column, 'INIT_C', [1.0,]) - set_double(column, 'INIT_VOLUME', obj.volume) - set_double(column, 'FLOWRATE_FILTER', [0.0,]) - set_int(column, 'NCOMP', 1) - set_string(column, 'UNIT_TYPE', 'CSTR') - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 1) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_BULK', 0) - set_int(column, 'WRITE_SENS_INLET', 1) - set_int(column, 'WRITE_SENS_OUTLET', 1) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_BULK', 0) - set_int(column, 'WRITE_SOLUTION_INLET', 1) - set_int(column, 'WRITE_SOLUTION_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - - - column = ret.create_group('unit_002') - - set_int(column, 'WRITE_SENS_BULK', 0) - set_int(column, 'WRITE_SENS_INLET', 1) - set_int(column, 'WRITE_SENS_OUTLET', 1) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_BULK', 0) - set_int(column, 'WRITE_SOLUTION_INLET', 1) - set_int(column, 'WRITE_SOLUTION_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 1) - set_double(solver, 'USER_SOLUTION_TIMES', range(0, 101)) - set_int(solver, 'CONSISTENT_INIT_MODE', 1) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 1) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0, 100.0]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-8) - set_double(time_integrator, 'ALGTOL', 1e-12) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-6) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def createSensitivity(input): - sensitivity = input.create_group('sensitivity') - set_int(sensitivity, 'NSENS', 4) - set_string(sensitivity, 'SENS_METHOD', 'ad1') - - param = sensitivity.create_group('param_000') - set_int(param, 'SENS_UNIT', 0) - set_string(param, 'SENS_NAME', 'CONST_COEFF') - set_int(param, 'SENS_COMP', 0) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', 0) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_001') - set_int(param, 'SENS_UNIT', 0) - set_string(param, 'SENS_NAME', 'LIN_COEFF') - set_int(param, 'SENS_COMP', 0) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', 0) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_002') - set_int(param, 'SENS_UNIT', 1) - set_string(param, 'SENS_NAME', 'POROSITY') - set_int(param, 'SENS_COMP', -1) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', -1) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_003') - set_int(param, 'SENS_UNIT', 1) - set_string(param, 'SENS_NAME', 'FLOWRATE_FILTER') - set_int(param, 'SENS_COMP', -1) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', -1) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - -def runSimulation(obj): - proc = subprocess.Popen([cadet_location, obj.filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(obj): - with h5py.File(obj.filename, 'r') as h5: - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - f.tight_layout() - plt.show() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet = np.array(h5['/output/solution/unit_001/SOLUTION_INLET_COMP_000'].value) - - - axis.set_title("Inlet") - axis.plot(solution_times, inlet, 'b-', label="Product") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol', color='b') - axis.tick_params('y', colors='b') - - lines, labels = axis.get_legend_handles_labels() - axis.legend(lines, labels, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet = np.array(h5['/output/solution/unit_001/SOLUTION_OUTLET_COMP_000'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet, 'b-', label="Product") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol', color='b') - axis.tick_params('y', colors='b') - - lines, labels = axis.get_legend_handles_labels() - axis.legend(lines, labels, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Cadet_Python.py b/examples/Cadet_Python.py deleted file mode 100644 index 5559b91..0000000 --- a/examples/Cadet_Python.py +++ /dev/null @@ -1,266 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "f:/temp/LWE.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - - #read sensitivity data - - #sensitivity of protein 1 to nu 1 - #s1 = simulation.root.output.sensitivity.param_000.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 2 to nu 1 - #s2 = simulation.root.output.sensitivity.param_000.unit_001.sens_column_outlet_comp_002 - - #sensitivity of protein 1 to nu 2 - #s3 = simulation.root.output.sensitivity.param_001.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 1 to sigma 1 - #s4 = simulation.root.output.sensitivity.param_002.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 1 to lambda - #s5 = simulation.root.output.sensitivity.param_003.unit_001.sens_column_outlet_comp_001 - - #Sensitvity of first species to loads of all species (except salt) - #s6 = simulation.root.output.sensitivity.param_004.unit_001.sens_column_outlet_comp_001 - - #Sensitvity of first species to velocity - #s7 = simulation.root.output.sensitivity.param_005.unit_001.sens_column_outlet_comp_001 - - - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - root.input.model.unit_001.discretization.use_analytic_jacobian = 1 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 4 - - root.input.model.unit_000.sec_000.const_coeff = [50.0, 1.0, 1.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [100.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.2, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 5.75e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.37 - root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6] - root.input.model.unit_001.init_c = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.init_q = [1200.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.ncomp = 4 - root.input.model.unit_001.par_diffusion = [7e-9, 6.07e-9, 6.07e-9, 6.07e-9] - root.input.model.unit_001.par_porosity = 0.75 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES' - - #root.input.model.unit_001.velocity = 1 - #root.input.model.unit_001.cross_section_area = 4700.352526439483 - root.input.model.unit_001.velocity = 5.75e-4 - - - root.input.model.unit_001.adsorption.is_kinetic = 0 - root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5, 1.59, 7.7] - root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0, 1000.0, 1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 1200.0 - root.input.model.unit_001.adsorption.sma_nu = [0.0, 4.7, 5.29, 3.7] - root.input.model.unit_001.adsorption.sma_sigma = [0.0, 11.83, 10.6, 10.0] - - root.input.model.unit_001.discretization.nbound = [1, 1, 1, 1] - root.input.model.unit_001.discretization.ncol = 100 - root.input.model.unit_001.discretization.npar = 5 - - root.input.model.unit_002.ncomp = 4 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 1500, 1500) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0] - - #sensitivities - root.input.sensitivity.nsens = 0 - root.input.sensitivity.sens_method = 'AD1' - - #nu1 (sensitivity to protein 1 nu) - #root.input.sensitivity.param_000.sens_abstol = 1e-8 - #root.input.sensitivity.param_000.sens_boundphase = 0 - #root.input.sensitivity.param_000.sens_comp = 1 - #root.input.sensitivity.param_000.sens_factor = 1.0 - #root.input.sensitivity.param_000.sens_name = "SMA_NU" - #root.input.sensitivity.param_000.sens_reaction = -1 - #root.input.sensitivity.param_000.sens_section = -1 - #root.input.sensitivity.param_000.sens_unit = 1 - - #nu2 - #root.input.sensitivity.param_001.sens_abstol = 1e-8 - #root.input.sensitivity.param_001.sens_boundphase = 0 - #root.input.sensitivity.param_001.sens_comp = 2 - #root.input.sensitivity.param_001.sens_factor = 1.0 - #root.input.sensitivity.param_001.sens_name = "SMA_NU" - #root.input.sensitivity.param_001.sens_reaction = -1 - #root.input.sensitivity.param_001.sens_section = -1 - #root.input.sensitivity.param_001.sens_unit = 1 - - #sigma1 - #root.input.sensitivity.param_002.sens_abstol = 1e-8 - #root.input.sensitivity.param_002.sens_boundphase = 0 - #root.input.sensitivity.param_002.sens_comp = 1 - #root.input.sensitivity.param_002.sens_factor = 1.0 - #root.input.sensitivity.param_002.sens_name = "SMA_SIGMA" - #root.input.sensitivity.param_002.sens_reaction = -1 - #root.input.sensitivity.param_002.sens_section = -1 - #root.input.sensitivity.param_002.sens_unit = 1 - - #lambda - #root.input.sensitivity.param_003.sens_abstol = 1e-8 - #root.input.sensitivity.param_003.sens_boundphase = -1 - #root.input.sensitivity.param_003.sens_comp = -1 - #root.input.sensitivity.param_003.sens_factor = 1.0 - #root.input.sensitivity.param_003.sens_name = "SMA_LAMBDA" - #root.input.sensitivity.param_003.sens_reaction = -1 - #root.input.sensitivity.param_003.sens_section = -1 - #root.input.sensitivity.param_003.sens_unit = 1 - - #Sensitvity of first species to loads of all species (except salt) - #root.input.sensitivity.param_004.sens_abstol = [1e-8, 1e-8, 1e-8] - #root.input.sensitivity.param_004.sens_boundphase = [-1, -1, -1] - #root.input.sensitivity.param_004.sens_comp = [1, 2, 3] - #root.input.sensitivity.param_004.sens_factor = [1.0, 1.0, 1.0] - #root.input.sensitivity.param_004.sens_name = ["CONST_COEFF", "CONST_COEFF", "CONST_COEFF"] - #root.input.sensitivity.param_004.sens_reaction = [-1, -1, -1] - #root.input.sensitivity.param_004.sens_section = [0, 0, 0] - #root.input.sensitivity.param_004.sens_unit = [0, 0, 0] - - #Sensitvity of first species to velocity - #root.input.sensitivity.param_005.sens_abstol = 1e-8 - #root.input.sensitivity.param_005.sens_boundphase = -1 - #root.input.sensitivity.param_005.sens_comp = -1 - #root.input.sensitivity.param_005.sens_factor = 1.0 - #root.input.sensitivity.param_005.sens_name = "VELOCITY" - #root.input.sensitivity.param_005.sens_reaction = -1 - #root.input.sensitivity.param_005.sens_section = -1 - #root.input.sensitivity.param_005.sens_unit = 1 - - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 - inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 - outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Cadet_Python_Scoop.py b/examples/Cadet_Python_Scoop.py deleted file mode 100644 index 7c68007..0000000 --- a/examples/Cadet_Python_Scoop.py +++ /dev/null @@ -1,361 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt - -#parallelization -from scoop import futures - -import itertools -import hashlib - -#location of the cadet binary -cadet_location = "C:\\Users\\kosh_000\\cadet_build\\cadet_3\\cadet\\MS_SMKL_RELEASE\\bin\\cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(args): - ka, nu, sigma = args - #filename = "%.2g_%.2g_%.2g.h5" % (ka, nu, sigma) - - filename = hashlib.md5(str(args).encode('utf-8','ignore')).hexdigest() + '.h5' - - createSimulation(filename, ka, nu, sigma) - runSimulation(filename) - plotSimulation(filename) - return "Completed Ka: %s Nu: %s Sigma: %s" % (ka, nu, sigma) - -def createSimulation(filename, ka, nu, sigma): - with h5py.File(filename, 'w') as f: - createInput(f, ka, nu, sigma) - -def createInput(f, ka, nu, sigma): - input = f.create_group("input") - createModel(input, ka, nu, sigma) - createReturn(input) - createSolver(input) - -def createModel(input, ka, nu, sigma): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - # Part of CADET 3.1 - set_int(model, 'GS_TYPE', 1) - set_int(model, 'MAX_KRYLOV', 0) - set_int(model, 'MAX_RESTARTS', 0) - set_double(model, 'SCHUR_SAFETY', 1.0e-8) - - createConnections(model) - createInlet(model) - createColumn(model, ka, nu, sigma) - createOutlet(model) - -def createConnections(model): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections) - -def createSwitch(connections): - """Create a switch in the system. In 3.0 these are very limited but in 3.1 they allow arbitrary connections between unit operations. - The format is a set of 4 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID. If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, 1, 2, -1, -1]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 4) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - # CADET 3.1 - set_double(inlet, 'FLOW', 1.0) - - createLoad(inlet) - createWash(inlet) - createElute(inlet) - -def createLoad(inlet): - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [50.0, 1.0, 1.0, 1.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createWash(inlet): - sec = inlet.create_group('sec_001') - - set_double(sec, 'CONST_COEFF', [50.0, 0.0, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createElute(inlet): - sec = inlet.create_group('sec_002') - - set_double(sec, 'CONST_COEFF', [100.0, 0.0, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.2, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createColumn(model, ka, nu, sigma): - column = model.create_group('unit_001') - - set_string(column, 'ADSORPTION_MODEL', 'STERIC_MASS_ACTION') - set_double(column, 'COL_DISPERSION', 5.75e-8) - set_double(column, 'COL_LENGTH', 0.014) - set_double(column, 'COL_POROSITY', 0.37) - set_double(column, 'FILM_DIFFUSION', [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6]) - set_double(column, 'INIT_C', [50.0, 0.0, 0.0, 0.0]) - set_double(column, 'INIT_Q', [1200.0, 0.0, 0.0, 0.0]) - set_int(column, 'NCOMP', 4) - set_double(column, 'PAR_DIFFUSION', [7e-10, 6.07e-11, 6.07e-11, 6.07e-11]) - set_double(column, 'PAR_POROSITY', 0.75) - set_double(column, 'PAR_RADIUS', 4.5e-5) - set_double(column, 'PAR_SURFDIFFUSION', [0.0, 0.0, 0.0, 0.0]) - set_string(column, 'UNIT_TYPE', 'GENERAL_RATE_MODEL') - set_double(column, 'VELOCITY', 5.75e-4) - - # CADET 3.1 - set_double(column, 'FLOW', 1.0) - - createAdsorption(column, ka, nu, sigma) - createDiscretization(column) - -def createAdsorption(column, ka, nu, sigma): - ads = column.create_group('adsorption') - - set_int(ads, 'IS_KINETIC', 0) - set_double(ads, 'SMA_KA', [0.0, ka, 1.59, 7.7]) - set_double(ads, 'SMA_KD', [0, 1000, 1000, 1000]) - set_double(ads, 'SMA_LAMBDA', 1200) - set_double(ads, 'SMA_NU', [0.0, nu, 5.29, 3.7]) - set_double(ads, 'SMA_SIGMA', [0, sigma, 10.6, 10.0]) - -def createDiscretization(column): - disc = column.create_group('discretization') - - set_int(disc, 'GS_TYPE', 1) - set_int(disc, 'MAX_KRYLOV', 0) - set_int(disc, 'MAX_RESTARTS', 10) - set_int(disc, 'NBOUND', [1, 1, 1, 1]) - set_int(disc, 'NCOL', 50) - set_int(disc, 'NPAR', 10) - set_string(disc, 'PAR_DISC_TYPE', 'EQUIDISTANT_PAR') - set_double(disc, 'SCHUR_SAFETY', 1.0e-8) - set_int(disc, 'USE_ANALYTIC_JACOBIAN', 1) - - createWENO(disc) - -def createWENO(disc): - weno = disc.create_group("weno") - set_int(weno, 'BOUNDARY_MODEL', 0) - set_double(weno, 'WENO_EPS', 1e-10) - set_int(weno, 'WENO_ORDER', 3) - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 4) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - - # CADET 3.1 - set_double(outlet, 'FLOW', 1.0) - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_COLUMN', 0) - set_int(column, 'WRITE_SENS_COLUMN_INLET', 0) - set_int(column, 'WRITE_SENS_COLUMN_OUTLET', 0) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_COLUMN', 0) - set_int(column, 'WRITE_SOLUTION_COLUMN_INLET', 1) - set_int(column, 'WRITE_SOLUTION_COLUMN_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 1) - set_double(solver, 'USER_SOLUTION_TIMES', range(0, 1500+1)) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 3) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0, 10, 90, 1500]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-8) - set_double(time_integrator, 'ALGTOL', 1e-12) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-6) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def runSimulation(filename): - proc = subprocess.Popen([cadet_location, filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(filename): - with h5py.File(filename, 'r') as h5: - fig = plt.figure(figsize=[20, 10]) - - ax1 = fig.add_subplot(1, 2, 1) - ax2 = fig.add_subplot(1, 2, 2) - - #f, (ax1, ax2) = fig.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - fig.tight_layout() - #plt.show() - plt.savefig(filename.replace('h5', 'png'), dpi=100) - plt.close() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_000'].value) - inlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_001'].value) - inlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_002'].value) - inlet_p3 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_003'].value) - - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_000'].value) - outlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_001'].value) - outlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_002'].value) - outlet_p3 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_003'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - - ka = [30, 35, 40] - nu = [4.6, 4.7, 4.8] - sigma = [11.7, 11.8, 11.9] - - settings = itertools.product(ka, nu, sigma) - - results = futures.map(main, settings) - - for result in results: - print(result) diff --git a/examples/HICWang.py b/examples/HICWang.py deleted file mode 100644 index 53bce64..0000000 --- a/examples/HICWang.py +++ /dev/null @@ -1,172 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -#Cadet.cadet_path = r"C:\Users\kosh_000\cadet_build\CADET-fran\MS_SMKL_DEBUG\bin\cadet-cli.exe" -Cadet.cadet_path = r"C:\Users\kosh_000\cadet_build\CADET-fran\MS_SMKL_RELEASE\bin\cadet-cli.exe" - - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "F:/temp/HICWang.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - - root.input.solver.time_integrator.abstol = 1e-6 - root.input.solver.time_integrator.reltol = 0.0 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 2 - - root.input.model.unit_000.sec_000.const_coeff = [4000.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [4000.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [-4000/15000, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'HICWANG' - root.input.model.unit_001.col_dispersion = 5.6e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.4 - root.input.model.unit_001.film_diffusion = [1.4e-7, 1.4e-7] - root.input.model.unit_001.init_c = [4000.0, 0.0] - root.input.model.unit_001.init_q = [0.0] - root.input.model.unit_001.ncomp = 2 - root.input.model.unit_001.par_diffusion = [4e-12, 4e-12] - root.input.model.unit_001.par_porosity = 0.98 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0] - root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - - #root.input.model.unit_001.velocity = 1 - #root.input.model.unit_001.cross_section_area = 4700.352526439483 - root.input.model.unit_001.velocity = 60/(3600*100) - - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.hicwang_kkin = [1e2, 1.0/0.16] - root.input.model.unit_001.adsorption.hicwang_keq = [0.0, 34] - root.input.model.unit_001.adsorption.hicwang_nu = [1.0, 9.5] - root.input.model.unit_001.adsorption.hicwang_qmax = [1.0, 1.3e-2*1000] - root.input.model.unit_001.adsorption.hicwang_beta0 = 3.6e-2 - root.input.model.unit_001.adsorption.hicwang_beta1 = 1.0/1e3 - - root.input.model.unit_001.discretization.nbound = [0, 1] - root.input.model.unit_001.discretization.ncol = 3 - root.input.model.unit_001.discretization.npar = 1 - root.input.model.unit_001.discretization.use_analytic_jacobian = 0 - - - root.input.model.unit_002.ncomp = 2 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 15000, 15000) - root.input.solver.sections.nsec = 2 - root.input.solver.sections.section_continuity = [0] - root.input.solver.sections.section_times = [0.0, 10.0, 15000.0] - - root.input.solver.time_integrator.init_step_size = 1e-2 - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_column_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_column_inlet_comp_001 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_column_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_column_outlet_comp_001 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() - diff --git a/examples/MSSMA.py b/examples/MSSMA.py deleted file mode 100644 index a6ce018..0000000 --- a/examples/MSSMA.py +++ /dev/null @@ -1,166 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -import numpy - -#use to render results -import matplotlib.pyplot as plt - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/cadet3.1-win7-x64/bin/cadet-cli.exe" - -import common - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "MSSMA.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 2 - - root.input.model.unit_000.sec_000.const_coeff = [100.0, 0.15] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [75, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [75, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.05, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'MULTISTATE_STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 1.5E-7 - root.input.model.unit_001.col_length = 0.25 - root.input.model.unit_001.col_porosity = 0.3 - root.input.model.unit_001.film_diffusion = [2.14E-4, 2.1e-5] - root.input.model.unit_001.init_c = [75, 0.0] - root.input.model.unit_001.init_q = [225, 0.0, 0.0] - root.input.model.unit_001.ncomp = 2 - root.input.model.unit_001.par_diffusion = [4.08E-10, 9.0E-12] - root.input.model.unit_001.par_porosity = 0.4 - root.input.model.unit_001.par_radius = 3.25E-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - root.input.model.unit_001.velocity = 0.001 - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.mssma_ka = [0.0, 1E11, 8E6] - root.input.model.unit_001.adsorption.mssma_kd = [0.0, 6E11, 2E16] - root.input.model.unit_001.adsorption.mssma_lambda = 225 - root.input.model.unit_001.adsorption.mssma_nu = [0.0, 10, 25] - root.input.model.unit_001.adsorption.mssma_sigma = [0.0, 48, 66] - root.input.model.unit_001.adsorption.mssma_refq = 225 - root.input.model.unit_001.adsorption.mssma_refc0 = 520.0 - root.input.model.unit_001.adsorption.mssma_rates = [0.0, 0.0, 1e20, 10, 0.0] - - root.input.model.unit_001.discretization.nbound = [1, 2] - root.input.model.unit_001.discretization.ncol = 50 - root.input.model.unit_001.discretization.npar = 5 - - root.input.model.unit_002.ncomp = 2 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 15000, 1000) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 4500, 6100, 15000] - - root.input.solver.consistent_init_mode = 3 - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - - data = numpy.vstack([solution_times, outlet_p1]).transpose() - numpy.savetxt('comp_1.csv', data, delimiter=',') - - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/MSSMA2.py b/examples/MSSMA2.py deleted file mode 100644 index 36d633a..0000000 --- a/examples/MSSMA2.py +++ /dev/null @@ -1,323 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET_31.pdf in the same directory -# - -#Designed for CADET 3.1 - -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -#use to render results -import matplotlib.pyplot as plt - -#location of the cadet binary -cadet_location = "C:/Users/kosh_000/cadet_build/CADET/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(): - filename = "MSSMA2.h5" - createSimulation(filename) - runSimulation(filename) - plotSimulation(filename) - -def createSimulation(filename): - with h5py.File(filename, 'w') as f: - createInput(f) - -def createInput(f): - input = f.create_group("input") - createModel(input) - createReturn(input) - createSolver(input) - -def createModel(input): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - solver = model.create_group('solver') - set_int(solver, 'GS_TYPE', 1) - set_int(solver, 'MAX_KRYLOV', 0) - set_int(solver, 'MAX_RESTARTS', 10) - set_double(solver, 'SCHUR_SAFETY', 1e-8) - - createConnections(model) - createInlet(model) - createColumn(model) - createOutlet(model) - -def createConnections(model): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections) - -def createSwitch(connections): - """Create a switch in the system. In 3.1 they allow arbitrary connections between unit operations. - The format is a set of 5 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID and Volumetric flow rate (m^3/s). If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 3) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - createLoad(inlet) - createWash(inlet) - createElute(inlet) - -def createLoad(inlet): - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [92.0, 0.10631294584377825, 0.10631294584377825/2]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createWash(inlet): - sec = inlet.create_group('sec_001') - - set_double(sec, 'CONST_COEFF', [69.97439960989882, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createElute(inlet): - sec = inlet.create_group('sec_002') - - set_double(sec, 'CONST_COEFF', [69.97439960989882, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.053, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createColumn(model): - column = model.create_group('unit_001') - - set_string(column, 'ADSORPTION_MODEL', 'MULTISTATE_STERIC_MASS_ACTION') - set_double(column, 'COL_DISPERSION', 1.5E-7) - set_double(column, 'COL_LENGTH', 0.215) - set_double(column, 'COL_POROSITY', 0.33999999999999997) - set_double(column, 'FILM_DIFFUSION', [2.14E-4, 2.1e-5, 2.1e-5]) - set_double(column, 'INIT_C', [69.9743996098988, 0.0, 0.0]) - set_double(column, 'INIT_Q', [223.547, 0.0, 0.0, 0.0, 0.0]) - set_int(column, 'NCOMP', 3) - set_double(column, 'PAR_DIFFUSION', [4.08E-10, 9.0E-12, 9.0e-12]) - set_double(column, 'PAR_POROSITY', 0.39) - set_double(column, 'PAR_RADIUS', 3.25E-5) - set_double(column, 'PAR_SURFDIFFUSION', [0.0, 0.0, 0.0, 0.0, 0.0]) - set_string(column, 'UNIT_TYPE', 'GENERAL_RATE_MODEL') - set_double(column, 'VELOCITY', 0.0011437908496732027) - - createAdsorption(column) - createDiscretization(column) - -def createAdsorption(column): - ads = column.create_group('adsorption') - - set_int(ads, 'IS_KINETIC', 1) - set_double(ads, 'MSSMA_KA', [0.0, 1.0652004307518004E31, 7.724553149425915E26, 1.969122487513422E30, 1.1177522067458229E27]) - set_double(ads, 'MSSMA_KD', [0.0, 5.88452172578919E31, 1.955092026422206E36, 9.923200169245614E32, 8.083909678639826E38]) - set_double(ads, 'MSSMA_LAMBDA', 223.547) - set_double(ads, 'MSSMA_NU', [0.0, 9.618977853171593, 24.75290977103934, 6.058214688510013, 20.20231695871297]) - set_double(ads, 'MSSMA_SIGMA', [0.0, 47.82861669713074, 65.93967947378826, 40.22617141805257, 63.71221340773053 -]) - set_double(ads, 'MSSMA_REFQ', 223.547) - set_double(ads, 'MSSMA_REFC0', 520.0) - set_double(ads, 'MSSMA_RATES', [0.0, 0.0, 9.39710359947847E39, 9.503195767335168, 0.0, 0.0, 5.571477811298548E30, 82427.80960452619, 0.0]) - -def createDiscretization(column): - disc = column.create_group('discretization') - - set_int(disc, 'GS_TYPE', 1) - set_int(disc, 'MAX_KRYLOV', 0) - set_int(disc, 'MAX_RESTARTS', 0) - set_int(disc, 'NBOUND', [1, 2, 2]) - set_int(disc, 'NCOL', 100) - set_int(disc, 'NPAR', 5) - set_string(disc, 'PAR_DISC_TYPE', 'EQUIDISTANT_PAR') - set_double(disc, 'SCHUR_SAFETY', 1.0e-8) - set_int(disc, 'USE_ANALYTIC_JACOBIAN', 1) - - createWENO(disc) - -def createWENO(disc): - weno = disc.create_group("weno") - set_int(weno, 'BOUNDARY_MODEL', 0) - set_double(weno, 'WENO_EPS', 1e-10) - set_int(weno, 'WENO_ORDER', 3) - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 3) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_COLUMN', 0) - set_int(column, 'WRITE_SENS_COLUMN_INLET', 0) - set_int(column, 'WRITE_SENS_COLUMN_OUTLET', 0) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_COLUMN', 0) - set_int(column, 'WRITE_SOLUTION_COLUMN_INLET', 1) - set_int(column, 'WRITE_SOLUTION_COLUMN_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 0) - set_int(solver, 'CONSISTENT_INIT_MODE', 3) - set_double(solver, 'USER_SOLUTION_TIMES', np.linspace(0, 14731.2, 1000)) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 3) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0.0, 4445.422740524782, 6103.9941690962105, 14731.2]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-10) - set_double(time_integrator, 'ALGTOL', 1e-10) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-8) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def runSimulation(filename): - proc = subprocess.Popen([cadet_location, filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(filename): - with h5py.File(filename, 'r') as h5: - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - f.tight_layout() - plt.show() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_000'].value) - inlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_001'].value) - inlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_002'].value) - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'r-', label="P2") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_000'].value) - outlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_001'].value) - outlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_002'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'r-', label="P2") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Pilot_Bypass_noBT_noC.py b/examples/Pilot_Bypass_noBT_noC.py deleted file mode 100644 index 1a053b8..0000000 --- a/examples/Pilot_Bypass_noBT_noC.py +++ /dev/null @@ -1,198 +0,0 @@ - -import matplotlib.pyplot as plt - -import numpy - -import pandas - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -common = Cadet() -root = common.root - -root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR' -root.input.model.unit_001.discretization.schur_safety = 1.0e-8 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.weno.boundary_model = 0 -root.input.model.unit_001.discretization.weno.weno_eps = 1e-10 -root.input.model.unit_001.discretization.weno.weno_order = 3 -root.input.model.unit_001.discretization.gs_type = 1 -root.input.model.unit_001.discretization.max_krylov = 0 -root.input.model.unit_001.discretization.max_restarts = 10 - -root.input.solver.time_integrator.abstol = 1e-8 -root.input.solver.time_integrator.algtol = 1e-12 -root.input.solver.time_integrator.init_step_size = 1e-6 -root.input.solver.time_integrator.max_steps = 1000000 -root.input.solver.time_integrator.reltol = 1e-6 - -root.input.model.solver.gs_type = 1 -root.input.model.solver.max_krylov = 0 -root.input.model.solver.max_restarts = 10 -root.input.model.solver.schur_safety = 1e-8 - -#CADET 3.1 and CADET-dev flags are in here so that it works with both -#CADET-dev removed column from the name on the inputs and outputs since for many -#operations it no longer makes sense -root.input['return'].write_solution_times = 1 -root.input['return'].split_components_data = 1 -root.input['return'].unit_000.write_sens_bulk = 0 -root.input['return'].unit_000.write_sens_flux = 0 -root.input['return'].unit_000.write_sens_inlet = 0 -root.input['return'].unit_000.write_sens_outlet = 0 -root.input['return'].unit_000.write_sens_particle = 0 -root.input['return'].unit_000.write_solution_bulk = 0 -root.input['return'].unit_000.write_solution_flux = 0 -root.input['return'].unit_000.write_solution_inlet = 1 -root.input['return'].unit_000.write_solution_outlet = 1 -root.input['return'].unit_000.write_solution_particle = 0 -root.input['return'].unit_000.write_sens_column = 0 -root.input['return'].unit_000.write_sens_column_inlet = 0 -root.input['return'].unit_000.write_sens_column_outlet = 0 -root.input['return'].unit_000.write_solution_column = 0 -root.input['return'].unit_000.write_solution_column_inlet = 1 -root.input['return'].unit_000.write_solution_column_outlet = 1 - -root.input['return'].unit_001 = root.input['return'].unit_000 -root.input['return'].unit_002 = root.input['return'].unit_000 -root.input['return'].unit_003 = root.input['return'].unit_000 - -root.input.solver.nthreads = 1 - - -cstr = Cadet() -root = cstr.root - -root.input.model.unit_002.unit_type = 'CSTR' -root.input.model.unit_002.ncomp = 1 -root.input.model.unit_002.nbound =[0] -root.input.model.unit_002.init_c = [0] -root.input.model.unit_002.porosity = 1.0 -root.input.model.unit_002.init_volume = 5e-6 -root.input.model.unit_002.flowrate_filter = 0.0 - -dpfr = Cadet() -root = dpfr.root - -root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES' -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.adsorption_model= 'NONE' -root.input.model.unit_001.init_c = [0] -root.input.model.unit_001.init_q = [] -root.input.model.unit_001.col_dispersion = 4e-8 -root.input.model.unit_001.col_length = 1.0 -root.input.model.unit_001.total_porosity = 1.0 -root.input.model.unit_001.velocity = 1.0 -root.input.model.unit_001.cross_section_area = ((3.75e-5)**2)*numpy.pi - -root.input.model.unit_001.discretization.nbound = [0] -root.input.model.unit_001.discretization.ncol = 50 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.reconstruction = 'WENO' - - -io = Cadet() -root = io.root - -root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' -root.input.model.unit_000.unit_type = 'INLET' -root.input.model.unit_000.ncomp = 1 - -root.input.model.unit_000.sec_000.const_coeff = [5.0] -root.input.model.unit_000.sec_000.lin_coeff = [0] -root.input.model.unit_000.sec_000.quad_coeff = [0] -root.input.model.unit_000.sec_000.cube_coeff = [0] - -root.input.model.unit_000.sec_001.const_coeff = [1000.0] -root.input.model.unit_000.sec_001.lin_coeff = [0] -root.input.model.unit_000.sec_001.quad_coeff = [0] -root.input.model.unit_000.sec_001.cube_coeff = [0] - -root.input.model.unit_000.sec_002.const_coeff = [0] -root.input.model.unit_000.sec_002.lin_coeff = [0] -root.input.model.unit_000.sec_002.quad_coeff = [0] -root.input.model.unit_000.sec_002.cube_coeff = [0] - -root.input.model.unit_003.unit_type = 'OUTLET' -root.input.model.unit_003.ncomp = 1 - -connectivity = Cadet() -root = connectivity.root - -root.input.model.nunits = 4 - -root.input.model.connections.nswitches = 1 -root.input.model.connections.switch_000.section = 0 -root.input.model.connections.switch_000.connections = [0, 2, -1, -1, (100/(60*1e6)), - 2, 1, -1, -1, (100/(60*1e6)), - 1, 3, -1, -1, (100/(60*1e6))] - -root.input.solver.user_solution_times = numpy.linspace(0, 153.87, 1000) -root.input.solver.sections.nsec = 3 -root.input.solver.sections.section_continuity = [0] -root.input.solver.sections.section_times = [0.0, 15.006, 100.0, 153.87] - - -def main(): - sim = Cadet(common.root, dpfr.root, io.root, cstr.root, connectivity.root) - sim.filename = r"F:\jurgen\Pilot_test.h5" - #createSimulation(sim) - sim.save() - sim.run() - sim.load() - plotSimulation(sim) - - #sim = Cadet(common.root, cstr.root, io.root, connectivity.root) - #sim.root.input['return'].unit_001.write_solution_volume = 1 - #sim.root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.5, - # 1, 2, -1, -1, 1.0] - #sim.filename = r"C:\Users\Kohl\Desktop\Cadet\test_file_cstr.h5" - #sim.save() - #sim.run() - #sim.load() - #plotSimulation(sim) - #plotVolume(sim) - - writer = pandas.ExcelWriter(r'F:\jurgen\test_file_cstr.xlsx') - - inputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_000.solution_inlet_comp_000)]) - outputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_outlet_comp_000)]) - #volumes = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Volume', numpy.squeeze(sim.root.output.solution.unit_001.solution_volume))]) - - inputs.to_excel(writer, 'Input', index=False) - outputs.to_excel(writer, 'Output', index=False) - #volumes.to_excel(writer, 'Volume', index=False) - - writer.save() - - -def plotSimulation(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_002.solution_outlet_comp_000) - plt.show() - -def plotVolume(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_001.solution_volume) - plt.show() - -def createSimulation(sim): - root = sim.root - - #SMA Model - root.input.model.unit_001.adsorption_model= 'STERIC_MASS_ACTION' - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.sma_ka = [0,35.5] - root.input.model.unit_001.adsorption.sma_kd = [0,1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 800.0 - root.input.model.unit_001.adsorption.sma_nu = [1,7.0] - root.input.model.unit_001.adsorption.sma_sigma = [0,10.0] - - #Linear isotherm - # root.input.model.unit_001.adsorption_model = 'LINEAR' - - #root.input.model.unit_001.adsorption.is_kinetic = 1 - #root.input.model.unit_001.adsorption.lin_ka = [5e-3] - #root.input.model.unit_001.adsorption.lin_kd = [1e-3] - -if __name__ == "__main__": - main() diff --git a/examples/Test.py b/examples/Test.py deleted file mode 100644 index 01583f4..0000000 --- a/examples/Test.py +++ /dev/null @@ -1,182 +0,0 @@ -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -import pandas - -common = Cadet() -root = common.root - -root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR' -root.input.model.unit_001.discretization.schur_safety = 1.0e-8 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.weno.boundary_model = 0 -root.input.model.unit_001.discretization.weno.weno_eps = 1e-10 -root.input.model.unit_001.discretization.weno.weno_order = 3 -root.input.model.unit_001.discretization.gs_type = 1 -root.input.model.unit_001.discretization.max_krylov = 0 -root.input.model.unit_001.discretization.max_restarts = 10 - -root.input.solver.time_integrator.abstol = 1e-8 -root.input.solver.time_integrator.algtol = 1e-12 -root.input.solver.time_integrator.init_step_size = 1e-6 -root.input.solver.time_integrator.max_steps = 1000000 -root.input.solver.time_integrator.reltol = 1e-6 - -root.input.model.solver.gs_type = 1 -root.input.model.solver.max_krylov = 0 -root.input.model.solver.max_restarts = 10 -root.input.model.solver.schur_safety = 1e-8 - -#CADET 3.1 and CADET-dev flags are in here so that it works with both -#CADET-dev removed column from the name on the inputs and outputs since for many -#operations it no longer makes sense -root.input['return'].write_solution_times = 1 -root.input['return'].split_components_data = 1 -root.input['return'].unit_000.write_sens_bulk = 0 -root.input['return'].unit_000.write_sens_flux = 0 -root.input['return'].unit_000.write_sens_inlet = 0 -root.input['return'].unit_000.write_sens_outlet = 0 -root.input['return'].unit_000.write_sens_particle = 0 -root.input['return'].unit_000.write_solution_bulk = 0 -root.input['return'].unit_000.write_solution_flux = 0 -root.input['return'].unit_000.write_solution_inlet = 1 -root.input['return'].unit_000.write_solution_outlet = 1 -root.input['return'].unit_000.write_solution_particle = 0 -root.input['return'].unit_000.write_sens_column = 0 -root.input['return'].unit_000.write_sens_column_inlet = 0 -root.input['return'].unit_000.write_sens_column_outlet = 0 -root.input['return'].unit_000.write_solution_column = 0 -root.input['return'].unit_000.write_solution_column_inlet = 1 -root.input['return'].unit_000.write_solution_column_outlet = 1 - -root.input['return'].unit_001 = root.input['return'].unit_000 -root.input['return'].unit_002 = root.input['return'].unit_000 - -root.input.solver.nthreads = 1 - -column_setup = Cadet() -root = column_setup.root - -root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' -root.input.model.unit_001.col_dispersion = 5.75e-8 -root.input.model.unit_001.col_length = 0.014 -root.input.model.unit_001.col_porosity = 0.37 -root.input.model.unit_001.film_diffusion = [6.9e-6] -root.input.model.unit_001.init_c = [0.0] -root.input.model.unit_001.init_q = [0.0] -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.par_diffusion = [7e-10] -root.input.model.unit_001.par_porosity = 0.75 -root.input.model.unit_001.par_radius = 4.5e-5 -root.input.model.unit_001.par_surfdiffusion = [0.0] -root.input.model.unit_001.velocity = 1 -root.input.model.unit_001.cross_section_area = 4700.352526439483 - -root.input.model.unit_001.discretization.nbound = [1] -root.input.model.unit_001.discretization.ncol = 50 -root.input.model.unit_001.discretization.npar = 4 - -cstr = Cadet() -root = cstr.root - -root.input.model.unit_001.unit_type = 'CSTR' -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.nbound = 0 -root.input.model.unit_001.init_c = [1.0] -root.input.model.unit_001.porosity = 1.0 -root.input.model.unit_001.init_volume = 10.0 -root.input.model.unit_001.flowrate_filter = 0.0 - -io = Cadet() -root = io.root - -root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' -root.input.model.unit_000.unit_type = 'INLET' -root.input.model.unit_000.ncomp = 1 - -root.input.model.unit_000.sec_000.const_coeff = [1.0] -root.input.model.unit_000.sec_000.lin_coeff = [0.0] -root.input.model.unit_000.sec_000.quad_coeff = [0.0] -root.input.model.unit_000.sec_000.cube_coeff = [0.0] - -root.input.model.unit_000.sec_001.const_coeff = [0.0] -root.input.model.unit_000.sec_001.lin_coeff = [0.0] -root.input.model.unit_000.sec_001.quad_coeff = [0.0] -root.input.model.unit_000.sec_001.cube_coeff = [0.0] - -root.input.model.unit_002.unit_type = 'OUTLET' -root.input.model.unit_002.ncomp = 1 - -connectivity = Cadet() -root = connectivity.root - -root.input.model.nunits = 3 - -root.input.model.connections.nswitches = 1 -root.input.model.connections.switch_000.section = 0 -root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - -root.input.solver.user_solution_times = numpy.linspace(0, 500, 501) -root.input.solver.sections.nsec = 2 -root.input.solver.sections.section_continuity = [0] -root.input.solver.sections.section_times = [0.0, 100.0, 500.0] - - -def main(): - #sim = Cadet(common.root, column_setup.root, io.root, connectivity.root) - #sim.filename = "F:/temp/test_file.h5" - #createSimulation(sim) - #sim.save() - #sim.run() - #sim.load() - #plotSimulation(sim) - - sim = Cadet(common.root, cstr.root, io.root, connectivity.root) - sim.root.input['return'].unit_001.write_solution_volume = 1 - sim.root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.5, - 1, 2, -1, -1, 1.0] - sim.filename = "F:/temp/test_file_cstr.h5" - sim.save() - sim.run() - sim.load() - plotSimulation(sim) - plotVolume(sim) - - writer = pandas.ExcelWriter('F:/temp/test_file_cstr.xlsx') - - inputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_inlet_comp_000)]) - outputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_outlet_comp_000)]) - volumes = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Volume', numpy.squeeze(sim.root.output.solution.unit_001.solution_volume))]) - - inputs.to_excel(writer, 'Input', index=False) - outputs.to_excel(writer, 'Output', index=False) - volumes.to_excel(writer, 'Volume', index=False) - - writer.save() - - -def plotSimulation(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_002.solution_outlet_comp_000) - plt.show() - -def plotVolume(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_001.solution_volume) - plt.show() - -def createSimulation(sim): - root = sim.root - - - root.input.model.unit_001.adsorption_model = 'LINEAR' - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.lin_ka = [5e-3] - root.input.model.unit_001.adsorption.lin_kd = [1e-3] - -if __name__ == "__main__": - main() diff --git a/examples/__pycache__/common.cpython-37.pyc b/examples/__pycache__/common.cpython-37.pyc deleted file mode 100644 index 186ce56..0000000 Binary files a/examples/__pycache__/common.cpython-37.pyc and /dev/null differ diff --git a/examples/cadet_json.py b/examples/cadet_json.py deleted file mode 100644 index 576c129..0000000 --- a/examples/cadet_json.py +++ /dev/null @@ -1,29 +0,0 @@ -import cadet - -data = cadet.H5() -data.filename = r"F:\MCMC\Synthetic\Bypass\Bypass_MCMC_pump_delay\mcmc_refine\mcmc\mcmc.h5" -data.load(paths=["/bounds_change/json",]) - - -print(data) - -print("bounds json", data.root.bounds_change) - - -print("json", data.root.bounds_change.json) - - - -data.load(paths=["/bounds_change",]) - - - -print("bounds_change", data) - - - -data.load() - - - -print("all", data) \ No newline at end of file diff --git a/examples/cadet_library.py b/examples/cadet_library.py deleted file mode 100644 index c2a4517..0000000 --- a/examples/cadet_library.py +++ /dev/null @@ -1,291 +0,0 @@ -#Python 3.5+ -#Depends on addict https://github.com/mewwts/addict -#Depends on h5py, numpy - -from addict import Dict - -import warnings -with warnings.catch_warnings(): - warnings.filterwarnings("ignore",category=FutureWarning) - import h5py -import numpy -import subprocess -import pprint -import copy -import ctypes - -class Cadet(): - - #cadet_path must be set in order for simulations to run - cadet_path = None - cadet_library_path = None - return_information = None - - pp = pprint.PrettyPrinter(indent=4) - - def __init__(self, *data): - self.root = Dict() - self.filename = None - for i in data: - self.root.update(copy.deepcopy(i)) - - def load(self): - if self.filename is not None: - with h5py.File(self.filename, 'r') as h5file: - self.root = Dict(recursively_load(h5file, '/')) - else: - print('Filename must be set before load can be used') - - load_file = load - - def load_memory(self): - self.root = recursively_load_memory(self.memory) - - def save(self): - if self.filename is not None: - with h5py.File(self.filename, 'w') as h5file: - recursively_save(h5file, '/', self.root) - else: - print("Filename must be set before save can be used") - - save_file = save - - def save_memory(self): - self.memory, self.struct_class = recursively_save_memory(self.root) - - def run(self, timeout = None, check=None): - if self.filename is not None: - data = subprocess.run([self.cadet_path, self.filename], timeout = timeout, check=check, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - self.return_information = data - return data - else: - print("Filename must be set before run can be used") - - run_file = run - - def run_memory(self): - cadet_dll = ctypes.cdll.LoadLibrary(self.cadet_library_path) - print(self.memory) - - cadet_dll.cadetGetLibraryVersion.restype = ctypes.c_char_p - print(cadet_dll.cadetGetLibraryVersion().decode()) - - def __str__(self): - temp = [] - temp.append('Filename = %s' % self.filename) - temp.append(self.pp.pformat(self.root)) - return '\n'.join(temp) - - def update(self, merge): - self.root.update(merge.root) - - def __getitem__(self, key): - key = key.lower() - obj = self.root - for i in key.split('/'): - if i: - obj = obj[i] - return obj - - def __setitem__(self, key, value): - key = key.lower() - obj = self.root - parts = key.split('/') - for i in parts[:-1]: - if i: - obj = obj[i] - obj[parts[-1]] = value - -def recursively_load( h5file, path): - - ans = {} - for key, item in h5file[path].items(): - key = key.lower() - if isinstance(item, h5py._hl.dataset.Dataset): - ans[key] = item.value - elif isinstance(item, h5py._hl.group.Group): - ans[key] = recursively_load(h5file, path + key + '/') - return ans - -def recursively_save( h5file, path, dic): - - # argument type checking - if not isinstance(dic, dict): - raise ValueError("must provide a dictionary") - - if not isinstance(path, str): - raise ValueError("path must be a string") - if not isinstance(h5file, h5py._hl.files.File): - raise ValueError("must be an open h5py file") - # save items to the hdf5 file - for key, item in dic.items(): - key = str(key) - if not isinstance(key, str): - raise ValueError("dict keys must be strings to save to hdf5") - #handle int, float, string and ndarray of int32, int64, float64 - if isinstance(item, str): - h5file[path + key.upper()] = numpy.array(item, dtype='S') - - elif isinstance(item, int): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.int32) - - elif isinstance(item, float): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float64: - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float32: - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int32: - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int64: - h5file[path + key.upper()] = item.astype(numpy.int32) - - elif isinstance(item, numpy.ndarray) and item.dtype.kind == 'S': - h5file[path + key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, int) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.int32) - - elif isinstance(item, list) and any(isinstance(i, float) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.int32): - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.float64): - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.float32): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.bytes_): - h5file[path + key.upper()] = item - - elif isinstance(item, bytes): - h5file[path + key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, str) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype="S") - - # save dictionaries - elif isinstance(item, dict): - recursively_save(h5file, path + key + '/', item) - # other types cannot be saved and will result in an error - else: - raise ValueError('Cannot save %s/%s key with %s type.' % (path, key.upper(), type(item))) - - -def recursively_save_memory(dic): - - # argument type checking - if not isinstance(dic, dict): - raise ValueError("must provide a dictionary") - - class TempClass(ctypes.Structure): - pass - - # save items to the hdf5 file - _fields_ = [] - values = {} - for key, item in dic.items(): - key = str(key) - if not isinstance(key, str): - raise ValueError("dict keys must be strings to save to hdf5") - #handle int, float, string and ndarray of int32, int64, float64 - if isinstance(item, str): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, int): - _fields_.append( (key.upper(), ctypes.c_int ) ) - - values[key.upper()] = item - - elif isinstance(item, float): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float64: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float32: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.float64).ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int32: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int64: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] =item.astype(numpy.int32).ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, numpy.ndarray) and item.dtype.kind == 'S': - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.c_char_p) - - elif isinstance(item, list) and all(isinstance(i, int) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.int32).ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, list) and any(isinstance(i, float) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.float64).ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.int32): - _fields_.append( (key.upper(), ctypes.c_int ) ) - - values[key.upper()] = int(item) - - elif isinstance(item, numpy.float64): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = float(item) - - elif isinstance(item, numpy.float32): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = float(item) - - elif isinstance(item, numpy.bytes_): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, bytes): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, str) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_char_p)) ) - - values[key.upper()] = numpy.array(item, dtype="S").ctypes.data_as(ctypes.POINTER(ctypes.c_char_p)) - - # save dictionaries - elif isinstance(item, dict): - temp_recur, TempClass_recur = recursively_save_memory(item) - _fields_.append( (key, TempClass_recur ) ) - values[key] = temp_recur - # other types cannot be saved and will result in an error - else: - raise ValueError('Cannot save %s/%s key with %s type.' % (path, key.upper(), type(item))) - - TempClass._fields_ = _fields_ - temp = TempClass(**values) - #a = None - return temp, TempClass diff --git a/examples/dll_test.py b/examples/dll_test.py deleted file mode 100644 index 380689a..0000000 --- a/examples/dll_test.py +++ /dev/null @@ -1,40 +0,0 @@ -import cadet.cadet -import cadet.cadet_dll -import time -import matplotlib.pyplot as plt - -print('here') - -cadet.Cadet.cadet_path = r"C:\Users\kosh\cadet_build\CADET\VCPKG_42_dll\bin\cadet.dll" - -sim = cadet.Cadet() -sim.filename = r"C:\Users\kosh\cadet_build\CADET\VCPKG_42_dll\bin\LWE.h5" -sim.load() - -sim.root.input['return'].unit_000.write_solution_inlet = 1 -sim.root.input['return'].unit_000.write_solution_bulk = 1 -sim.root.input['return'].unit_000.single_as_multi_port = 0 -sim.root.input['return'].unit_001.write_solution_bulk = 0 - -print('before run') - -sim.run_load() - -print('after run') - -plt.figure(figsize=[15,15]) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_001) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_002) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_003) -plt.show() - - -plt.figure(figsize=[15,15]) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_001) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_002) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_003) -plt.show() - -print(sim.root.output.solution.unit_000.solution_bulk.shape) - -a = None diff --git a/examples/fix_h5.py b/examples/fix_h5.py deleted file mode 100644 index 25bbaed..0000000 --- a/examples/fix_h5.py +++ /dev/null @@ -1,10 +0,0 @@ -from cadet import Cadet - -fix = Cadet() -fix.filename = r"F:\output\dextran\dextran_pulse.h5" -fix.load() -fix.root.input.model.unit_001.adsorption_model = "LINEAR" -fix.root.input.model.unit_001.adsorption.is_kinetic = 1 -fix.root.input.model.unit_001.adsorption.lin_ka = [0.0] -fix.root.input.model.unit_001.adsorption.lin_kd = [1e3] -fix.save() \ No newline at end of file diff --git a/examples/library_test.py b/examples/library_test.py deleted file mode 100644 index da989ca..0000000 --- a/examples/library_test.py +++ /dev/null @@ -1,11 +0,0 @@ -import cadet_library - -data = cadet_library.Cadet() -data.filename = r"C:\Users\kosh_000\Documents\Visual Studio 2017\Projects\CADETMatch\Examples\Example1\Dextran\dextran_pulse.h5" -data.cadet_library_path = r"C:\Users\kosh_000\cadet_build\CADET-dev\MS_SMKL_RELEASE\bin\cadet.dll" -data.load() - -data.save_memory() - - -data.run_memory() \ No newline at end of file diff --git a/examples/ron.py b/examples/ron.py deleted file mode 100644 index ed2b64a..0000000 --- a/examples/ron.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "ron.h5" - createSimulation(simulation) - simulation.save() - out = simulation.run() - print(out) - simulation.load() - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - root.input.model.unit_001.discretization.use_analytic_jacobian = 1 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 4 - - root.input.model.unit_000.sec_000.const_coeff = [50.0, 1.0, 1.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [100.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.2, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 5.75e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.37 - root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6] - root.input.model.unit_001.init_c = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.init_q = [1200.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.ncomp = 4 - root.input.model.unit_001.par_diffusion = [7e-10, 6.07e-11, 6.07e-11, 6.07e-11] - root.input.model.unit_001.par_porosity = 0.75 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES' - # root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - - root.input.model.unit_001.velocity = 5.75e-4 - - root.input.model.unit_001.adsorption.is_kinetic = 0 - root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5, 1.59, 7.7] - root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0, 1000.0, 1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 1200.0 - root.input.model.unit_001.adsorption.sma_nu = [0.0, 4.7, 5.29, 3.7] - root.input.model.unit_001.adsorption.sma_sigma = [0.0, 11.83, 10.6, 10.0] - - root.input.model.unit_001.discretization.nbound = [1, 1, 1, 1] - root.input.model.unit_001.discretization.ncol = 3 - root.input.model.unit_001.discretization.npar = 1 - - root.input.model.unit_002.ncomp = 4 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 1500, 1500) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0] - - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 - inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 - outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() - diff --git a/setup.cfg b/setup.cfg index 9d5f797..23d72ac 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,7 @@ # Inside of setup.cfg [metadata] -description-file = README.md \ No newline at end of file +description-file = README.md + +[tool:pytest] +markers = + slow: marks tests as slow (deselect with '-m "not slow"') diff --git a/setup.py b/setup.py index 4d11f1d..54cf7e2 100644 --- a/setup.py +++ b/setup.py @@ -24,5 +24,12 @@ "License :: OSI Approved :: BSD License", "Operating System :: OS Independent", ], + extras_require={ + "testing": [ + "pytest", + "matplotlib", + "pandas", + ], + }, python_requires='>=3.7', ) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/common.py b/tests/common.py similarity index 90% rename from examples/common.py rename to tests/common.py index 0caccaa..db6a9fd 100644 --- a/examples/common.py +++ b/tests/common.py @@ -7,11 +7,9 @@ root.input.model.solver.gs_type = 1 root.input.model.solver.max_krylov = 0 root.input.model.solver.max_restarts = 10 -root.input.model.solver.schur_safety = 1e-8 +root.input.model.solver.schur_safety = 1e-8 -#CADET 3.1 and CADET-dev flags are in here so that it works with both -#CADET-dev removed column from the name on the inputs and outputs since for many -#operations it no longer makes sense +# CADET 3.1 and CADET-dev flags are in here so that it works with both root.input['return'].write_solution_times = 1 root.input['return'].split_components_data = 1 root.input['return'].unit_000.write_sens_bulk = 0 diff --git a/examples/MSSMA_2comp.py b/tests/test_MSSMA_2comp.py similarity index 83% rename from examples/MSSMA_2comp.py rename to tests/test_MSSMA_2comp.py index ab3f5fa..06f6713 100644 --- a/examples/MSSMA_2comp.py +++ b/tests/test_MSSMA_2comp.py @@ -1,31 +1,17 @@ -#!/usr/bin/env python3.6 +import os -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -import numpy - -#use to render results import matplotlib.pyplot as plt +import numpy +import pandas +import pytest from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/cadet3.1-win7-x64/bin/cadet-cli.exe" +from tests import common -import common -import pandas def gen_fraction_times(start, stop, bins): - return numpy.linspace(start, stop, bins+1, endpoint=True) + return numpy.linspace(start, stop, bins + 1, endpoint=True) + def gen_fractions(fraction_times, sim): nComp = sim.root.input.model.unit_000.ncomp @@ -36,24 +22,29 @@ def gen_fractions(fraction_times, sim): for idx, (start, stop) in enumerate(zip(fraction_times[:-1], fraction_times[1:])): selected = (times >= start) & (times <= stop) - temp = {'Start':start, 'Stop':stop} + temp = {'Start': start, 'Stop': stop} for comp in range(nComp): local_times = times[selected] local_values = sim.root.output.solution.unit_001["solution_outlet_comp_%03d" % comp][selected] - + temp[str(comp)] = numpy.trapz(local_values, local_times) / (stop - start) df = df.append(temp, ignore_index=True) return df + def main(): + if not os.path.exists("tmp"): + os.makedirs("tmp") simulation = Cadet(common.common.root) - simulation.filename = "MSSMA_2comp.h5" + simulation.filename = "tmp/MSSMA_2comp.h5" createSimulation(simulation) simulation.save() data = simulation.run() print(data) simulation.load() plotSimulation(simulation) + return simulation + def createSimulation(simulation): root = simulation.root @@ -97,18 +88,18 @@ def createSimulation(simulation): root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0, 0.0] root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' root.input.model.unit_001.velocity = 0.001 - + root.input.model.unit_001.adsorption.is_kinetic = 1 root.input.model.unit_001.adsorption.mssma_ka = [0.0, 1E11, 8E6, 1E11, 8E6] root.input.model.unit_001.adsorption.mssma_kd = [0.0, 6E11, 2E16, 6E11, 2E16] root.input.model.unit_001.adsorption.mssma_lambda = 225 root.input.model.unit_001.adsorption.mssma_nu = [0.0, 10, 25, 20, 50] - root.input.model.unit_001.adsorption.mssma_sigma = [0.0, 48, 66, 48*2, 66*2] + root.input.model.unit_001.adsorption.mssma_sigma = [0.0, 48, 66, 48 * 2, 66 * 2] root.input.model.unit_001.adsorption.mssma_refq = 225 root.input.model.unit_001.adsorption.mssma_refc0 = 520.0 root.input.model.unit_001.adsorption.mssma_rates = [0.0, 0.0, 1e20, 10, 0.0, 0.0, 1e20, 10, 0.0] - root.input.model.unit_001.discretization.nbound = [1, 2,2] + root.input.model.unit_001.discretization.nbound = [1, 2, 2] root.input.model.unit_001.discretization.ncol = 50 root.input.model.unit_001.discretization.npar = 5 @@ -122,12 +113,14 @@ def createSimulation(simulation): root.input.solver.consistent_init_mode = 3 + def plotSimulation(simulation): f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) plotInlet(ax1, simulation) plotOutlet(ax2, simulation) f.tight_layout() - plt.show() + plt.savefig("tmp/MSSMA.png") + def plotInlet(axis, simulation): solution_times = simulation.root.output.solution.solution_times @@ -139,7 +132,7 @@ def plotInlet(axis, simulation): axis.set_title("Inlet") axis.plot(solution_times, inlet_salt, 'b-', label="Salt") axis.set_xlabel('time (s)') - + # Make the y-axis label, ticks and tick labels match the line color. axis.set_ylabel('mMol Salt', color='b') axis.tick_params('y', colors='b') @@ -150,7 +143,6 @@ def plotInlet(axis, simulation): axis2.set_ylabel('mMol Protein', color='r') axis2.tick_params('y', colors='r') - lines, labels = axis.get_legend_handles_labels() lines2, labels2 = axis2.get_legend_handles_labels() axis2.legend(lines + lines2, labels + labels2, loc=0) @@ -164,22 +156,22 @@ def plotOutlet(axis, simulation): outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 data = numpy.vstack([solution_times, outlet_p1]).transpose() - numpy.savetxt('comp2_1.csv', data, delimiter=',') + numpy.savetxt('tmp/comp2_1.csv', data, delimiter=',') data = numpy.vstack([solution_times, outlet_p2]).transpose() - numpy.savetxt('comp2_2.csv', data, delimiter=',') + numpy.savetxt('tmp/comp2_2.csv', data, delimiter=',') data = numpy.vstack([solution_times, outlet_p1 + outlet_p2]).transpose() - numpy.savetxt('comp2_comb.csv', data, delimiter=',') + numpy.savetxt('tmp/comp2_comb.csv', data, delimiter=',') fraction_times = gen_fraction_times(6000, 14000, 10) df = gen_fractions(fraction_times, simulation) - df.to_csv("comp2_fraction.csv", columns=('Start', 'Stop', '1'), index=False) + df.to_csv("tmp/comp2_fraction.csv", columns=('Start', 'Stop', '1'), index=False) axis.set_title("Output") axis.plot(solution_times, outlet_salt, 'b-', label="Salt") axis.set_xlabel('time (s)') - + # Make the y-axis label, ticks and tick labels match the line color. axis.set_ylabel('mMol Salt', color='b') axis.tick_params('y', colors='b') @@ -190,13 +182,20 @@ def plotOutlet(axis, simulation): axis2.set_ylabel('mMol Protein', color='r') axis2.tick_params('y', colors='r') - lines, labels = axis.get_legend_handles_labels() lines2, labels2 = axis2.get_legend_handles_labels() axis2.legend(lines + lines2, labels + labels2, loc=0) +@pytest.mark.slow +def test_MSSMA_2comp(): + sim = main() + assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_001, numpy.ndarray) + assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_000, numpy.ndarray) + + if __name__ == "__main__": import sys + print(sys.version) main() diff --git a/examples/SMB.py b/tests/test_SMB.py similarity index 72% rename from examples/SMB.py rename to tests/test_SMB.py index df13e12..07a7e22 100644 --- a/examples/SMB.py +++ b/tests/test_SMB.py @@ -1,65 +1,67 @@ -#!/usr/bin/env python3.6 +# Basic 4-zone SMB setup -#Basic 4-zone SMB setup - -import numpy as np import math +import os.path -from cadet import Cadet -#Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET/VCPKG/bin/cadet-cli.exe" - -#use to render results import matplotlib.pyplot as plt +import numpy +import pytest + +from cadet import Cadet -#number of columns in a cycle +# number of columns in a cycle cycle_size = 8 -#number of cycles +# number of cycles cycles = 4 -#number of times flows have to be expanded for a 4-zone model -repeat_size = int(cycle_size/4) +# number of times flows have to be expanded for a 4-zone model +repeat_size = int(cycle_size / 4) + def gen_connections(units, cycle_size, size, step, flows, flows_static): temp = [] - connections = list(zip(units, np.roll(units,-1), flows)) - io = np.roll(units, step)[[0, size*2, size-1, size*3-1]] + connections = list(zip(units, numpy.roll(units, -1), flows)) + io = numpy.roll(units, step)[[0, size * 2, size - 1, size * 3 - 1]] ios = list(zip([0, 1, 2, 3], io)) for connection in connections: temp.append([connection[0], connection[1], -1, -1, connection[2]]) - #inputs + # inputs idx = 0 for io in ios[:2]: temp.append([io[0], io[1], -1, -1, flows_static[idx]]) - idx+=1; - #outputs + idx += 1 + # outputs for io in ios[2:]: temp.append([io[1], io[0], -1, -1, flows_static[idx]]) - idx+=1; - return np.array(temp) + idx += 1 + return numpy.array(temp) + def expand_flow(seq, inlet1, inlet2, number): "expand the flows for smb, this is more complex since we need link values" temp = [] - temp.extend([seq[3] + inlet1] * (number-1)) + temp.extend([seq[3] + inlet1] * (number - 1)) temp.append(seq[0]) - temp.extend([seq[0]] * (number-1)) + temp.extend([seq[0]] * (number - 1)) temp.append(seq[1]) - temp.extend([seq[1] + inlet2] * (number-1)) + temp.extend([seq[1] + inlet2] * (number - 1)) temp.append(seq[2]) - temp.extend([seq[2]] * (number-1)) + temp.extend([seq[2]] * (number - 1)) temp.append(seq[3]) return temp + def main(): + if not os.path.exists("tmp"): + os.makedirs("tmp") smb = Cadet() - smb.filename = 'F:/temp/SMB.h5' + smb.filename = 'tmp/SMB.h5' createSimulation(smb) print("Simulated Created") smb.save() @@ -67,6 +69,8 @@ def main(): smb.load() print("Simulation Run") plotSimulation(smb) + return smb + def createSimulation(simulation): simulation.root.input.model.nunits = 4 + cycle_size @@ -76,27 +80,31 @@ def createSimulation(simulation): simulation.root.input.model.solver.max_restarts = 0 simulation.root.input.model.solver.schur_safety = 1e-8 - - #setup connections + # setup connections simulation.root.input.model.connections.nswitches = cycle_size - units = range(4, 4+cycle_size) + units = range(4, 4 + cycle_size) flows = expand_flow([7.66E-07, 7.66E-07, 8.08E-07, 8.08E-07], 0.98e-7, 1.96e-07, repeat_size) - flows_static = np.array([0.98e-7, 1.96e-7, 1.4e-7, 1.54e-7]) + flows_static = numpy.array([0.98e-7, 1.96e-7, 1.4e-7, 1.54e-7]) for i in range(cycle_size): simulation.root.input.model.connections["switch_%03d" % i].section = i - simulation.root.input.model.connections["switch_%03d" % i].connections = gen_connections(units, cycle_size, repeat_size, -i, np.array(list(np.roll(flows, i))), flows_static ) - - #setup inlets + simulation.root.input.model.connections["switch_%03d" % i].connections = gen_connections(units, cycle_size, + repeat_size, -i, + numpy.array(list( + numpy.roll(flows, + i))), + flows_static) + + # setup inlets simulation.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' simulation.root.input.model.unit_000.ncomp = 2 simulation.root.input.model.unit_000.unit_type = 'INLET' for i in range(cycle_size): - #section - simulation.root.input.model.unit_000["sec_%03d" % i].const_coeff = [0.55/180.16, 0.55/180.16] + # section + simulation.root.input.model.unit_000["sec_%03d" % i].const_coeff = [0.55 / 180.16, 0.55 / 180.16] simulation.root.input.model.unit_000["sec_%03d" % i].lin_coeff = [0.0, 0.0] simulation.root.input.model.unit_000["sec_%03d" % i].quad_coeff = [0.0, 0.0] simulation.root.input.model.unit_000["sec_%03d" % i].cube_coeff = [0.0, 0.0] @@ -106,23 +114,22 @@ def createSimulation(simulation): simulation.root.input.model.unit_001.unit_type = 'INLET' for i in range(cycle_size): - #section + # section simulation.root.input.model.unit_001["sec_%03d" % i].const_coeff = [0.0, 0.0] simulation.root.input.model.unit_001["sec_%03d" % i].lin_coeff = [0.0, 0.0] simulation.root.input.model.unit_001["sec_%03d" % i].quad_coeff = [0.0, 0.0] simulation.root.input.model.unit_001["sec_%03d" % i].cube_coeff = [0.0, 0.0] - #create columns + # create columns for unit in range(4, 4 + cycle_size): - simulation.root.input.model["unit_%03d" % unit].unit_type = 'GENERAL_RATE_MODEL' col = simulation.root.input.model["unit_%03d" % unit] col.ncomp = 2 - col.cross_section_area = math.pi * (0.02**2)/4.0 + col.cross_section_area = math.pi * (0.02 ** 2) / 4.0 col.col_dispersion = 3.8148e-20 - col.col_length = 0.25/repeat_size + col.col_length = 0.25 / repeat_size col.col_porosity = 0.83 col.init_c = [0.0, 0.0] col.init_q = [0.0, 0.0] @@ -154,14 +161,14 @@ def createSimulation(simulation): col.discretization.weno.weno_eps = 1e-12 col.discretization.weno.weno_order = 3 - #create outlets + # create outlets simulation.root.input.model.unit_002.ncomp = 2 simulation.root.input.model.unit_002.unit_type = 'OUTLET' simulation.root.input.model.unit_003.ncomp = 2 simulation.root.input.model.unit_003.unit_type = 'OUTLET' - #create output information + # create output information simulation.root.input['return'].write_solution_times = 1 @@ -192,12 +199,13 @@ def createSimulation(simulation): ret.unit_003.write_solution_column_outlet = 1 ret.unit_003.write_solution_flux = 0 ret.unit_003.write_solution_particle = 0 - + simulation.root.input.solver.nthreads = 0 - simulation.root.input.solver.user_solution_times = np.linspace(0, cycles*180*4, 1000*cycle_size*cycles) - simulation.root.input.solver.sections.nsec = cycle_size*cycles - simulation.root.input.solver.sections.section_continuity = [0] * (cycle_size*cycles -1) - simulation.root.input.solver.sections.section_times = [float(i) * 180*4.0/cycle_size for i in range(cycle_size*cycles+1)] + simulation.root.input.solver.user_solution_times = numpy.linspace(0, cycles * 180 * 4, 1000 * cycle_size * cycles) + simulation.root.input.solver.sections.nsec = cycle_size * cycles + simulation.root.input.solver.sections.section_continuity = [0] * (cycle_size * cycles - 1) + simulation.root.input.solver.sections.section_times = [float(i) * 180 * 4.0 / cycle_size for i in + range(cycle_size * cycles + 1)] simulation.root.input.solver.time_integrator.abstol = 1e-10 simulation.root.input.solver.time_integrator.algtol = 1e-10 @@ -215,8 +223,7 @@ def plotSimulation(simulation): r_0 = simulation.root.output.solution.unit_003.solution_outlet_comp_000 r_1 = simulation.root.output.solution.unit_003.solution_outlet_comp_001 - - fig = plt.figure(figsize=[10, 2*10]) + fig = plt.figure(figsize=[10, 2 * 10]) graph = fig.add_subplot(2, 1, 1) graph.set_title("Extract") @@ -224,17 +231,28 @@ def plotSimulation(simulation): graph.plot(solution_times, e_0, 'r', label='1') graph.plot(solution_times, e_1, 'g', label='2') graph.legend() - + graph = fig.add_subplot(2, 1, 2) graph.set_title("Raffinate") graph.plot(solution_times, r_0, 'r', label='1') graph.plot(solution_times, r_1, 'g', label='2') graph.legend() - plt.show() + plt.savefig("tmp/SMB.png") + + +@pytest.mark.slow +def test_SMB(): + from datetime import datetime + start = datetime.now() + sim = main() + print(datetime.now() - start) + assert isinstance(sim.root.output.solution.unit_003.solution_outlet_comp_001, numpy.ndarray) + assert isinstance(sim.root.output.solution.unit_003.solution_outlet_comp_000, numpy.ndarray) if __name__ == "__main__": import sys + print(sys.version) main() diff --git a/tests/test_autodetection.py b/tests/test_autodetection.py new file mode 100644 index 0000000..8401049 --- /dev/null +++ b/tests/test_autodetection.py @@ -0,0 +1,6 @@ +from cadet import Cadet + + +def test_autodetection(): + sim = Cadet() + assert sim.cadet_runner is not None diff --git a/tests/test_duplicate_keys.py b/tests/test_duplicate_keys.py index b148b87..53ed45e 100644 --- a/tests/test_duplicate_keys.py +++ b/tests/test_duplicate_keys.py @@ -1,6 +1,7 @@ -import pytest import tempfile +import pytest + from cadet import Cadet diff --git a/examples/test_lwe.py b/tests/test_lwe.py similarity index 78% rename from examples/test_lwe.py rename to tests/test_lwe.py index e134cf5..bd0acad 100644 --- a/examples/test_lwe.py +++ b/tests/test_lwe.py @@ -1,27 +1,13 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that +# Normally the main function is placed at the bottom of the file, but I have placed it at the top so that # This interface is more like a tutorial and can be read from the top down and any given function # can be drilled into to learn more about it. +import os -#use to render results import matplotlib.pyplot as plt - import numpy from cadet import Cadet -import common - -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" +from tests import common # Helper functions that make it easier to set the values in the HDF5 file @@ -29,14 +15,18 @@ # below match those types. def main(): + if not os.path.exists("tmp"): + os.makedirs("tmp") simulation = Cadet(common.common.root) - simulation.filename = "f:/temp/LWE.h5" + simulation.filename = "tmp/LWE.h5" createSimulation(simulation) simulation.save() simulation.run() simulation.load() plotSimulation(simulation) + return simulation + def createSimulation(simulation): root = simulation.root @@ -81,11 +71,10 @@ def createSimulation(simulation): root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0] root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - #root.input.model.unit_001.velocity = 1 - #root.input.model.unit_001.cross_section_area = 4700.352526439483 + # root.input.model.unit_001.velocity = 1 + # root.input.model.unit_001.cross_section_area = 4700.352526439483 root.input.model.unit_001.velocity = 5.75e-4 - root.input.model.unit_001.adsorption.is_kinetic = 0 root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5] root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0] @@ -104,38 +93,39 @@ def createSimulation(simulation): root.input.solver.sections.nsec = 3 root.input.solver.sections.section_continuity = [0, 0] root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0] - + + def plotSimulation(simulation): f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) plotInlet(ax1, simulation) plotOutlet(ax2, simulation) f.tight_layout() - plt.show() + plt.savefig("tmp/lwe.png") + def plotInlet(axis, simulation): solution_times = simulation.root.output.solution.solution_times inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - #inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 - #inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 + # inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 + # inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 axis.set_title("Inlet") axis.plot(solution_times, inlet_salt, 'b-', label="Salt") axis.set_xlabel('time (s)') - + # Make the y-axis label, ticks and tick labels match the line color. axis.set_ylabel('mMol Salt', color='b') axis.tick_params('y', colors='b') axis2 = axis.twinx() axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - #axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - #axis2.plot(solution_times, inlet_p3, 'k-', label="P3") + # axis2.plot(solution_times, inlet_p2, 'g-', label="P2") + # axis2.plot(solution_times, inlet_p3, 'k-', label="P3") axis2.set_ylabel('mMol Protein', color='r') axis2.tick_params('y', colors='r') - lines, labels = axis.get_legend_handles_labels() lines2, labels2 = axis2.get_legend_handles_labels() axis2.legend(lines + lines2, labels + labels2, loc=0) @@ -146,32 +136,30 @@ def plotOutlet(axis, simulation): outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - #outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 - #outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 + # outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 + # outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 axis.set_title("Output") axis.plot(solution_times, outlet_salt, 'b-', label="Salt") axis.set_xlabel('time (s)') - + # Make the y-axis label, ticks and tick labels match the line color. axis.set_ylabel('mMol Salt', color='b') axis.tick_params('y', colors='b') axis2 = axis.twinx() axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - #axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - #axis2.plot(solution_times, outlet_p3, 'k-', label="P3") + # axis2.plot(solution_times, outlet_p2, 'g-', label="P2") + # axis2.plot(solution_times, outlet_p3, 'k-', label="P3") axis2.set_ylabel('mMol Protein', color='r') axis2.tick_params('y', colors='r') - lines, labels = axis.get_legend_handles_labels() lines2, labels2 = axis2.get_legend_handles_labels() axis2.legend(lines + lines2, labels + labels2, loc=0) -if __name__ == "__main__": - import sys - print(sys.version) - main() - +def test_lwe(): + sim = main() + assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_001, numpy.ndarray) + assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_000, numpy.ndarray)