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

Multiple Grid Cells in `musica`#

As we saw in the last tutorial, gas-phase systems can be solved in musica using a solver and a state. For box models, that’s all we need. For more complex applications (columns, 3D grids, etc.), we’ll need to manage the state of a number of independent well-mixed air masses, or “grid cells”.

NOTE: We use “grid cell” throughout these tutorials to mean a well-mixed air mass. A vertical stack of “grid cells” would be used for a column-model. Collections of vertical stacks of grid cells would be used for a 3D grid or mesh. However, micm does not make assumptions about the shape, size, or arragement of the well-mixed air mass(es) a state represents, so the state can be used to represent air masses of any size/shape in any arrangement.

Why multiple grid cell states?

If we can create as many states as we want and use the same solver on them, maybe you’re wondering why we even need multiple-grid-cell states. The answer, once again, is for performance. In HPC applications, micm can be used to solve hundreds of thousands of independent grid cells simultaneously. To do this efficiently, the data needs to be efficiently arranaged in memory. Multiple-grid-cell states allow us to do this.

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) # This is done to make the arrays more readable
np.set_printoptions(suppress=True) # This is done to make the arrays more readable

2. Defining a System#

In the first tutorial, we used an example configuration file to define our chemical system. This time, we’ll use the musica API to define a chemical system in code. We’ll use this apporach for future tutorials as well.

Let’s set up a system analogous to the simple analytically solvable system we used in the first tutorial

[ ]:

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

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

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

Similar to the last tutorial, we have defined a simple system with:

Three species (A, B, and C) in the gas-phase
Two Arrhenius reactions
- A -> B
- B -> C

We try to keep the python API consistent with the format of the configuration files, and have both structured around the science.

The documentation can be used to clarify important details, like the parameter names and units for rate constant types

[ ]:

mc.Arrhenius?

In the last chapter, we used the musica parser to turn our configuration file into a mechanism. Here, we’ll create the mechanism out of the chemistry objects we’ve defined above

[ ]:

mechanism = mc.Mechanism(
    name="musica_micm_example",
    species=species,
    phases=[gas],
    reactions=[r1, r2]
)
for reaction in mechanism.reactions:
    print(f"[{reaction.type.name}] {reaction.to_equation()}")

3. Creating the Solver#

We create a solver for our mechanism the same way we did in the last tutorial. One minor difference here, is that we specifcy the solver_type. For more information on the types of solvers available, see here. The standard-ordered Rosenbrock solver is actually the default solver_type, so including it here has no real effect.

NOTE: If you’re wondering what “standard order” means, it refers to how the matrix data in the state in arranged in memory. micm is designed to be performant enough to be included in large-scale climate and weather models. To achieve optimal performance we have advanced matrix ordering strategies we can apply to encourage vectorization of the solver instructions. For most Python-based applications, the standard ordering will be sufficient.

[ ]:

solver = musica.MICM(mechanism = mechanism, solver_type = musica.SolverType.rosenbrock_standard_order)

4. Creating the State#

In the last tutorial, we created a single-grid-cell state. This time, let’s double the complexity of our system and create a two-grid-cell state!

[ ]:

num_grid_cells = 2
state = solver.create_state(num_grid_cells)

5. Setting Conditions#

Next, let’s set the conditions for our two-grid-cell state. In the last tutorial, you may have wondered why the temperature, pressure, and species concentrations were in array format. Hopefully, now it’s clear - one element for each grid cell!

[ ]:

state.set_conditions(temperatures=[300, 100], pressures=[101253.3, 11253.3])
state.set_concentrations({
    'A': [5, 20],
    'B': [5, 3],
    'C': [5, 7]
})

7. Solving the System#

This time around, we’ll advance the state for 60 s, collecting intermediate states each second.

[ ]:

concentrations_solved = []
time_step = 10 # s
sim_length = 600 # s
curr_time = 0 # s
while curr_time <= sim_length:
    solver.solve(state, time_step)
    concentrations_solved.append(state.get_concentrations())
    curr_time += time_step

8. Preparing the Results (Advanced; Optional Read)#

We’re going to create a short helper function to put the results into a Pandas DataFrame for easier plotting. The function will extract results for a single grid cell, and we call it twice, once for each of the grid cells in our state. In this simulation, the temperature, pressure, and air density are all constant, so numpy’s repeat() function is used to repeat their respective values for every time step.

[ ]:

def convert_results_single_cell(cell_index):
    concentrations_solved_pd = []
    for i in range(0, sim_length + 1, time_step):
        concentrations_solved_pd.append({species: concentration[cell_index] for species, concentration in concentrations_solved[int(i/time_step)].items()})
    df = pd.DataFrame(concentrations_solved_pd)
    df = df.rename(columns = {'A' : 'CONC.A.mol m-3', 'B' : 'CONC.B.mol m-3', 'C' : 'CONC.C.mol m-3'})
    df['time.s'] = list(map(float, range(0, sim_length + 1, time_step)))
    df['ENV.temperature.K'] = np.repeat(state.get_conditions()['temperature'][cell_index], sim_length/time_step + 1.0)
    df['ENV.pressure.Pa'] = np.repeat(state.get_conditions()['pressure'][cell_index], sim_length/time_step + 1.0)
    df['ENV.air number density.mol m-3'] = np.repeat(state.get_conditions()['air_density'][cell_index], sim_length/time_step + 1.0)
    df = df[['time.s', 'ENV.temperature.K', 'ENV.pressure.Pa', 'ENV.air number density.mol m-3', 'CONC.A.mol m-3', 'CONC.B.mol m-3', 'CONC.C.mol m-3']]
    return df

[ ]:

df_0 = convert_results_single_cell(0)
df_1 = convert_results_single_cell(1)

9. Viewing the Results#

Finally, let’s desplay and plot the DataFrames to show the evolution of both of the independent systems over time.

[ ]:

display(df_0)
display(df_1)
df_0.plot(x='time.s', y=['CONC.A.mol m-3', 'CONC.B.mol m-3', 'CONC.C.mol m-3'], title='Concentration over time (grid cell 0)', ylabel='Concentration (mol m-3)', xlabel='Time (s)')
df_1.plot(x='time.s', y=['CONC.A.mol m-3', 'CONC.B.mol m-3', 'CONC.C.mol m-3'], title='Concentration over time (grid cell 1)', ylabel='Concentration (mol m-3)', xlabel='Time (s)')
plt.show()

10. Going Full Circle#

We’ve now seen the two ways to configure musica in practice: using configuration files, and in-code using the musica API. What if we want to create a configuration code from our in-code mechanism?

Let’s give it a try!

[ ]:

mechanism.export("my_cool_musica_mechanism.json")

Let’s take a look at what it looks like

[ ]:

import json

with open("my_cool_musica_mechanism.json", "r") as f:
    config = json.load(f)

print(json.dumps(config, indent=2))

Why is this useful? The musica library than underpins the musica python package is also used in large 3D models. For those applications, text-based configuration files allow run-time specification of the chemistry system without the need to modify source code and recompile the models. So, you can develop and test your mechanism in Python and when you’re ready, turn it into a configuration file that you can use in a model like CheMPAS-A or CAM-SIMA!

Multiple Grid Cells in musica#