run_mmc_sim.py

import os

import matplotlib.pyplot as plt
import numpy

from cadet import Cadet
from tests import common

CADET_PATH = "XXXXXXXXXXXXXX/cadet-cli.exe"

def main(mode="iex"):
    if not os.path.exists("tmp"):
        os.makedirs("tmp")
    simulation = Cadet(common.common.root)
    Cadet.cadet_path = CADET_PATH
    simulation.filename = "tmp/LWE.h5"
    createSimulation(simulation, mode=mode)
    simulation.save()
    results = simulation.run()
    print(results)
    print(results.stderr)
    print(results.stdout)
    simulation.load()

    plotSimulation(simulation)
    return simulation


def createSimulation(simulation, mode="iex"):

    root = simulation.root

    root.input.model.nunits = 3
    root.input.model.unit_001.discretization.use_analytic_jacobian = 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

    if mode == "hic":
        root.input.model.unit_000.sec_000.const_coeff = [2000.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 = [2000.0, 0.0]
        root.input.model.unit_000.sec_001.lin_coeff = [-1.5, 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 = [10.0, 0.0]
        root.input.model.unit_000.sec_002.lin_coeff = [0.0, 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.init_c = [2000.0, 0.0]
        root.input.model.unit_001.init_q = [1200.0, 0.0]
        root.input.model.unit_001.adsorption.is_kinetic = 0
        root.input.model.unit_001.adsorption.mmcnfor_ka = [0.0, 1e-2]
        root.input.model.unit_001.adsorption.mmcnfor_kd = [0.0, 1.0]
        root.input.model.unit_001.adsorption.mmcnfor_ks = [0.0, 0.0007]
        root.input.model.unit_001.adsorption.mmcnfor_kp = [0.0, 0.0]
        root.input.model.unit_001.adsorption.mmcnfor_lambda = 1200.0
        root.input.model.unit_001.adsorption.mmcnfor_nu = [0.0, 0]
        root.input.model.unit_001.adsorption.mmcnfor_sigma = [0.0, 0]
        root.input.model.unit_001.adsorption.mmcnfor_n = [0.0, 1]
        root.input.model.unit_001.adsorption.mmcnfor_s = [0.0, 1]
        root.input.model.unit_001.adsorption.mmcnfor_refc = 1
        root.input.model.unit_001.adsorption.mmcnfor_refq = 1
    elif mode == "iex":
        root.input.model.unit_000.sec_000.const_coeff = [50.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 = [50.0, 0.0]
        root.input.model.unit_000.sec_001.lin_coeff = [0.2, 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 = [400.0, 0.0]
        root.input.model.unit_000.sec_002.lin_coeff = [0.0, 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.init_c = [50.0, 0.0]
        root.input.model.unit_001.init_q = [1200.0, 0.0]
        root.input.model.unit_001.adsorption.is_kinetic = 0
        root.input.model.unit_001.adsorption.mmcnfor_ka = [0.0, 1e-2]
        root.input.model.unit_001.adsorption.mmcnfor_kd = [0.0, 1.0]
        root.input.model.unit_001.adsorption.mmcnfor_ks = [0.0, 0.0]
        root.input.model.unit_001.adsorption.mmcnfor_kp = [0.0, 0.0]
        root.input.model.unit_001.adsorption.mmcnfor_lambda = 1200.0
        root.input.model.unit_001.adsorption.mmcnfor_nu = [0.0, 4]
        root.input.model.unit_001.adsorption.mmcnfor_sigma = [0.0, 11.83]
        root.input.model.unit_001.adsorption.mmcnfor_n = [0.0, 0]
        root.input.model.unit_001.adsorption.mmcnfor_s = [0.0, 0]
        root.input.model.unit_001.adsorption.mmcnfor_refc = 1
        root.input.model.unit_001.adsorption.mmcnfor_refq = 1

    root.input.model.unit_001.adsorption_model = 'MMC_NFOR'
    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-5, 6.9e-5]
    root.input.model.unit_001.ncomp = 2
    root.input.model.unit_001.par_diffusion = [7e-10, 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]
    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 = 5.75e-4

    root.input.model.unit_001.discretization.nbound = [1, 1]
    root.input.model.unit_001.discretization.ncol = 100
    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, 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, 1300.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.savefig("tmp/lwe.png")


def plotInlet(axis, simulation):
    solution_times = simulation.root.output.solution.solution_times

    inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000
    inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001
    # inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002
    # inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003

    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__':
    sim_iex = main(mode="iex")
    sim_hic = main(mode="hic")