diff --git a/examples/CSTR.py b/examples/CSTR.py deleted file mode 100644 index 9b35d96..0000000 --- a/examples/CSTR.py +++ /dev/null @@ -1,483 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -#use to render results -import matplotlib.pyplot as plt - -import types - -#location of the cadet binary -cadet_location = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(): - obj = types.SimpleNamespace() - obj.filename = "CSTR_V1_C1.h5" - obj.volume = 1 - obj.conc = 1 - obj.linear = 0 - obj.flowIn = 1 - obj.flowOut = 1 - obj.porosity = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V1_C2.h5" - obj.volume = 1 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V1_C2_F2_F1.h5" - obj.volume = 1 - obj.conc = 2 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V10_C1_F2_F1.h5" - obj.volume = 10 - obj.conc = 1 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V10_C2.h5" - obj.volume = 10 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C1_F2_F2.h5" - obj.volume = 100 - obj.conc = 1 - obj.flowIn = 2 - obj.flowOut = 2 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C2.h5" - obj.volume = 100 - obj.conc = 2 - obj.flowIn = 1 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_V100_C2_F2_F1.h5" - obj.volume = 100 - obj.conc = 2 - obj.flowIn = 2 - obj.flowOut = 1 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_sens.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_con_down.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - 1e-6 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_con_up.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 + 1e-6 - obj.linear = 0.01 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_lin_down.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_lin_up.h5" - obj.porosity = 0.5 - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 + 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_por_up.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 + 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_por_down.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - 1e-6 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_filter_up.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 + 1e-6 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - - obj.filename = "CSTR_filter_down.h5" - obj.volume = 1 - obj.conc = 2 - obj.linear = 0.01 - obj.porosity = 0.5 - obj.flowIn = 2 - obj.flowOut = 1 - obj.filter = 0.01 - 1e-6 - createSimulation(obj) - runSimulation(obj) - plotSimulation(obj) - -def createSimulation(obj): - with h5py.File(obj.filename, 'w') as f: - createInput(f, obj) - -def createInput(f, obj): - input = f.create_group("input") - createModel(input, obj) - createReturn(input) - createSolver(input) - createSensitivity(input) - -def createModel(input, obj): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - # Part of CADET 3.1 - solver = model.create_group('solver') - set_int(solver, 'GS_TYPE', 1) - set_int(solver, 'MAX_KRYLOV', 0) - set_int(solver, 'MAX_RESTARTS', 0) - set_double(solver, 'SCHUR_SAFETY', 1.0e-8) - - createConnections(model, obj) - createInlet(model, obj) - createColumn(model, obj) - createOutlet(model) - -def createConnections(model, obj): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections, obj) - -def createSwitch(connections, obj): - """Create a switch in the system. In 3.0 these are very limited but in 3.1 they allow arbitrary connections between unit operations. - The format is a set of 4 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID. If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, obj.flowIn, - 1, 2, -1, -1, obj.flowOut]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model, obj): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 1) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [obj.conc,]) - set_double(sec, 'LIN_COEFF', [obj.linear,]) - set_double(sec, 'QUAD_COEFF', [0.0,]) - set_double(sec, 'CUBE_COEFF', [0.0,]) - -def createColumn(model, obj): - column = model.create_group('unit_001') - - set_double(column, 'INIT_C', [1.0,]) - set_double(column, 'INIT_VOLUME', obj.volume) - set_double(column, 'FLOWRATE_FILTER', [0.0,]) - set_int(column, 'NCOMP', 1) - set_string(column, 'UNIT_TYPE', 'CSTR') - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 1) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_BULK', 0) - set_int(column, 'WRITE_SENS_INLET', 1) - set_int(column, 'WRITE_SENS_OUTLET', 1) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_BULK', 0) - set_int(column, 'WRITE_SOLUTION_INLET', 1) - set_int(column, 'WRITE_SOLUTION_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - - - column = ret.create_group('unit_002') - - set_int(column, 'WRITE_SENS_BULK', 0) - set_int(column, 'WRITE_SENS_INLET', 1) - set_int(column, 'WRITE_SENS_OUTLET', 1) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_BULK', 0) - set_int(column, 'WRITE_SOLUTION_INLET', 1) - set_int(column, 'WRITE_SOLUTION_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 1) - set_double(solver, 'USER_SOLUTION_TIMES', range(0, 101)) - set_int(solver, 'CONSISTENT_INIT_MODE', 1) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 1) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0, 100.0]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-8) - set_double(time_integrator, 'ALGTOL', 1e-12) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-6) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def createSensitivity(input): - sensitivity = input.create_group('sensitivity') - set_int(sensitivity, 'NSENS', 4) - set_string(sensitivity, 'SENS_METHOD', 'ad1') - - param = sensitivity.create_group('param_000') - set_int(param, 'SENS_UNIT', 0) - set_string(param, 'SENS_NAME', 'CONST_COEFF') - set_int(param, 'SENS_COMP', 0) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', 0) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_001') - set_int(param, 'SENS_UNIT', 0) - set_string(param, 'SENS_NAME', 'LIN_COEFF') - set_int(param, 'SENS_COMP', 0) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', 0) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_002') - set_int(param, 'SENS_UNIT', 1) - set_string(param, 'SENS_NAME', 'POROSITY') - set_int(param, 'SENS_COMP', -1) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', -1) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - - param = sensitivity.create_group('param_003') - set_int(param, 'SENS_UNIT', 1) - set_string(param, 'SENS_NAME', 'FLOWRATE_FILTER') - set_int(param, 'SENS_COMP', -1) - set_int(param, 'SENS_REACTION', -1) - set_int(param, 'SENS_BOUNDPHASE', -1) - set_int(param, 'SENS_SECTION', -1) - - set_double(param, 'SENS_ABSTOL', 1e-6) - set_double(param, 'SENS_FACTOR', 1) - -def runSimulation(obj): - proc = subprocess.Popen([cadet_location, obj.filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(obj): - with h5py.File(obj.filename, 'r') as h5: - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - f.tight_layout() - plt.show() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet = np.array(h5['/output/solution/unit_001/SOLUTION_INLET_COMP_000'].value) - - - axis.set_title("Inlet") - axis.plot(solution_times, inlet, 'b-', label="Product") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol', color='b') - axis.tick_params('y', colors='b') - - lines, labels = axis.get_legend_handles_labels() - axis.legend(lines, labels, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet = np.array(h5['/output/solution/unit_001/SOLUTION_OUTLET_COMP_000'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet, 'b-', label="Product") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol', color='b') - axis.tick_params('y', colors='b') - - lines, labels = axis.get_legend_handles_labels() - axis.legend(lines, labels, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Cadet_Python.py b/examples/Cadet_Python.py deleted file mode 100644 index 5559b91..0000000 --- a/examples/Cadet_Python.py +++ /dev/null @@ -1,266 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "f:/temp/LWE.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - - #read sensitivity data - - #sensitivity of protein 1 to nu 1 - #s1 = simulation.root.output.sensitivity.param_000.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 2 to nu 1 - #s2 = simulation.root.output.sensitivity.param_000.unit_001.sens_column_outlet_comp_002 - - #sensitivity of protein 1 to nu 2 - #s3 = simulation.root.output.sensitivity.param_001.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 1 to sigma 1 - #s4 = simulation.root.output.sensitivity.param_002.unit_001.sens_column_outlet_comp_001 - - #sensitivity of protein 1 to lambda - #s5 = simulation.root.output.sensitivity.param_003.unit_001.sens_column_outlet_comp_001 - - #Sensitvity of first species to loads of all species (except salt) - #s6 = simulation.root.output.sensitivity.param_004.unit_001.sens_column_outlet_comp_001 - - #Sensitvity of first species to velocity - #s7 = simulation.root.output.sensitivity.param_005.unit_001.sens_column_outlet_comp_001 - - - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - root.input.model.unit_001.discretization.use_analytic_jacobian = 1 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 4 - - root.input.model.unit_000.sec_000.const_coeff = [50.0, 1.0, 1.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [100.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.2, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 5.75e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.37 - root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6] - root.input.model.unit_001.init_c = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.init_q = [1200.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.ncomp = 4 - root.input.model.unit_001.par_diffusion = [7e-9, 6.07e-9, 6.07e-9, 6.07e-9] - root.input.model.unit_001.par_porosity = 0.75 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES' - - #root.input.model.unit_001.velocity = 1 - #root.input.model.unit_001.cross_section_area = 4700.352526439483 - root.input.model.unit_001.velocity = 5.75e-4 - - - root.input.model.unit_001.adsorption.is_kinetic = 0 - root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5, 1.59, 7.7] - root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0, 1000.0, 1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 1200.0 - root.input.model.unit_001.adsorption.sma_nu = [0.0, 4.7, 5.29, 3.7] - root.input.model.unit_001.adsorption.sma_sigma = [0.0, 11.83, 10.6, 10.0] - - root.input.model.unit_001.discretization.nbound = [1, 1, 1, 1] - root.input.model.unit_001.discretization.ncol = 100 - root.input.model.unit_001.discretization.npar = 5 - - root.input.model.unit_002.ncomp = 4 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 1500, 1500) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0] - - #sensitivities - root.input.sensitivity.nsens = 0 - root.input.sensitivity.sens_method = 'AD1' - - #nu1 (sensitivity to protein 1 nu) - #root.input.sensitivity.param_000.sens_abstol = 1e-8 - #root.input.sensitivity.param_000.sens_boundphase = 0 - #root.input.sensitivity.param_000.sens_comp = 1 - #root.input.sensitivity.param_000.sens_factor = 1.0 - #root.input.sensitivity.param_000.sens_name = "SMA_NU" - #root.input.sensitivity.param_000.sens_reaction = -1 - #root.input.sensitivity.param_000.sens_section = -1 - #root.input.sensitivity.param_000.sens_unit = 1 - - #nu2 - #root.input.sensitivity.param_001.sens_abstol = 1e-8 - #root.input.sensitivity.param_001.sens_boundphase = 0 - #root.input.sensitivity.param_001.sens_comp = 2 - #root.input.sensitivity.param_001.sens_factor = 1.0 - #root.input.sensitivity.param_001.sens_name = "SMA_NU" - #root.input.sensitivity.param_001.sens_reaction = -1 - #root.input.sensitivity.param_001.sens_section = -1 - #root.input.sensitivity.param_001.sens_unit = 1 - - #sigma1 - #root.input.sensitivity.param_002.sens_abstol = 1e-8 - #root.input.sensitivity.param_002.sens_boundphase = 0 - #root.input.sensitivity.param_002.sens_comp = 1 - #root.input.sensitivity.param_002.sens_factor = 1.0 - #root.input.sensitivity.param_002.sens_name = "SMA_SIGMA" - #root.input.sensitivity.param_002.sens_reaction = -1 - #root.input.sensitivity.param_002.sens_section = -1 - #root.input.sensitivity.param_002.sens_unit = 1 - - #lambda - #root.input.sensitivity.param_003.sens_abstol = 1e-8 - #root.input.sensitivity.param_003.sens_boundphase = -1 - #root.input.sensitivity.param_003.sens_comp = -1 - #root.input.sensitivity.param_003.sens_factor = 1.0 - #root.input.sensitivity.param_003.sens_name = "SMA_LAMBDA" - #root.input.sensitivity.param_003.sens_reaction = -1 - #root.input.sensitivity.param_003.sens_section = -1 - #root.input.sensitivity.param_003.sens_unit = 1 - - #Sensitvity of first species to loads of all species (except salt) - #root.input.sensitivity.param_004.sens_abstol = [1e-8, 1e-8, 1e-8] - #root.input.sensitivity.param_004.sens_boundphase = [-1, -1, -1] - #root.input.sensitivity.param_004.sens_comp = [1, 2, 3] - #root.input.sensitivity.param_004.sens_factor = [1.0, 1.0, 1.0] - #root.input.sensitivity.param_004.sens_name = ["CONST_COEFF", "CONST_COEFF", "CONST_COEFF"] - #root.input.sensitivity.param_004.sens_reaction = [-1, -1, -1] - #root.input.sensitivity.param_004.sens_section = [0, 0, 0] - #root.input.sensitivity.param_004.sens_unit = [0, 0, 0] - - #Sensitvity of first species to velocity - #root.input.sensitivity.param_005.sens_abstol = 1e-8 - #root.input.sensitivity.param_005.sens_boundphase = -1 - #root.input.sensitivity.param_005.sens_comp = -1 - #root.input.sensitivity.param_005.sens_factor = 1.0 - #root.input.sensitivity.param_005.sens_name = "VELOCITY" - #root.input.sensitivity.param_005.sens_reaction = -1 - #root.input.sensitivity.param_005.sens_section = -1 - #root.input.sensitivity.param_005.sens_unit = 1 - - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 - inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 - outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Cadet_Python_Scoop.py b/examples/Cadet_Python_Scoop.py deleted file mode 100644 index 7c68007..0000000 --- a/examples/Cadet_Python_Scoop.py +++ /dev/null @@ -1,361 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt - -#parallelization -from scoop import futures - -import itertools -import hashlib - -#location of the cadet binary -cadet_location = "C:\\Users\\kosh_000\\cadet_build\\cadet_3\\cadet\\MS_SMKL_RELEASE\\bin\\cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(args): - ka, nu, sigma = args - #filename = "%.2g_%.2g_%.2g.h5" % (ka, nu, sigma) - - filename = hashlib.md5(str(args).encode('utf-8','ignore')).hexdigest() + '.h5' - - createSimulation(filename, ka, nu, sigma) - runSimulation(filename) - plotSimulation(filename) - return "Completed Ka: %s Nu: %s Sigma: %s" % (ka, nu, sigma) - -def createSimulation(filename, ka, nu, sigma): - with h5py.File(filename, 'w') as f: - createInput(f, ka, nu, sigma) - -def createInput(f, ka, nu, sigma): - input = f.create_group("input") - createModel(input, ka, nu, sigma) - createReturn(input) - createSolver(input) - -def createModel(input, ka, nu, sigma): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - # Part of CADET 3.1 - set_int(model, 'GS_TYPE', 1) - set_int(model, 'MAX_KRYLOV', 0) - set_int(model, 'MAX_RESTARTS', 0) - set_double(model, 'SCHUR_SAFETY', 1.0e-8) - - createConnections(model) - createInlet(model) - createColumn(model, ka, nu, sigma) - createOutlet(model) - -def createConnections(model): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections) - -def createSwitch(connections): - """Create a switch in the system. In 3.0 these are very limited but in 3.1 they allow arbitrary connections between unit operations. - The format is a set of 4 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID. If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, 1, 2, -1, -1]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 4) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - # CADET 3.1 - set_double(inlet, 'FLOW', 1.0) - - createLoad(inlet) - createWash(inlet) - createElute(inlet) - -def createLoad(inlet): - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [50.0, 1.0, 1.0, 1.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createWash(inlet): - sec = inlet.create_group('sec_001') - - set_double(sec, 'CONST_COEFF', [50.0, 0.0, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createElute(inlet): - sec = inlet.create_group('sec_002') - - set_double(sec, 'CONST_COEFF', [100.0, 0.0, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.2, 0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0, 0.0]) - -def createColumn(model, ka, nu, sigma): - column = model.create_group('unit_001') - - set_string(column, 'ADSORPTION_MODEL', 'STERIC_MASS_ACTION') - set_double(column, 'COL_DISPERSION', 5.75e-8) - set_double(column, 'COL_LENGTH', 0.014) - set_double(column, 'COL_POROSITY', 0.37) - set_double(column, 'FILM_DIFFUSION', [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6]) - set_double(column, 'INIT_C', [50.0, 0.0, 0.0, 0.0]) - set_double(column, 'INIT_Q', [1200.0, 0.0, 0.0, 0.0]) - set_int(column, 'NCOMP', 4) - set_double(column, 'PAR_DIFFUSION', [7e-10, 6.07e-11, 6.07e-11, 6.07e-11]) - set_double(column, 'PAR_POROSITY', 0.75) - set_double(column, 'PAR_RADIUS', 4.5e-5) - set_double(column, 'PAR_SURFDIFFUSION', [0.0, 0.0, 0.0, 0.0]) - set_string(column, 'UNIT_TYPE', 'GENERAL_RATE_MODEL') - set_double(column, 'VELOCITY', 5.75e-4) - - # CADET 3.1 - set_double(column, 'FLOW', 1.0) - - createAdsorption(column, ka, nu, sigma) - createDiscretization(column) - -def createAdsorption(column, ka, nu, sigma): - ads = column.create_group('adsorption') - - set_int(ads, 'IS_KINETIC', 0) - set_double(ads, 'SMA_KA', [0.0, ka, 1.59, 7.7]) - set_double(ads, 'SMA_KD', [0, 1000, 1000, 1000]) - set_double(ads, 'SMA_LAMBDA', 1200) - set_double(ads, 'SMA_NU', [0.0, nu, 5.29, 3.7]) - set_double(ads, 'SMA_SIGMA', [0, sigma, 10.6, 10.0]) - -def createDiscretization(column): - disc = column.create_group('discretization') - - set_int(disc, 'GS_TYPE', 1) - set_int(disc, 'MAX_KRYLOV', 0) - set_int(disc, 'MAX_RESTARTS', 10) - set_int(disc, 'NBOUND', [1, 1, 1, 1]) - set_int(disc, 'NCOL', 50) - set_int(disc, 'NPAR', 10) - set_string(disc, 'PAR_DISC_TYPE', 'EQUIDISTANT_PAR') - set_double(disc, 'SCHUR_SAFETY', 1.0e-8) - set_int(disc, 'USE_ANALYTIC_JACOBIAN', 1) - - createWENO(disc) - -def createWENO(disc): - weno = disc.create_group("weno") - set_int(weno, 'BOUNDARY_MODEL', 0) - set_double(weno, 'WENO_EPS', 1e-10) - set_int(weno, 'WENO_ORDER', 3) - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 4) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - - # CADET 3.1 - set_double(outlet, 'FLOW', 1.0) - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_COLUMN', 0) - set_int(column, 'WRITE_SENS_COLUMN_INLET', 0) - set_int(column, 'WRITE_SENS_COLUMN_OUTLET', 0) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_COLUMN', 0) - set_int(column, 'WRITE_SOLUTION_COLUMN_INLET', 1) - set_int(column, 'WRITE_SOLUTION_COLUMN_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 1) - set_double(solver, 'USER_SOLUTION_TIMES', range(0, 1500+1)) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 3) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0, 10, 90, 1500]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-8) - set_double(time_integrator, 'ALGTOL', 1e-12) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-6) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def runSimulation(filename): - proc = subprocess.Popen([cadet_location, filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(filename): - with h5py.File(filename, 'r') as h5: - fig = plt.figure(figsize=[20, 10]) - - ax1 = fig.add_subplot(1, 2, 1) - ax2 = fig.add_subplot(1, 2, 2) - - #f, (ax1, ax2) = fig.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - fig.tight_layout() - #plt.show() - plt.savefig(filename.replace('h5', 'png'), dpi=100) - plt.close() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_000'].value) - inlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_001'].value) - inlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_002'].value) - inlet_p3 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_003'].value) - - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_000'].value) - outlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_001'].value) - outlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_002'].value) - outlet_p3 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_003'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - - ka = [30, 35, 40] - nu = [4.6, 4.7, 4.8] - sigma = [11.7, 11.8, 11.9] - - settings = itertools.product(ka, nu, sigma) - - results = futures.map(main, settings) - - for result in results: - print(result) diff --git a/examples/HICWang.py b/examples/HICWang.py deleted file mode 100644 index 53bce64..0000000 --- a/examples/HICWang.py +++ /dev/null @@ -1,172 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -#Cadet.cadet_path = r"C:\Users\kosh_000\cadet_build\CADET-fran\MS_SMKL_DEBUG\bin\cadet-cli.exe" -Cadet.cadet_path = r"C:\Users\kosh_000\cadet_build\CADET-fran\MS_SMKL_RELEASE\bin\cadet-cli.exe" - - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "F:/temp/HICWang.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - - root.input.solver.time_integrator.abstol = 1e-6 - root.input.solver.time_integrator.reltol = 0.0 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 2 - - root.input.model.unit_000.sec_000.const_coeff = [4000.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [4000.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [-4000/15000, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'HICWANG' - root.input.model.unit_001.col_dispersion = 5.6e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.4 - root.input.model.unit_001.film_diffusion = [1.4e-7, 1.4e-7] - root.input.model.unit_001.init_c = [4000.0, 0.0] - root.input.model.unit_001.init_q = [0.0] - root.input.model.unit_001.ncomp = 2 - root.input.model.unit_001.par_diffusion = [4e-12, 4e-12] - root.input.model.unit_001.par_porosity = 0.98 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0] - root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - - #root.input.model.unit_001.velocity = 1 - #root.input.model.unit_001.cross_section_area = 4700.352526439483 - root.input.model.unit_001.velocity = 60/(3600*100) - - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.hicwang_kkin = [1e2, 1.0/0.16] - root.input.model.unit_001.adsorption.hicwang_keq = [0.0, 34] - root.input.model.unit_001.adsorption.hicwang_nu = [1.0, 9.5] - root.input.model.unit_001.adsorption.hicwang_qmax = [1.0, 1.3e-2*1000] - root.input.model.unit_001.adsorption.hicwang_beta0 = 3.6e-2 - root.input.model.unit_001.adsorption.hicwang_beta1 = 1.0/1e3 - - root.input.model.unit_001.discretization.nbound = [0, 1] - root.input.model.unit_001.discretization.ncol = 3 - root.input.model.unit_001.discretization.npar = 1 - root.input.model.unit_001.discretization.use_analytic_jacobian = 0 - - - root.input.model.unit_002.ncomp = 2 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 15000, 15000) - root.input.solver.sections.nsec = 2 - root.input.solver.sections.section_continuity = [0] - root.input.solver.sections.section_times = [0.0, 10.0, 15000.0] - - root.input.solver.time_integrator.init_step_size = 1e-2 - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_column_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_column_inlet_comp_001 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_column_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_column_outlet_comp_001 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() - diff --git a/examples/MSSMA.py b/examples/MSSMA.py deleted file mode 100644 index a6ce018..0000000 --- a/examples/MSSMA.py +++ /dev/null @@ -1,166 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -import numpy - -#use to render results -import matplotlib.pyplot as plt - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/cadet3.1-win7-x64/bin/cadet-cli.exe" - -import common - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "MSSMA.h5" - createSimulation(simulation) - simulation.save() - simulation.run() - simulation.load() - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 2 - - root.input.model.unit_000.sec_000.const_coeff = [100.0, 0.15] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [75, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [75, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.05, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'MULTISTATE_STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 1.5E-7 - root.input.model.unit_001.col_length = 0.25 - root.input.model.unit_001.col_porosity = 0.3 - root.input.model.unit_001.film_diffusion = [2.14E-4, 2.1e-5] - root.input.model.unit_001.init_c = [75, 0.0] - root.input.model.unit_001.init_q = [225, 0.0, 0.0] - root.input.model.unit_001.ncomp = 2 - root.input.model.unit_001.par_diffusion = [4.08E-10, 9.0E-12] - root.input.model.unit_001.par_porosity = 0.4 - root.input.model.unit_001.par_radius = 3.25E-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - root.input.model.unit_001.velocity = 0.001 - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.mssma_ka = [0.0, 1E11, 8E6] - root.input.model.unit_001.adsorption.mssma_kd = [0.0, 6E11, 2E16] - root.input.model.unit_001.adsorption.mssma_lambda = 225 - root.input.model.unit_001.adsorption.mssma_nu = [0.0, 10, 25] - root.input.model.unit_001.adsorption.mssma_sigma = [0.0, 48, 66] - root.input.model.unit_001.adsorption.mssma_refq = 225 - root.input.model.unit_001.adsorption.mssma_refc0 = 520.0 - root.input.model.unit_001.adsorption.mssma_rates = [0.0, 0.0, 1e20, 10, 0.0] - - root.input.model.unit_001.discretization.nbound = [1, 2] - root.input.model.unit_001.discretization.ncol = 50 - root.input.model.unit_001.discretization.npar = 5 - - root.input.model.unit_002.ncomp = 2 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 15000, 1000) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 4500, 6100, 15000] - - root.input.solver.consistent_init_mode = 3 - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - - data = numpy.vstack([solution_times, outlet_p1]).transpose() - numpy.savetxt('comp_1.csv', data, delimiter=',') - - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/MSSMA2.py b/examples/MSSMA2.py deleted file mode 100644 index 36d633a..0000000 --- a/examples/MSSMA2.py +++ /dev/null @@ -1,323 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET_31.pdf in the same directory -# - -#Designed for CADET 3.1 - -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -# The core part of python with CADET is HDF5 and numpy -import h5py -import numpy as np - -#used to run CADET -import subprocess -import io - -#use to render results -import matplotlib.pyplot as plt - -#location of the cadet binary -cadet_location = "C:/Users/kosh_000/cadet_build/CADET/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def set_int(node, nameH5, value): - "set one or more integers in the hdf5 file" - data = np.array(value, dtype="i4") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_double(node, nameH5, value): - "set one or more doubles in the hdf5 file" - data = np.array(value, dtype="f8") - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data, maxshape=tuple(None for i in range(data.ndim)), fillvalue=[0]) - -def set_string(node, nameH5, value): - "set a string value in the hdf5 file" - if isinstance(value, list): - dtype = 'S' + str(len(value[0])+1) - else: - dtype = 'S' + str(len(value)+1) - data = np.array(value, dtype=dtype) - if node.get(nameH5, None) is not None: - del node[nameH5] - node.create_dataset(nameH5, data=data) - - -def main(): - filename = "MSSMA2.h5" - createSimulation(filename) - runSimulation(filename) - plotSimulation(filename) - -def createSimulation(filename): - with h5py.File(filename, 'w') as f: - createInput(f) - -def createInput(f): - input = f.create_group("input") - createModel(input) - createReturn(input) - createSolver(input) - -def createModel(input): - model = input.create_group("model") - set_int(model, 'NUNITS', 3) - - solver = model.create_group('solver') - set_int(solver, 'GS_TYPE', 1) - set_int(solver, 'MAX_KRYLOV', 0) - set_int(solver, 'MAX_RESTARTS', 10) - set_double(solver, 'SCHUR_SAFETY', 1e-8) - - createConnections(model) - createInlet(model) - createColumn(model) - createOutlet(model) - -def createConnections(model): - connections = model.create_group("connections") - set_int(connections, 'NSWITCHES', 1) - createSwitch(connections) - -def createSwitch(connections): - """Create a switch in the system. In 3.1 they allow arbitrary connections between unit operations. - The format is a set of 5 numbers (Source Unit Operation ID, Destination Unit Operation ID, Source Component ID, Destination Component ID and Volumetric flow rate (m^3/s). If the - Source and Destination Component IDs are set to -1 that means connect all to all and this is the most common setup.""" - - # CADET uses 0 based numbers and 3 digits to identify anything there are multiples of like unit operations, sections, switches etc - switch_000 = connections.create_group("switch_000") - - #Connect all of Inlet [0] to Column [1] and Column [1] to Outlet [2] - set_int(switch_000, 'CONNECTIONS', [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0]) - set_int(switch_000, 'SECTION', 0) - -def createInlet(model): - inlet = model.create_group("unit_000") - - set_string(inlet, 'INLET_TYPE', 'PIECEWISE_CUBIC_POLY') - set_int(inlet, 'NCOMP', 3) - set_string(inlet, 'UNIT_TYPE', 'INLET') - - createLoad(inlet) - createWash(inlet) - createElute(inlet) - -def createLoad(inlet): - sec = inlet.create_group('sec_000') - - set_double(sec, 'CONST_COEFF', [92.0, 0.10631294584377825, 0.10631294584377825/2]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createWash(inlet): - sec = inlet.create_group('sec_001') - - set_double(sec, 'CONST_COEFF', [69.97439960989882, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createElute(inlet): - sec = inlet.create_group('sec_002') - - set_double(sec, 'CONST_COEFF', [69.97439960989882, 0.0, 0.0]) - set_double(sec, 'LIN_COEFF', [0.053, 0.0, 0.0]) - set_double(sec, 'QUAD_COEFF', [0.0, 0.0, 0.0]) - set_double(sec, 'CUBE_COEFF', [0.0, 0.0, 0.0]) - -def createColumn(model): - column = model.create_group('unit_001') - - set_string(column, 'ADSORPTION_MODEL', 'MULTISTATE_STERIC_MASS_ACTION') - set_double(column, 'COL_DISPERSION', 1.5E-7) - set_double(column, 'COL_LENGTH', 0.215) - set_double(column, 'COL_POROSITY', 0.33999999999999997) - set_double(column, 'FILM_DIFFUSION', [2.14E-4, 2.1e-5, 2.1e-5]) - set_double(column, 'INIT_C', [69.9743996098988, 0.0, 0.0]) - set_double(column, 'INIT_Q', [223.547, 0.0, 0.0, 0.0, 0.0]) - set_int(column, 'NCOMP', 3) - set_double(column, 'PAR_DIFFUSION', [4.08E-10, 9.0E-12, 9.0e-12]) - set_double(column, 'PAR_POROSITY', 0.39) - set_double(column, 'PAR_RADIUS', 3.25E-5) - set_double(column, 'PAR_SURFDIFFUSION', [0.0, 0.0, 0.0, 0.0, 0.0]) - set_string(column, 'UNIT_TYPE', 'GENERAL_RATE_MODEL') - set_double(column, 'VELOCITY', 0.0011437908496732027) - - createAdsorption(column) - createDiscretization(column) - -def createAdsorption(column): - ads = column.create_group('adsorption') - - set_int(ads, 'IS_KINETIC', 1) - set_double(ads, 'MSSMA_KA', [0.0, 1.0652004307518004E31, 7.724553149425915E26, 1.969122487513422E30, 1.1177522067458229E27]) - set_double(ads, 'MSSMA_KD', [0.0, 5.88452172578919E31, 1.955092026422206E36, 9.923200169245614E32, 8.083909678639826E38]) - set_double(ads, 'MSSMA_LAMBDA', 223.547) - set_double(ads, 'MSSMA_NU', [0.0, 9.618977853171593, 24.75290977103934, 6.058214688510013, 20.20231695871297]) - set_double(ads, 'MSSMA_SIGMA', [0.0, 47.82861669713074, 65.93967947378826, 40.22617141805257, 63.71221340773053 -]) - set_double(ads, 'MSSMA_REFQ', 223.547) - set_double(ads, 'MSSMA_REFC0', 520.0) - set_double(ads, 'MSSMA_RATES', [0.0, 0.0, 9.39710359947847E39, 9.503195767335168, 0.0, 0.0, 5.571477811298548E30, 82427.80960452619, 0.0]) - -def createDiscretization(column): - disc = column.create_group('discretization') - - set_int(disc, 'GS_TYPE', 1) - set_int(disc, 'MAX_KRYLOV', 0) - set_int(disc, 'MAX_RESTARTS', 0) - set_int(disc, 'NBOUND', [1, 2, 2]) - set_int(disc, 'NCOL', 100) - set_int(disc, 'NPAR', 5) - set_string(disc, 'PAR_DISC_TYPE', 'EQUIDISTANT_PAR') - set_double(disc, 'SCHUR_SAFETY', 1.0e-8) - set_int(disc, 'USE_ANALYTIC_JACOBIAN', 1) - - createWENO(disc) - -def createWENO(disc): - weno = disc.create_group("weno") - set_int(weno, 'BOUNDARY_MODEL', 0) - set_double(weno, 'WENO_EPS', 1e-10) - set_int(weno, 'WENO_ORDER', 3) - -def createOutlet(model): - outlet = model.create_group('unit_002') - set_int(outlet, 'NCOMP', 3) - set_string(outlet, 'UNIT_TYPE', 'OUTLET') - -def createReturn(input): - ret = input.create_group('return') - - set_int(ret, 'WRITE_SOLUTION_TIMES', 1) - - createColumnOutput(ret) - -def createColumnOutput(ret): - column = ret.create_group('unit_001') - - set_int(column, 'WRITE_SENS_COLUMN', 0) - set_int(column, 'WRITE_SENS_COLUMN_INLET', 0) - set_int(column, 'WRITE_SENS_COLUMN_OUTLET', 0) - set_int(column, 'WRITE_SENS_FLUX', 0) - set_int(column, 'WRITE_SENS_PARTICLE', 0) - - set_int(column, 'WRITE_SOLUTION_COLUMN', 0) - set_int(column, 'WRITE_SOLUTION_COLUMN_INLET', 1) - set_int(column, 'WRITE_SOLUTION_COLUMN_OUTLET', 1) - set_int(column, 'WRITE_SOLUTION_FLUX', 0) - set_int(column, 'WRITE_SOLUTION_PARTICLE', 0) - -def createSolver(input): - solver = input.create_group("solver") - set_int(solver, 'NTHREADS', 0) - set_int(solver, 'CONSISTENT_INIT_MODE', 3) - set_double(solver, 'USER_SOLUTION_TIMES', np.linspace(0, 14731.2, 1000)) - - createSections(solver) - createTimeIntegrator(solver) - -def createSections(solver): - sections = solver.create_group("sections") - set_int(sections, 'NSEC', 3) - set_int(sections, 'SECTION_CONTINUITY', [0,0]) - set_double(sections, 'SECTION_TIMES', [0.0, 4445.422740524782, 6103.9941690962105, 14731.2]) - -def createTimeIntegrator(solver): - time_integrator = solver.create_group("time_integrator") - set_double(time_integrator, 'ABSTOL', 1e-10) - set_double(time_integrator, 'ALGTOL', 1e-10) - set_double(time_integrator, 'INIT_STEP_SIZE', 1e-8) - set_int(time_integrator, 'MAX_STEPS', 10000) - set_double(time_integrator, 'RELTOL', 1e-6) - -def runSimulation(filename): - proc = subprocess.Popen([cadet_location, filename], bufsize=0, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout, stderr = proc.communicate() - proc.wait() - print("CADET Output") - print(stdout) - print("CADET Errors") - print(stderr) - -def plotSimulation(filename): - with h5py.File(filename, 'r') as h5: - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, h5) - plotOutlet(ax2, h5) - f.tight_layout() - plt.show() - -def plotInlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - inlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_000'].value) - inlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_001'].value) - inlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_INLET_COMP_002'].value) - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'r-', label="P2") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, h5): - solution_times = np.array(h5['/output/solution/SOLUTION_TIMES'].value) - - outlet_salt = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_000'].value) - outlet_p1 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_001'].value) - outlet_p2 = np.array(h5['/output/solution/unit_001/SOLUTION_COLUMN_OUTLET_COMP_002'].value) - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'r-', label="P2") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() diff --git a/examples/Pilot_Bypass_noBT_noC.py b/examples/Pilot_Bypass_noBT_noC.py deleted file mode 100644 index 1a053b8..0000000 --- a/examples/Pilot_Bypass_noBT_noC.py +++ /dev/null @@ -1,198 +0,0 @@ - -import matplotlib.pyplot as plt - -import numpy - -import pandas - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -common = Cadet() -root = common.root - -root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR' -root.input.model.unit_001.discretization.schur_safety = 1.0e-8 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.weno.boundary_model = 0 -root.input.model.unit_001.discretization.weno.weno_eps = 1e-10 -root.input.model.unit_001.discretization.weno.weno_order = 3 -root.input.model.unit_001.discretization.gs_type = 1 -root.input.model.unit_001.discretization.max_krylov = 0 -root.input.model.unit_001.discretization.max_restarts = 10 - -root.input.solver.time_integrator.abstol = 1e-8 -root.input.solver.time_integrator.algtol = 1e-12 -root.input.solver.time_integrator.init_step_size = 1e-6 -root.input.solver.time_integrator.max_steps = 1000000 -root.input.solver.time_integrator.reltol = 1e-6 - -root.input.model.solver.gs_type = 1 -root.input.model.solver.max_krylov = 0 -root.input.model.solver.max_restarts = 10 -root.input.model.solver.schur_safety = 1e-8 - -#CADET 3.1 and CADET-dev flags are in here so that it works with both -#CADET-dev removed column from the name on the inputs and outputs since for many -#operations it no longer makes sense -root.input['return'].write_solution_times = 1 -root.input['return'].split_components_data = 1 -root.input['return'].unit_000.write_sens_bulk = 0 -root.input['return'].unit_000.write_sens_flux = 0 -root.input['return'].unit_000.write_sens_inlet = 0 -root.input['return'].unit_000.write_sens_outlet = 0 -root.input['return'].unit_000.write_sens_particle = 0 -root.input['return'].unit_000.write_solution_bulk = 0 -root.input['return'].unit_000.write_solution_flux = 0 -root.input['return'].unit_000.write_solution_inlet = 1 -root.input['return'].unit_000.write_solution_outlet = 1 -root.input['return'].unit_000.write_solution_particle = 0 -root.input['return'].unit_000.write_sens_column = 0 -root.input['return'].unit_000.write_sens_column_inlet = 0 -root.input['return'].unit_000.write_sens_column_outlet = 0 -root.input['return'].unit_000.write_solution_column = 0 -root.input['return'].unit_000.write_solution_column_inlet = 1 -root.input['return'].unit_000.write_solution_column_outlet = 1 - -root.input['return'].unit_001 = root.input['return'].unit_000 -root.input['return'].unit_002 = root.input['return'].unit_000 -root.input['return'].unit_003 = root.input['return'].unit_000 - -root.input.solver.nthreads = 1 - - -cstr = Cadet() -root = cstr.root - -root.input.model.unit_002.unit_type = 'CSTR' -root.input.model.unit_002.ncomp = 1 -root.input.model.unit_002.nbound =[0] -root.input.model.unit_002.init_c = [0] -root.input.model.unit_002.porosity = 1.0 -root.input.model.unit_002.init_volume = 5e-6 -root.input.model.unit_002.flowrate_filter = 0.0 - -dpfr = Cadet() -root = dpfr.root - -root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES' -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.adsorption_model= 'NONE' -root.input.model.unit_001.init_c = [0] -root.input.model.unit_001.init_q = [] -root.input.model.unit_001.col_dispersion = 4e-8 -root.input.model.unit_001.col_length = 1.0 -root.input.model.unit_001.total_porosity = 1.0 -root.input.model.unit_001.velocity = 1.0 -root.input.model.unit_001.cross_section_area = ((3.75e-5)**2)*numpy.pi - -root.input.model.unit_001.discretization.nbound = [0] -root.input.model.unit_001.discretization.ncol = 50 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.reconstruction = 'WENO' - - -io = Cadet() -root = io.root - -root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' -root.input.model.unit_000.unit_type = 'INLET' -root.input.model.unit_000.ncomp = 1 - -root.input.model.unit_000.sec_000.const_coeff = [5.0] -root.input.model.unit_000.sec_000.lin_coeff = [0] -root.input.model.unit_000.sec_000.quad_coeff = [0] -root.input.model.unit_000.sec_000.cube_coeff = [0] - -root.input.model.unit_000.sec_001.const_coeff = [1000.0] -root.input.model.unit_000.sec_001.lin_coeff = [0] -root.input.model.unit_000.sec_001.quad_coeff = [0] -root.input.model.unit_000.sec_001.cube_coeff = [0] - -root.input.model.unit_000.sec_002.const_coeff = [0] -root.input.model.unit_000.sec_002.lin_coeff = [0] -root.input.model.unit_000.sec_002.quad_coeff = [0] -root.input.model.unit_000.sec_002.cube_coeff = [0] - -root.input.model.unit_003.unit_type = 'OUTLET' -root.input.model.unit_003.ncomp = 1 - -connectivity = Cadet() -root = connectivity.root - -root.input.model.nunits = 4 - -root.input.model.connections.nswitches = 1 -root.input.model.connections.switch_000.section = 0 -root.input.model.connections.switch_000.connections = [0, 2, -1, -1, (100/(60*1e6)), - 2, 1, -1, -1, (100/(60*1e6)), - 1, 3, -1, -1, (100/(60*1e6))] - -root.input.solver.user_solution_times = numpy.linspace(0, 153.87, 1000) -root.input.solver.sections.nsec = 3 -root.input.solver.sections.section_continuity = [0] -root.input.solver.sections.section_times = [0.0, 15.006, 100.0, 153.87] - - -def main(): - sim = Cadet(common.root, dpfr.root, io.root, cstr.root, connectivity.root) - sim.filename = r"F:\jurgen\Pilot_test.h5" - #createSimulation(sim) - sim.save() - sim.run() - sim.load() - plotSimulation(sim) - - #sim = Cadet(common.root, cstr.root, io.root, connectivity.root) - #sim.root.input['return'].unit_001.write_solution_volume = 1 - #sim.root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.5, - # 1, 2, -1, -1, 1.0] - #sim.filename = r"C:\Users\Kohl\Desktop\Cadet\test_file_cstr.h5" - #sim.save() - #sim.run() - #sim.load() - #plotSimulation(sim) - #plotVolume(sim) - - writer = pandas.ExcelWriter(r'F:\jurgen\test_file_cstr.xlsx') - - inputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_000.solution_inlet_comp_000)]) - outputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_outlet_comp_000)]) - #volumes = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Volume', numpy.squeeze(sim.root.output.solution.unit_001.solution_volume))]) - - inputs.to_excel(writer, 'Input', index=False) - outputs.to_excel(writer, 'Output', index=False) - #volumes.to_excel(writer, 'Volume', index=False) - - writer.save() - - -def plotSimulation(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_002.solution_outlet_comp_000) - plt.show() - -def plotVolume(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_001.solution_volume) - plt.show() - -def createSimulation(sim): - root = sim.root - - #SMA Model - root.input.model.unit_001.adsorption_model= 'STERIC_MASS_ACTION' - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.sma_ka = [0,35.5] - root.input.model.unit_001.adsorption.sma_kd = [0,1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 800.0 - root.input.model.unit_001.adsorption.sma_nu = [1,7.0] - root.input.model.unit_001.adsorption.sma_sigma = [0,10.0] - - #Linear isotherm - # root.input.model.unit_001.adsorption_model = 'LINEAR' - - #root.input.model.unit_001.adsorption.is_kinetic = 1 - #root.input.model.unit_001.adsorption.lin_ka = [5e-3] - #root.input.model.unit_001.adsorption.lin_kd = [1e-3] - -if __name__ == "__main__": - main() diff --git a/examples/Test.py b/examples/Test.py deleted file mode 100644 index 01583f4..0000000 --- a/examples/Test.py +++ /dev/null @@ -1,182 +0,0 @@ -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" - -import pandas - -common = Cadet() -root = common.root - -root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR' -root.input.model.unit_001.discretization.schur_safety = 1.0e-8 -root.input.model.unit_001.discretization.use_analytic_jacobian = 1 -root.input.model.unit_001.discretization.weno.boundary_model = 0 -root.input.model.unit_001.discretization.weno.weno_eps = 1e-10 -root.input.model.unit_001.discretization.weno.weno_order = 3 -root.input.model.unit_001.discretization.gs_type = 1 -root.input.model.unit_001.discretization.max_krylov = 0 -root.input.model.unit_001.discretization.max_restarts = 10 - -root.input.solver.time_integrator.abstol = 1e-8 -root.input.solver.time_integrator.algtol = 1e-12 -root.input.solver.time_integrator.init_step_size = 1e-6 -root.input.solver.time_integrator.max_steps = 1000000 -root.input.solver.time_integrator.reltol = 1e-6 - -root.input.model.solver.gs_type = 1 -root.input.model.solver.max_krylov = 0 -root.input.model.solver.max_restarts = 10 -root.input.model.solver.schur_safety = 1e-8 - -#CADET 3.1 and CADET-dev flags are in here so that it works with both -#CADET-dev removed column from the name on the inputs and outputs since for many -#operations it no longer makes sense -root.input['return'].write_solution_times = 1 -root.input['return'].split_components_data = 1 -root.input['return'].unit_000.write_sens_bulk = 0 -root.input['return'].unit_000.write_sens_flux = 0 -root.input['return'].unit_000.write_sens_inlet = 0 -root.input['return'].unit_000.write_sens_outlet = 0 -root.input['return'].unit_000.write_sens_particle = 0 -root.input['return'].unit_000.write_solution_bulk = 0 -root.input['return'].unit_000.write_solution_flux = 0 -root.input['return'].unit_000.write_solution_inlet = 1 -root.input['return'].unit_000.write_solution_outlet = 1 -root.input['return'].unit_000.write_solution_particle = 0 -root.input['return'].unit_000.write_sens_column = 0 -root.input['return'].unit_000.write_sens_column_inlet = 0 -root.input['return'].unit_000.write_sens_column_outlet = 0 -root.input['return'].unit_000.write_solution_column = 0 -root.input['return'].unit_000.write_solution_column_inlet = 1 -root.input['return'].unit_000.write_solution_column_outlet = 1 - -root.input['return'].unit_001 = root.input['return'].unit_000 -root.input['return'].unit_002 = root.input['return'].unit_000 - -root.input.solver.nthreads = 1 - -column_setup = Cadet() -root = column_setup.root - -root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' -root.input.model.unit_001.col_dispersion = 5.75e-8 -root.input.model.unit_001.col_length = 0.014 -root.input.model.unit_001.col_porosity = 0.37 -root.input.model.unit_001.film_diffusion = [6.9e-6] -root.input.model.unit_001.init_c = [0.0] -root.input.model.unit_001.init_q = [0.0] -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.par_diffusion = [7e-10] -root.input.model.unit_001.par_porosity = 0.75 -root.input.model.unit_001.par_radius = 4.5e-5 -root.input.model.unit_001.par_surfdiffusion = [0.0] -root.input.model.unit_001.velocity = 1 -root.input.model.unit_001.cross_section_area = 4700.352526439483 - -root.input.model.unit_001.discretization.nbound = [1] -root.input.model.unit_001.discretization.ncol = 50 -root.input.model.unit_001.discretization.npar = 4 - -cstr = Cadet() -root = cstr.root - -root.input.model.unit_001.unit_type = 'CSTR' -root.input.model.unit_001.ncomp = 1 -root.input.model.unit_001.nbound = 0 -root.input.model.unit_001.init_c = [1.0] -root.input.model.unit_001.porosity = 1.0 -root.input.model.unit_001.init_volume = 10.0 -root.input.model.unit_001.flowrate_filter = 0.0 - -io = Cadet() -root = io.root - -root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' -root.input.model.unit_000.unit_type = 'INLET' -root.input.model.unit_000.ncomp = 1 - -root.input.model.unit_000.sec_000.const_coeff = [1.0] -root.input.model.unit_000.sec_000.lin_coeff = [0.0] -root.input.model.unit_000.sec_000.quad_coeff = [0.0] -root.input.model.unit_000.sec_000.cube_coeff = [0.0] - -root.input.model.unit_000.sec_001.const_coeff = [0.0] -root.input.model.unit_000.sec_001.lin_coeff = [0.0] -root.input.model.unit_000.sec_001.quad_coeff = [0.0] -root.input.model.unit_000.sec_001.cube_coeff = [0.0] - -root.input.model.unit_002.unit_type = 'OUTLET' -root.input.model.unit_002.ncomp = 1 - -connectivity = Cadet() -root = connectivity.root - -root.input.model.nunits = 3 - -root.input.model.connections.nswitches = 1 -root.input.model.connections.switch_000.section = 0 -root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - -root.input.solver.user_solution_times = numpy.linspace(0, 500, 501) -root.input.solver.sections.nsec = 2 -root.input.solver.sections.section_continuity = [0] -root.input.solver.sections.section_times = [0.0, 100.0, 500.0] - - -def main(): - #sim = Cadet(common.root, column_setup.root, io.root, connectivity.root) - #sim.filename = "F:/temp/test_file.h5" - #createSimulation(sim) - #sim.save() - #sim.run() - #sim.load() - #plotSimulation(sim) - - sim = Cadet(common.root, cstr.root, io.root, connectivity.root) - sim.root.input['return'].unit_001.write_solution_volume = 1 - sim.root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.5, - 1, 2, -1, -1, 1.0] - sim.filename = "F:/temp/test_file_cstr.h5" - sim.save() - sim.run() - sim.load() - plotSimulation(sim) - plotVolume(sim) - - writer = pandas.ExcelWriter('F:/temp/test_file_cstr.xlsx') - - inputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_inlet_comp_000)]) - outputs = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Concentration', sim.root.output.solution.unit_002.solution_outlet_comp_000)]) - volumes = pandas.DataFrame.from_items([('Time', sim.root.output.solution.solution_times), ('Volume', numpy.squeeze(sim.root.output.solution.unit_001.solution_volume))]) - - inputs.to_excel(writer, 'Input', index=False) - outputs.to_excel(writer, 'Output', index=False) - volumes.to_excel(writer, 'Volume', index=False) - - writer.save() - - -def plotSimulation(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_002.solution_outlet_comp_000) - plt.show() - -def plotVolume(sim): - plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_001.solution_volume) - plt.show() - -def createSimulation(sim): - root = sim.root - - - root.input.model.unit_001.adsorption_model = 'LINEAR' - - root.input.model.unit_001.adsorption.is_kinetic = 1 - root.input.model.unit_001.adsorption.lin_ka = [5e-3] - root.input.model.unit_001.adsorption.lin_kd = [1e-3] - -if __name__ == "__main__": - main() diff --git a/examples/__pycache__/common.cpython-37.pyc b/examples/__pycache__/common.cpython-37.pyc deleted file mode 100644 index 186ce56..0000000 Binary files a/examples/__pycache__/common.cpython-37.pyc and /dev/null differ diff --git a/examples/cadet_json.py b/examples/cadet_json.py deleted file mode 100644 index 576c129..0000000 --- a/examples/cadet_json.py +++ /dev/null @@ -1,29 +0,0 @@ -import cadet - -data = cadet.H5() -data.filename = r"F:\MCMC\Synthetic\Bypass\Bypass_MCMC_pump_delay\mcmc_refine\mcmc\mcmc.h5" -data.load(paths=["/bounds_change/json",]) - - -print(data) - -print("bounds json", data.root.bounds_change) - - -print("json", data.root.bounds_change.json) - - - -data.load(paths=["/bounds_change",]) - - - -print("bounds_change", data) - - - -data.load() - - - -print("all", data) \ No newline at end of file diff --git a/examples/cadet_library.py b/examples/cadet_library.py deleted file mode 100644 index c2a4517..0000000 --- a/examples/cadet_library.py +++ /dev/null @@ -1,291 +0,0 @@ -#Python 3.5+ -#Depends on addict https://github.com/mewwts/addict -#Depends on h5py, numpy - -from addict import Dict - -import warnings -with warnings.catch_warnings(): - warnings.filterwarnings("ignore",category=FutureWarning) - import h5py -import numpy -import subprocess -import pprint -import copy -import ctypes - -class Cadet(): - - #cadet_path must be set in order for simulations to run - cadet_path = None - cadet_library_path = None - return_information = None - - pp = pprint.PrettyPrinter(indent=4) - - def __init__(self, *data): - self.root = Dict() - self.filename = None - for i in data: - self.root.update(copy.deepcopy(i)) - - def load(self): - if self.filename is not None: - with h5py.File(self.filename, 'r') as h5file: - self.root = Dict(recursively_load(h5file, '/')) - else: - print('Filename must be set before load can be used') - - load_file = load - - def load_memory(self): - self.root = recursively_load_memory(self.memory) - - def save(self): - if self.filename is not None: - with h5py.File(self.filename, 'w') as h5file: - recursively_save(h5file, '/', self.root) - else: - print("Filename must be set before save can be used") - - save_file = save - - def save_memory(self): - self.memory, self.struct_class = recursively_save_memory(self.root) - - def run(self, timeout = None, check=None): - if self.filename is not None: - data = subprocess.run([self.cadet_path, self.filename], timeout = timeout, check=check, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - self.return_information = data - return data - else: - print("Filename must be set before run can be used") - - run_file = run - - def run_memory(self): - cadet_dll = ctypes.cdll.LoadLibrary(self.cadet_library_path) - print(self.memory) - - cadet_dll.cadetGetLibraryVersion.restype = ctypes.c_char_p - print(cadet_dll.cadetGetLibraryVersion().decode()) - - def __str__(self): - temp = [] - temp.append('Filename = %s' % self.filename) - temp.append(self.pp.pformat(self.root)) - return '\n'.join(temp) - - def update(self, merge): - self.root.update(merge.root) - - def __getitem__(self, key): - key = key.lower() - obj = self.root - for i in key.split('/'): - if i: - obj = obj[i] - return obj - - def __setitem__(self, key, value): - key = key.lower() - obj = self.root - parts = key.split('/') - for i in parts[:-1]: - if i: - obj = obj[i] - obj[parts[-1]] = value - -def recursively_load( h5file, path): - - ans = {} - for key, item in h5file[path].items(): - key = key.lower() - if isinstance(item, h5py._hl.dataset.Dataset): - ans[key] = item.value - elif isinstance(item, h5py._hl.group.Group): - ans[key] = recursively_load(h5file, path + key + '/') - return ans - -def recursively_save( h5file, path, dic): - - # argument type checking - if not isinstance(dic, dict): - raise ValueError("must provide a dictionary") - - if not isinstance(path, str): - raise ValueError("path must be a string") - if not isinstance(h5file, h5py._hl.files.File): - raise ValueError("must be an open h5py file") - # save items to the hdf5 file - for key, item in dic.items(): - key = str(key) - if not isinstance(key, str): - raise ValueError("dict keys must be strings to save to hdf5") - #handle int, float, string and ndarray of int32, int64, float64 - if isinstance(item, str): - h5file[path + key.upper()] = numpy.array(item, dtype='S') - - elif isinstance(item, int): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.int32) - - elif isinstance(item, float): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float64: - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float32: - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int32: - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int64: - h5file[path + key.upper()] = item.astype(numpy.int32) - - elif isinstance(item, numpy.ndarray) and item.dtype.kind == 'S': - h5file[path + key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, int) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.int32) - - elif isinstance(item, list) and any(isinstance(i, float) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.int32): - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.float64): - h5file[path + key.upper()] = item - - elif isinstance(item, numpy.float32): - h5file[path + key.upper()] = numpy.array(item, dtype=numpy.float64) - - elif isinstance(item, numpy.bytes_): - h5file[path + key.upper()] = item - - elif isinstance(item, bytes): - h5file[path + key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, str) for i in item): - h5file[path + key.upper()] = numpy.array(item, dtype="S") - - # save dictionaries - elif isinstance(item, dict): - recursively_save(h5file, path + key + '/', item) - # other types cannot be saved and will result in an error - else: - raise ValueError('Cannot save %s/%s key with %s type.' % (path, key.upper(), type(item))) - - -def recursively_save_memory(dic): - - # argument type checking - if not isinstance(dic, dict): - raise ValueError("must provide a dictionary") - - class TempClass(ctypes.Structure): - pass - - # save items to the hdf5 file - _fields_ = [] - values = {} - for key, item in dic.items(): - key = str(key) - if not isinstance(key, str): - raise ValueError("dict keys must be strings to save to hdf5") - #handle int, float, string and ndarray of int32, int64, float64 - if isinstance(item, str): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, int): - _fields_.append( (key.upper(), ctypes.c_int ) ) - - values[key.upper()] = item - - elif isinstance(item, float): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = item - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float64: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.float32: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.float64).ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int32: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, numpy.ndarray) and item.dtype == numpy.int64: - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] =item.astype(numpy.int32).ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, numpy.ndarray) and item.dtype.kind == 'S': - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item.ctypes.data_as(ctypes.c_char_p) - - elif isinstance(item, list) and all(isinstance(i, int) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_int) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.int32).ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - - elif isinstance(item, list) and any(isinstance(i, float) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_double) ) ) - - values[key.upper()] = numpy.array(item, dtype=numpy.float64).ctypes.data_as(ctypes.POINTER(ctypes.c_double)) - - elif isinstance(item, numpy.int32): - _fields_.append( (key.upper(), ctypes.c_int ) ) - - values[key.upper()] = int(item) - - elif isinstance(item, numpy.float64): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = float(item) - - elif isinstance(item, numpy.float32): - _fields_.append( (key.upper(), ctypes.c_double ) ) - - values[key.upper()] = float(item) - - elif isinstance(item, numpy.bytes_): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, bytes): - _fields_.append( (key.upper(), ctypes.c_char_p) ) - - values[key.upper()] = item - - elif isinstance(item, list) and all(isinstance(i, str) for i in item): - _fields_.append( (key.upper(), ctypes.POINTER(ctypes.c_char_p)) ) - - values[key.upper()] = numpy.array(item, dtype="S").ctypes.data_as(ctypes.POINTER(ctypes.c_char_p)) - - # save dictionaries - elif isinstance(item, dict): - temp_recur, TempClass_recur = recursively_save_memory(item) - _fields_.append( (key, TempClass_recur ) ) - values[key] = temp_recur - # other types cannot be saved and will result in an error - else: - raise ValueError('Cannot save %s/%s key with %s type.' % (path, key.upper(), type(item))) - - TempClass._fields_ = _fields_ - temp = TempClass(**values) - #a = None - return temp, TempClass diff --git a/examples/dll_test.py b/examples/dll_test.py deleted file mode 100644 index 380689a..0000000 --- a/examples/dll_test.py +++ /dev/null @@ -1,40 +0,0 @@ -import cadet.cadet -import cadet.cadet_dll -import time -import matplotlib.pyplot as plt - -print('here') - -cadet.Cadet.cadet_path = r"C:\Users\kosh\cadet_build\CADET\VCPKG_42_dll\bin\cadet.dll" - -sim = cadet.Cadet() -sim.filename = r"C:\Users\kosh\cadet_build\CADET\VCPKG_42_dll\bin\LWE.h5" -sim.load() - -sim.root.input['return'].unit_000.write_solution_inlet = 1 -sim.root.input['return'].unit_000.write_solution_bulk = 1 -sim.root.input['return'].unit_000.single_as_multi_port = 0 -sim.root.input['return'].unit_001.write_solution_bulk = 0 - -print('before run') - -sim.run_load() - -print('after run') - -plt.figure(figsize=[15,15]) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_001) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_002) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_outlet_comp_003) -plt.show() - - -plt.figure(figsize=[15,15]) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_001) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_002) -plt.plot(sim.root.output.solution.solution_times, sim.root.output.solution.unit_000.solution_inlet_comp_003) -plt.show() - -print(sim.root.output.solution.unit_000.solution_bulk.shape) - -a = None diff --git a/examples/fix_h5.py b/examples/fix_h5.py deleted file mode 100644 index 25bbaed..0000000 --- a/examples/fix_h5.py +++ /dev/null @@ -1,10 +0,0 @@ -from cadet import Cadet - -fix = Cadet() -fix.filename = r"F:\output\dextran\dextran_pulse.h5" -fix.load() -fix.root.input.model.unit_001.adsorption_model = "LINEAR" -fix.root.input.model.unit_001.adsorption.is_kinetic = 1 -fix.root.input.model.unit_001.adsorption.lin_ka = [0.0] -fix.root.input.model.unit_001.adsorption.lin_kd = [1e3] -fix.save() \ No newline at end of file diff --git a/examples/library_test.py b/examples/library_test.py deleted file mode 100644 index da989ca..0000000 --- a/examples/library_test.py +++ /dev/null @@ -1,11 +0,0 @@ -import cadet_library - -data = cadet_library.Cadet() -data.filename = r"C:\Users\kosh_000\Documents\Visual Studio 2017\Projects\CADETMatch\Examples\Example1\Dextran\dextran_pulse.h5" -data.cadet_library_path = r"C:\Users\kosh_000\cadet_build\CADET-dev\MS_SMKL_RELEASE\bin\cadet.dll" -data.load() - -data.save_memory() - - -data.run_memory() \ No newline at end of file diff --git a/examples/ron.py b/examples/ron.py deleted file mode 100644 index ed2b64a..0000000 --- a/examples/ron.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python3.6 - -# Everything in here is based on CADET3.pdf in the same directory -# - -# Basic Python CADET file based interface compatible with CADET 3.0 and 3.1 -# Some additional fields have been added so that the generated simulations will also -# work in 3.1 and where those differences are I have noted them. -# This whole file follows the CADET pdf documentation. I have broken the system up into many -# functions in order to make it simpler and to make code reuse easier. - -# Normally the main function is placed at the bottom of the file but I have placed it at the top so that -# This interface is more like a tutorial and can be read from the top down and any given function -# can be drilled into to learn more about it. - -#use to render results -import matplotlib.pyplot as plt - -import numpy - -from cadet import Cadet -import common - -Cadet.cadet_path = "C:/Users/kosh_000/cadet_build/CADET-dev/MS_SMKL_RELEASE/bin/cadet-cli.exe" -# Helper functions that make it easier to set the values in the HDF5 file -# In the CADET pdf file each value in the hdf5 file has a type. The functions -# below match those types. - -def main(): - simulation = Cadet(common.common.root) - simulation.filename = "ron.h5" - createSimulation(simulation) - simulation.save() - out = simulation.run() - print(out) - simulation.load() - - plotSimulation(simulation) - -def createSimulation(simulation): - root = simulation.root - - root.input.model.nunits = 3 - root.input.model.unit_001.discretization.use_analytic_jacobian = 1 - - root.input.model.connections.nswitches = 1 - root.input.model.connections.switch_000.section = 0 - root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0, - 1, 2, -1, -1, 1.0] - root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY' - root.input.model.unit_000.unit_type = 'INLET' - root.input.model.unit_000.ncomp = 4 - - root.input.model.unit_000.sec_000.const_coeff = [50.0, 1.0, 1.0, 1.0] - root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_001.const_coeff = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_000.sec_002.const_coeff = [100.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.lin_coeff = [0.2, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0, 0.0, 0.0] - - root.input.model.unit_001.adsorption_model = 'STERIC_MASS_ACTION' - root.input.model.unit_001.col_dispersion = 5.75e-8 - root.input.model.unit_001.col_length = 0.014 - root.input.model.unit_001.col_porosity = 0.37 - root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6, 6.9e-6, 6.9e-6] - root.input.model.unit_001.init_c = [50.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.init_q = [1200.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.ncomp = 4 - root.input.model.unit_001.par_diffusion = [7e-10, 6.07e-11, 6.07e-11, 6.07e-11] - root.input.model.unit_001.par_porosity = 0.75 - root.input.model.unit_001.par_radius = 4.5e-5 - root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0] - root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES' - # root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL' - - root.input.model.unit_001.velocity = 5.75e-4 - - root.input.model.unit_001.adsorption.is_kinetic = 0 - root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5, 1.59, 7.7] - root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0, 1000.0, 1000.0] - root.input.model.unit_001.adsorption.sma_lambda = 1200.0 - root.input.model.unit_001.adsorption.sma_nu = [0.0, 4.7, 5.29, 3.7] - root.input.model.unit_001.adsorption.sma_sigma = [0.0, 11.83, 10.6, 10.0] - - root.input.model.unit_001.discretization.nbound = [1, 1, 1, 1] - root.input.model.unit_001.discretization.ncol = 3 - root.input.model.unit_001.discretization.npar = 1 - - root.input.model.unit_002.ncomp = 4 - root.input.model.unit_002.unit_type = 'OUTLET' - - root.input.solver.user_solution_times = numpy.linspace(0, 1500, 1500) - root.input.solver.sections.nsec = 3 - root.input.solver.sections.section_continuity = [0, 0] - root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0] - - -def plotSimulation(simulation): - f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8]) - plotInlet(ax1, simulation) - plotOutlet(ax2, simulation) - f.tight_layout() - plt.show() - - -def plotInlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000 - inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001 - inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002 - inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003 - - axis.set_title("Inlet") - axis.plot(solution_times, inlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, inlet_p1, 'r-', label="P1") - axis2.plot(solution_times, inlet_p2, 'g-', label="P2") - axis2.plot(solution_times, inlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -def plotOutlet(axis, simulation): - solution_times = simulation.root.output.solution.solution_times - - outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000 - outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001 - outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002 - outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003 - - axis.set_title("Output") - axis.plot(solution_times, outlet_salt, 'b-', label="Salt") - axis.set_xlabel('time (s)') - - # Make the y-axis label, ticks and tick labels match the line color. - axis.set_ylabel('mMol Salt', color='b') - axis.tick_params('y', colors='b') - - axis2 = axis.twinx() - axis2.plot(solution_times, outlet_p1, 'r-', label="P1") - axis2.plot(solution_times, outlet_p2, 'g-', label="P2") - axis2.plot(solution_times, outlet_p3, 'k-', label="P3") - axis2.set_ylabel('mMol Protein', color='r') - axis2.tick_params('y', colors='r') - - - lines, labels = axis.get_legend_handles_labels() - lines2, labels2 = axis2.get_legend_handles_labels() - axis2.legend(lines + lines2, labels + labels2, loc=0) - - -if __name__ == "__main__": - import sys - print(sys.version) - main() -