This is part 4 of a tutorial series. We recommend following them in order, starting with Part 0: Welcome to ``musica` <0.%20Welcome%20to%20MUSICA.ipynb>`__.

Multiple Mechanisms in the Same Host Application#

In the previous tutorials, we ran a single chemical mechanism at a time. musica also lets you create multiple independent solver instances — each configured for a different mechanism — and run them simultaneously in the same host application.

This is useful when:

  • Multi-domain models where different regions of the atmosphere require different chemistry (e.g., urban vs. remote)

  • Coupled models where gas-phase chemistry and aerosol chemistry are handled by separate solvers

  • Mechanism comparison studies where you want to run two mechanisms side-by-side under the same conditions

In this tutorial, we’ll define two distinct chemical mechanisms, create an independent solver and state for each, and advance both through time in a single simulation loop.

1. Importing Libraries#

Below is a list of the required libraries for this tutorial:

[ ]:

import musica
import musica.mechanism_configuration as mc
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
pd.set_option('display.float_format', str)  # avoid scientific notation
np.set_printoptions(suppress=True)  # avoid scientific notation

2. Defining Mechanism 1#

Our first mechanism is the simple three-species system we have used throughout the tutorials. It has three species (A, B, and C) and two Arrhenius reactions:

  • A → B

  • B → C

[ ]:

# Mechanism 1: A simple A -> B -> C system
A = mc.Species(name="A")
B = mc.Species(name="B")
C = mc.Species(name="C")

species_1 = [A, B, C]
gas_1 = mc.Phase(name="gas", species=species_1)

r1 = mc.Arrhenius(
    name="A_to_B",
    A=4.0e-3,
    C=50,
    reactants=[A],
    products=[B],
    gas_phase=gas_1
)

r2 = mc.Arrhenius(
    name="B_to_C",
    A=4.0e-3,
    C=50,
    reactants=[B],
    products=[C],
    gas_phase=gas_1
)

mechanism_1 = mc.Mechanism(
    name="mechanism_1",
    species=species_1,
    phases=[gas_1],
    reactions=[r1, r2]
)
print("Mechanism 1 reactions:")
for reaction in mechanism_1.reactions:
    print(f"  [{reaction.type.name}] {reaction.to_equation()}")

3. Defining Mechanism 2#

Our second mechanism uses a different set of species (D, E, and F) with its own reactions:

  • D → E + F (a fast reaction)

  • E + F → D (a slow back-reaction)

This represents a reversible chemical system that will evolve toward equilibrium.

[ ]:

# Mechanism 2: A reversible D <=> E + F system
D = mc.Species(name="D")
E = mc.Species(name="E")
F = mc.Species(name="F")

species_2 = [D, E, F]
gas_2 = mc.Phase(name="gas", species=species_2)

r3 = mc.Arrhenius(
    name="D_to_EF",
    A=2.0e-2,
    C=100,
    reactants=[D],
    products=[E, F],
    gas_phase=gas_2
)

r4 = mc.Arrhenius(
    name="EF_to_D",
    A=1.0e-4,
    C=0,
    reactants=[E, F],
    products=[D],
    gas_phase=gas_2
)

mechanism_2 = mc.Mechanism(
    name="mechanism_2",
    species=species_2,
    phases=[gas_2],
    reactions=[r3, r4]
)
print("Mechanism 2 reactions:")
for reaction in mechanism_2.reactions:
    print(f"  [{reaction.type.name}] {reaction.to_equation()}")

4. Creating Independent Solvers and States#

Now we create a separate MICM solver instance for each mechanism. Because each solver is constructed from its own mechanism, they are completely independent. This means:

  • solver_1 only knows about species A, B, and C

  • solver_2 only knows about species D, E, and F

Each solver also gets its own state object to hold the current concentrations and environmental conditions.

[ ]:

# Create solver 1 for mechanism 1
solver_1 = musica.MICM(
    mechanism=mechanism_1,
    solver_type=musica.SolverType.rosenbrock_standard_order
)
state_1 = solver_1.create_state()

# Create solver 2 for mechanism 2
solver_2 = musica.MICM(
    mechanism=mechanism_2,
    solver_type=musica.SolverType.rosenbrock_standard_order
)
state_2 = solver_2.create_state()

print("Solver 1 species:", list(state_1.get_species_ordering().keys()))
print("Solver 2 species:", list(state_2.get_species_ordering().keys()))

5. Setting Initial Conditions#

Each state is configured independently. Mechanism 1 starts with all A and no products. Mechanism 2 starts with all D and no products.

[ ]:

temperature = 300.0  # K
pressure = 101253.3  # Pa

# Initial conditions for mechanism 1
state_1.set_conditions(temperatures=[temperature], pressures=[pressure])
state_1.set_concentrations({'A': [1.0], 'B': [0.0], 'C': [0.0]})

# Initial conditions for mechanism 2
state_2.set_conditions(temperatures=[temperature], pressures=[pressure])
state_2.set_concentrations({'D': [1.0], 'E': [0.0], 'F': [0.0]})

print("Initial concentrations (Mechanism 1):", state_1.get_concentrations())
print("Initial concentrations (Mechanism 2):", state_2.get_concentrations())

6. Running Both Mechanisms Simultaneously#

The key step: we advance both solvers inside the same time loop. Each call to solver.solve() advances its own independent state. The two mechanisms do not interact, but they both progress through the same simulated time.

[ ]:

time_step = 10    # s
sim_length = 600  # s
curr_time = 0     # s

concentrations_1 = []
concentrations_2 = []

concentrations_1.append(state_1.get_concentrations())
concentrations_2.append(state_2.get_concentrations())

while curr_time < sim_length:
    solver_1.solve(state_1, time_step)
    solver_2.solve(state_2, time_step)
    concentrations_1.append(state_1.get_concentrations())
    concentrations_2.append(state_2.get_concentrations())
    curr_time += time_step

print(f"Solved {int(sim_length / time_step)} time steps for both mechanisms.")

7. Visualizing the Results#

Let’s plot the concentration evolution for both mechanisms side-by-side.

[ ]:

times = list(range(0, sim_length + 1, time_step))

# Build DataFrames for each mechanism
def build_df(concentrations, time_list):
    rows = []
    for t, conc in zip(time_list, concentrations):
        row = {'time (s)': t}
        for species, values in conc.items():
            row[species] = values[0]
        rows.append(row)
    return pd.DataFrame(rows)

df_1 = build_df(concentrations_1, times)
df_2 = build_df(concentrations_2, times)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot Mechanism 1
for col in ['A', 'B', 'C']:
    axes[0].plot(df_1['time (s)'], df_1[col], label=col)
axes[0].set_title('Mechanism 1: A → B → C')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Concentration (mol m⁻³)')
axes[0].legend()
axes[0].grid(True)

# Plot Mechanism 2
for col in ['D', 'E', 'F']:
    axes[1].plot(df_2['time (s)'], df_2[col], label=col)
axes[1].set_title('Mechanism 2: D ⇌ E + F')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Concentration (mol m⁻³)')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

8. Key Takeaways#

This tutorial has shown how to:

  1. Define multiple independent mechanisms using the musica API

  2. Create separate solvers — one per mechanism — that coexist in the same Python session

  3. Advance both mechanisms in lock-step through the same simulation time loop

This pattern applies any time you need to couple independent chemistry modules — gas-phase, aerosol, or domain-specific — in one simulation.

Next, see the Local Parallelization Tutorial, which shows how to distribute multiple independent solvers across CPU cores using Dask.