From bf7f26daf65eb3b6e7d4fa94a5a3be59a8556ed0 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Mon, 18 Mar 2024 16:39:43 +0100 Subject: [PATCH 1/5] Update formatting to PEP8 --- cadet/cadet.py | 64 +++++++++++++++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 24 deletions(-) diff --git a/cadet/cadet.py b/cadet/cadet.py index 1a3530f..86322d0 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,10 +120,12 @@ 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): _cadet_runner_class = None _is_file_class = None @@ -144,7 +149,7 @@ def __init__(cls): 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,18 @@ 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.return_information = None self._is_file = None + self.cadet_dll_path = None + self.cadet_create_lwe_path = None + self.cadet_cli_path = None + self._install_path = None @property def is_file(self): @@ -173,8 +183,9 @@ 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() @property def cadet_path(self): @@ -188,7 +199,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) @@ -209,14 +220,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 +238,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 +257,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 +277,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 +288,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 +300,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 +322,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): From dc69ee0bbdbc47187618b8f821e258d1d0c64870 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Schm=C3=B6lder?= Date: Mon, 18 Mar 2024 16:40:00 +0100 Subject: [PATCH 2/5] Add autodetect_cadet functionality from CADET-Process --- cadet/cadet.py | 153 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 145 insertions(+), 8 deletions(-) diff --git a/cadet/cadet.py b/cadet/cadet.py index 86322d0..a4257fe 100644 --- a/cadet/cadet.py +++ b/cadet/cadet.py @@ -127,6 +127,10 @@ def is_dll(value): 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 @@ -141,10 +145,6 @@ 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 @@ -165,12 +165,17 @@ class Cadet(H5, metaclass=CadetMeta): def __init__(self, *data): super().__init__(*data) self._cadet_runner = None - self.return_information = None - self._is_file = 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.cadet_cli_path = None - self._install_path = None + + self.return_information = None @property def is_file(self): @@ -186,6 +191,7 @@ def cadet_runner(self): 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): @@ -209,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) From 209dbaadb1d7b8f02f337e34521f7119c24a1cb9 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Fri, 22 Mar 2024 12:33:40 +0100 Subject: [PATCH 3/5] clean up .gitignore --- .gitignore | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) 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 From 026f18e1c077b88c77c83728b9e38e97b891b5e7 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Fri, 22 Mar 2024 12:34:13 +0100 Subject: [PATCH 4/5] Update test suite fixup! Update test suite --- .github/workflows/pipeline.yml | 39 +++++- setup.cfg | 6 +- setup.py | 7 ++ tests/__init__.py | 0 {examples => tests}/common.py | 6 +- .../test_MSSMA_2comp.py | 71 ++++++----- examples/SMB.py => tests/test_SMB.py | 118 ++++++++++-------- tests/test_autodetection.py | 6 + tests/test_duplicate_keys.py | 3 +- {examples => tests}/test_lwe.py | 68 +++++----- 10 files changed, 186 insertions(+), 138 deletions(-) create mode 100644 tests/__init__.py rename {examples => tests}/common.py (90%) rename examples/MSSMA_2comp.py => tests/test_MSSMA_2comp.py (83%) rename examples/SMB.py => tests/test_SMB.py (72%) create mode 100644 tests/test_autodetection.py rename {examples => tests}/test_lwe.py (78%) 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/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) From 5ee5d7a967e06b1e4b2446283030ba9e3e856973 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Fri, 22 Mar 2024 12:34:23 +0100 Subject: [PATCH 5/5] Clean up repository --- examples/CSTR.py | 483 --------------------- examples/Cadet_Python.py | 266 ------------ examples/Cadet_Python_Scoop.py | 361 --------------- examples/HICWang.py | 172 -------- examples/MSSMA.py | 166 ------- examples/MSSMA2.py | 323 -------------- examples/Pilot_Bypass_noBT_noC.py | 198 --------- examples/Test.py | 182 -------- examples/__pycache__/common.cpython-37.pyc | Bin 1589 -> 0 bytes examples/cadet_json.py | 29 -- examples/cadet_library.py | 291 ------------- examples/dll_test.py | 40 -- examples/fix_h5.py | 10 - examples/library_test.py | 11 - examples/ron.py | 176 -------- 15 files changed, 2708 deletions(-) delete mode 100644 examples/CSTR.py delete mode 100644 examples/Cadet_Python.py delete mode 100644 examples/Cadet_Python_Scoop.py delete mode 100644 examples/HICWang.py delete mode 100644 examples/MSSMA.py delete mode 100644 examples/MSSMA2.py delete mode 100644 examples/Pilot_Bypass_noBT_noC.py delete mode 100644 examples/Test.py delete mode 100644 examples/__pycache__/common.cpython-37.pyc delete mode 100644 examples/cadet_json.py delete mode 100644 examples/cadet_library.py delete mode 100644 examples/dll_test.py delete mode 100644 examples/fix_h5.py delete mode 100644 examples/library_test.py delete mode 100644 examples/ron.py 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 186ce56374baa6b80c6f5fd80a07353386124750..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1589 zcmaJ=%TE+B9G=-_huvi#JOu>t{Xm7F2M9b=a0=hS}-Zc7P>bj92w) zym|FzOgtH{{tL$OV4^2XJezp3oflhToMzka`}Ox~`|U!xT-3z(u=Wn@7c}jMB*~W* zU{3tly^N+I4e5^7)Pb&RSbwU449?hDkVP4swL9#bZNN^Xr8tkWDK4On6c;0ZlpEGy zR~s%xvC?jrVMSrp6F;%nt*{=&vsYn#3hP(m1{8Ki@f=jxkm5P4u(OKih_W}TuygVT zqSGY^jwQ%0CkUQLhQ#7N$LtI6VtR9u#djC)U5aW)%)Sh-pw5#>!LOpc!mmaB#L*{j zB{QGAmSlxElVM+nH&7wWBJnsXruZi6O7X2W|J!j_2An{p4Q;%9EY&ow2f0NPVScO& zb}S{=n0_|({Y7heHZUmWK6OEL`N7u8(#po>{OTrLn_q8y`1twDx5r;+j-?PAJwGR& zyn6X)rt$pO``5oey`8n@7OL9n_+*f)%hXheR11tcX*P))bWlQgkh5IR7fX{M>;wiQ z&OWA~v%w&5c^DL%<{{jpEr;xblEf4W+I!hLQL7jmnK z85m%mV{s@7dc?&pXAqg(4D!A!tW#4{pc;o^m%(k{*;5QVj(-Tc+dM4S!5s9oS%}Z2 zsjtoCiSDergHvPDQkXySOpIZT)6GY92I8cNYg}valsT=}-d9}&HC>SS%Hu^Q808*n zYdGGB)1R&yXBs(GRmnE`@69M#2y315D|)mf~ryh;