import marimo
__generated_with = "0.19.0"
app = marimo.App(width="full")
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
# MATE 664 Assignment 01 Demo
## 0D reversible reaction and entropy production
Reaction:
\(
A \;\xrightleftharpoons[\kappa_-]{\kappa_+=1}\; B
\)
We integrate the kinetic ODEs and compute (exact, not linearized):
- \(c_{\mathrm{A}}(t), c_{\mathrm{B}}(t)\)
- \(\Delta\mu_{\mathrm{A}}^{\mathrm{eq}}(t), \Delta\mu_{\mathrm{B}}^{\mathrm{eq}}(t)\)
- entropy production rate \(\dot{\sigma}(t)\)
- cumulative entropy production \(\Delta s(t)=\int\dot{\sigma}\,dt\)
For this question we will assume $c_{\mathrm{T}} = 1$, $k_{+}=1$ and $k_{B}T = 1$
## How to run the notebook
This demo is run as a web-assembly marimo notebook.
- Click the run button at the bottom
to run all the cells.
- We have provided you an interactive UI 
- For the demonstration, you can simply drag the sliders to change the values for the initial values for the problem. If you'd like to play with the code, simply change the python cells and rerun all of them!
- To learn about interactive marimo notebooks, see [Marimo documentation](https://docs.marimo.io)
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
## Simulation Panel Will Show Below
""")
return
@app.cell(hide_code=True)
def _(
cA_0_slider,
cA_eq_slider,
evolve_system,
mo,
plt,
print_system_info,
switch,
):
# marimo plotting cell (all symbols prefixed with _)
_t, _out_exact, _out_lin = evolve_system()
_show_linear = switch.value
# unpack exact results
_cA_ex = _out_exact["cA"]
_cB_ex = _out_exact["cB"]
_dmuA_ex = _out_exact["dmuA"]
_dmuB_ex = _out_exact["dmuB"]
_dsdt_ex = _out_exact["dsdt"]
_s_ex = _out_exact["sol"][:, 2] # ODE-integrated entropy
# unpack linear results
_cA_li = _out_lin["cA"]
_cB_li = _out_lin["cB"]
_dmuA_li = _out_lin["dmuA"]
_dmuB_li = _out_lin["dmuB"]
_dsdt_li = _out_lin["dsdt"]
_s_li = _out_lin["sol"][:, 2]
plt.figure(figsize=(8, 8))
# --- (1) concentration profile ---
plt.subplot(2, 2, 1)
plt.title("Concentration profile")
plt.plot(
_t, _cA_ex, color="C0", linestyle="-", label=r"$c_{\mathrm{A}}$ (exact)"
)
plt.plot(
_t, _cB_ex, color="C1", linestyle="-", label=r"$c_{\mathrm{B}}$ (exact)"
)
if _show_linear:
plt.plot(
_t,
_cA_li,
color="C0",
linestyle="--",
label=r"$c_{\mathrm{A}}$ (linear)",
)
plt.plot(
_t,
_cB_li,
color="C1",
linestyle="--",
label=r"$c_{\mathrm{B}}$ (linear)",
)
plt.xlabel(r"$t$ (abs unit)")
plt.ylabel(r"Concentration / $c_{\mathrm{T}}$")
plt.legend()
# --- (2) chemical potential deviation ---
plt.subplot(2, 2, 2)
plt.title("Chemical potential deviation")
plt.plot(
_t,
_dmuA_ex,
color="C0",
linestyle="-",
label=r"$\Delta\mu_{\mathrm{A}}^{\mathrm{eq}}$ (exact)",
)
plt.plot(
_t,
_dmuB_ex,
color="C1",
linestyle="-",
label=r"$\Delta\mu_{\mathrm{B}}^{\mathrm{eq}}$ (exact)",
)
if _show_linear:
plt.plot(
_t,
_dmuA_li,
color="C0",
linestyle="--",
label=r"$\Delta\mu_{\mathrm{A}}^{\mathrm{eq}}$ (linear)",
)
plt.plot(
_t,
_dmuB_li,
color="C1",
linestyle="--",
label=r"$\Delta\mu_{\mathrm{B}}^{\mathrm{eq}}$ (linear)",
)
plt.xlabel(r"$t$ (abs unit)")
plt.ylabel(r"$\Delta\mu$ (abs unit)")
plt.legend()
# --- (3) entropy production rate ---
plt.subplot(2, 2, 3)
plt.title(r"Entropy production rate $\dot{\sigma}(t)$")
plt.plot(
_t, _dsdt_ex, color="C2", linestyle="-", label=r"$\dot{\sigma}$ (exact)"
)
if _show_linear:
plt.plot(
_t,
_dsdt_li,
color="C2",
linestyle="--",
label=r"$\dot{\sigma}$ (linear)",
)
plt.xlabel(r"$t$ (abs unit)")
plt.ylabel(r"$\dot{\sigma}$ (abs unit)")
plt.legend()
# --- (4) total entropy change ---
plt.subplot(2, 2, 4)
plt.title(r"Total entropy change $\Delta s(t)$")
plt.plot(_t, _s_ex, color="C3", linestyle="-", label=r"$\Delta s$ (exact)")
if _show_linear:
plt.plot(
_t, _s_li, color="C3", linestyle="--", label=r"$\Delta s$ (linear)"
)
plt.xlabel(r"$t$ (abs unit)")
plt.ylabel(r"$\Delta s$ (abs unit)")
plt.legend()
plt.tight_layout()
# _ax = mo.mpl.interactive(plt.gca())
_ax = plt.gca()
mo.hstack(
[
mo.vstack(
[
mo.md("Conditions"),
cA_eq_slider,
cA_0_slider,
switch,
print_system_info(),
]
),
_ax,
],
widths=[1, 2],
)
return
@app.cell(hide_code=True)
def _():
import marimo as mo
return (mo,)
@app.cell(hide_code=True)
def _(cA_0_slider, cA_eq_slider, mo, np, odeint):
def Delta_chem_pot_eq(c, ceq, kBT=1):
return kBT * np.log(c / ceq)
def Delta_chem_pot_lin(c, ceq, kBT=1):
# linearization of kBT*ln(c/ceq) about c=ceq
# ln(1+u) ~ u with u = (c-ceq)/ceq
return kBT * (c - ceq) / ceq
def system_cond():
c_T = 1
cA_0 = cA_0_slider.value
cB_0 = 1 - cA_0
cA_eq = cA_eq_slider.value
cB_eq = c_T - cA_eq
K = cB_eq / cA_eq
r_pos = 1
r_neg = r_pos / K
# cA_eq = 1 / (1 + K) * c_T
#cB_eq = 1 - cA_eq
return cA_0, cB_0, K, r_pos, r_neg, cA_eq, cB_eq
def print_system_info():
"""
Returns a markdown block containing initial system information
based on the current slider values.
"""
cA_0, cB_0, K, r_pos, r_neg, cA_eq, cB_eq = system_cond()
info_md = mo.md(
f"""
Current System Conditions:
- Init conc. of A ($c_{{A}}(t=0)$): `{cA_0:.3f}`
- Init conc. of B ($c_{{B}}(t=0)$): `{cB_0:.3f}`
- Eq. concentration of A ($c_{{A}}^{{eq}}$): `{cA_eq:.3f}`
- Eq. concentration of B ($c_{{B}}^{{eq}}$): `{cB_eq:.3f}`
- Eq. constant ($K$): `{K:.3f}`
- Forward rate ($\\kappa_+$): `{r_pos:.3f}`
- Reverse rate ($\\kappa_-$): `{r_neg:.3f}`
"""
)
return info_md
def evolve_system(tmax=3.5):
t = np.linspace(0, tmax, 1024)
dt = t[1] - t[0]
cA_0, cB_0, K, r_pos, r_neg, cA_eq, cB_eq = system_cond()
def _exact(y, t):
cA, cB, s = y
dcadt = -r_pos * cA + r_neg * cB
dcbdt = +r_pos * cA - r_neg * cB
delta_mu_a = Delta_chem_pot_eq(cA, cA_eq, kBT=1)
delta_mu_b = Delta_chem_pot_eq(cB, cB_eq, kBT=1)
dsdt = -dcadt * delta_mu_a - dcbdt * delta_mu_b
ret = [dcadt, dcbdt, dsdt]
return ret
def _linear(y, t):
cA, cB, s = y
# dcadt = -r_pos * cA + r_neg * cB
# dcbdt = +r_pos * cA - r_neg * cB
delta_mu_a = Delta_chem_pot_lin(cA, cA_eq, kBT=1)
delta_mu_b = Delta_chem_pot_lin(cB, cB_eq, kBT=1)
coeff = r_pos * r_neg / (r_pos + r_neg)
dcadt = coeff * (-delta_mu_a + delta_mu_b)
dcbdt = coeff * (delta_mu_a - delta_mu_b)
dsdt = -dcadt * delta_mu_a - dcbdt * delta_mu_b
return [dcadt, dcbdt, dsdt]
y0 = [cA_0, cB_0, 0]
sol_exact = odeint(_exact, y0, t)
sol_lin = odeint(_linear, y0, t)
# --- append derived outputs (aligned with t) ---
# concentrations
cA_ex, cB_ex = sol_exact[:, 0], sol_exact[:, 1]
cA_li, cB_li = sol_lin[:, 0], sol_lin[:, 1]
# chemical potentials (exact vs linear definitions)
dmuA_ex = Delta_chem_pot_eq(cA_ex, cA_eq, kBT=1)
dmuB_ex = Delta_chem_pot_eq(cB_ex, cB_eq, kBT=1)
dmuA_li = Delta_chem_pot_lin(cA_li, cA_eq, kBT=1)
dmuB_li = Delta_chem_pot_lin(cB_li, cB_eq, kBT=1)
# reconstruct dc/dt for both models (vectorized, same length as t)
dcA_ex = -r_pos * cA_ex + r_neg * cB_ex
dcB_ex = +r_pos * cA_ex - r_neg * cB_ex
coeff = r_pos * r_neg / (r_pos + r_neg)
dcA_li = coeff * (-dmuA_li + dmuB_li)
dcB_li = coeff * ( dmuA_li - dmuB_li)
# entropy production rate (aligned with t)
dsdt_ex = -(dcA_ex * dmuA_ex + dcB_ex * dmuB_ex)
dsdt_li = -(dcA_li * dmuA_li + dcB_li * dmuB_li)
# # cumulative entropy via trapezoid (optional, also aligned with t)
# s_ex = np.concatenate([[0.0], np.cumsum(0.5 * (dsdt_ex[1:] + dsdt_ex[:-1]) * dt)])
# s_li = np.concatenate([[0.0], np.cumsum(0.5 * (dsdt_li[1:] + dsdt_li[:-1]) * dt)])
# pack results as dicts for clarity
out_exact = {
"sol": sol_exact, # columns: cA, cB, s(from ODE)
"cA": cA_ex, "cB": cB_ex,
"dmuA": dmuA_ex, "dmuB": dmuB_ex,
"dsdt": dsdt_ex,
}
out_lin = {
"sol": sol_lin,
"cA": cA_li, "cB": cB_li,
"dmuA": dmuA_li, "dmuB": dmuB_li,
"dsdt": dsdt_li,
}
return t, out_exact, out_lin
return evolve_system, print_system_info
@app.cell(hide_code=True)
def _(mo):
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#K_slider = mo.ui.slider(start=0.1, stop=10, step=0.1, value=1.0, label="$K$", show_value=True)
cA_eq_slider = mo.ui.slider(start=0, stop=1.0, step=0.05, value=0.5, label="Equilibrium $c_{A}^{eq}$", show_value=True)
cA_0_slider = mo.ui.slider(start=0, stop=1.0, step=0.05, value=0.55, label="Initial $c_{A}(t=t0)$", show_value=True)
switch = mo.ui.switch(value=False, label="Show linear solution")
return cA_0_slider, cA_eq_slider, np, odeint, plt, switch
if __name__ == "__main__":
app.run()