From 5ee5d7a967e06b1e4b2446283030ba9e3e856973 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Fri, 22 Mar 2024 12:34:23 +0100 Subject: [PATCH] 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;