MUSICA Julia API#

The Julia API provides access to MUSICA’s MICM chemical kinetics solver through the Musica.jl package, which uses CxxWrap.jl to interface with the underlying C++ library.

Installation#

Prerequisites#

  • Julia 1.10 or 1.11

  • CMake 3.24 or later

  • A C++ compiler with C++20 support

Building from Source#

  1. Clone and build MUSICA with Julia support:

git clone https://github.com/NCAR/musica.git
cd musica
cmake -S . -B build -D MUSICA_ENABLE_JULIA=ON -D CMAKE_BUILD_TYPE=Release
cmake --build build

  1. Install the Julia package dependencies:

cd julia
julia --project=. -e 'using Pkg; Pkg.instantiate()'

Quick Start#

The following example runs a single-cell chemical solve using a Rosenbrock solver:

using Musica

println("MUSICA version: ", Musica.get_version())

micm = MICM(config_path = "path/to/config")
state = create_state(micm)

set_conditions!(state, temperatures = 298.0, pressures = 101325.0)
set_concentrations!(state, Dict{String,Any}("A" => 1.0, "B" => 0.0))
set_user_defined_rate_parameters!(state, Dict{String,Any}("USER.k1" => 0.001))

result = solve!(micm, state, 60.0)  # integrate 60 seconds
println("Solver state: ", result.state)   # Converged
println("Steps taken: ", result.stats.number_of_steps)

concs = get_concentrations(state)
println("Final [A]: ", concs["A"][1])

API Reference#

Core#

get_version() → String#

Returns the version string of the MUSICA library.

println(Musica.get_version())  # e.g. "0.15.0"

Constants#

The following physical constants are exported from the Musica module:

Name

Value

Description

AVOGADRO

6.02214076 × 10²³ mol⁻¹

Avogadro’s number

BOLTZMANN

1.380649 × 10⁻²³ J K⁻¹

Boltzmann constant

GAS_CONSTANT

AVOGADRO * BOLTZMANN

Universal gas constant (J K⁻¹ mol⁻¹)

Types#

SolverType#

type SolverType#

Enum controlling which MICM solver backend is used. Pass to the MICM constructor.

Value

Int

Description

Rosenbrock

1

Vector-ordered Rosenbrock solver

RosenbrockStandardOrder

2

Standard-ordered Rosenbrock solver (default)

BackwardEuler

3

Vector-ordered Backward Euler solver

BackwardEulerStandardOrder

4

Standard-ordered Backward Euler solver

CudaRosenbrock

5

GPU Rosenbrock solver (requires CUDA build)

Conditions#

type Conditions#

Environmental conditions for a single grid cell.

Fields

  • temperature::Float64 — Temperature in Kelvin

  • pressure::Float64 — Pressure in Pascals

  • air_density::Float64 — Air number density in mol m⁻³

Constructor

Conditions(; temperature=0.0, pressure=0.0, air_density=nothing)

If air_density is not provided and both temperature and pressure are positive, air density is calculated from the Ideal Gas Law: p / (R * T).

c = Conditions(temperature = 298.0, pressure = 101325.0)
# air_density is computed automatically

c2 = Conditions(temperature = 300.0, pressure = 100000.0, air_density = 42.0)

RosenbrockSolverParameters#

type RosenbrockSolverParameters#

Tuning parameters for the Rosenbrock solver family.

Fields

  • relative_tolerance::Float64 — Relative tolerance (default: 1e-6)

  • absolute_tolerances::Union{Vector{Float64}, Nothing} — Per-species absolute tolerances; nothing uses solver defaults

  • h_min::Float64 — Minimum step size in seconds (default: 0.0)

  • h_max::Float64 — Maximum step size in seconds (default: 0.0)

  • h_start::Float64 — Initial step size in seconds (default: 0.0)

  • max_number_of_steps::Int — Maximum number of internal steps (default: 1000)

params = RosenbrockSolverParameters(
    relative_tolerance = 1e-8,
    h_start = 1e-5,
    max_number_of_steps = 500,
)
micm = MICM(config_path = path, solver_parameters = params)

BackwardEulerSolverParameters#

type BackwardEulerSolverParameters#

Tuning parameters for the Backward Euler solver family.

Fields

  • relative_tolerance::Float64 — Relative tolerance (default: 1e-6)

  • absolute_tolerances::Union{Vector{Float64}, Nothing} — Per-species absolute tolerances; nothing uses solver defaults

  • max_number_of_steps::Int — Maximum number of internal steps (default: 11)

  • time_step_reductions::Vector{Float64} — Five reduction factors applied after failed solves (default: [0.5, 0.5, 0.5, 0.5, 0.1])

params = BackwardEulerSolverParameters(relative_tolerance = 1e-4, max_number_of_steps = 20)
micm = MICM(config_path = path, solver_type = BackwardEulerStandardOrder,
            solver_parameters = params)

SolverState#

type SolverState#

Enum representing the outcome of a solve! call.

Value

Int

Meaning

NotYetCalled

0

solve! has not been called yet

Running

1

Solver is in progress (internal use)

Converged

2

Solution accepted within tolerances

ConvergenceExceededMaxSteps

3

Maximum internal steps reached before convergence

StepSizeTooSmall

4

Step size fell below numerical limit

RepeatedlySingularMatrix

5

Jacobian factorisation failed repeatedly

NaNDetected

6

NaN appeared in the solution

InfDetected

7

Inf appeared in the solution

AcceptingUnconvergedIntegration

8

Solution accepted despite not fully converging

SolverStats#

type SolverStats#

Performance counters returned by solve!.

Fields

  • function_calls::Int — Number of right-hand-side evaluations

  • jacobian_updates::Int — Number of Jacobian computations

  • number_of_steps::Int — Total internal steps taken

  • accepted::Int — Number of accepted steps

  • rejected::Int — Number of rejected steps

  • decompositions::Int — Number of LU decompositions

  • solves::Int — Number of linear system solves

  • final_time::Float64 — Simulated time reached (seconds)

SolverResult#

type SolverResult#

Combined outcome of a solve! call.

Fields

  • state::SolverState — Convergence status

  • stats::SolverStats — Performance counters

result = solve!(micm, state, 60.0)
if result.state == Converged
    println("Took ", result.stats.number_of_steps, " steps")
end

MICM#

type MICM#

Wrapper around the C++ MICM chemical kinetics solver.

Constructor

MICM(; config_path, solver_type=RosenbrockStandardOrder, solver_parameters=nothing)

  • config_path::String — Path to the mechanism configuration file or directory

  • solver_type::SolverType — Solver backend to use (see SolverType)

  • solver_parameters — Optional RosenbrockSolverParameters or BackwardEulerSolverParameters

micm = MICM(config_path = "configs/v0/analytical")
micm_be = MICM(config_path = "configs/v0/analytical",
               solver_type = BackwardEulerStandardOrder)

State#

type State#

Chemical state for one or more grid cells. Created via create_state; do not construct directly.

Holds species concentrations, environmental conditions, and user-defined rate parameters. Internally uses the same vector-ordering layout as the C++ and Python APIs.

Functions#

MICM Functions#

create_state(micm: :MICM; number_of_grid_cells = 1) → State#

Create a new State object for the given solver.

  • number_of_grid_cells — Number of independent atmospheric columns (default: 1)

state = create_state(micm)
state3 = create_state(micm, number_of_grid_cells = 3)
solve!(micm::MICM, state::State, time_step::Real) -> SolverResult

Integrate the chemical system forward by time_step seconds. Species concentrations in state are updated in-place.

result = solve!(micm, state, 60.0)
@assert result.state == Converged
solver_type(micm: :MICM) → SolverType#

Return the solver type this MICM was created with.

set_solver_parameters!(micm::MICM, params)

Update solver tuning parameters. params must match the solver type (RosenbrockSolverParameters for Rosenbrock solvers, BackwardEulerSolverParameters for Backward Euler solvers).

set_solver_parameters!(micm, RosenbrockSolverParameters(relative_tolerance = 1e-8))
get_solver_parameters(micm: :MICM)#

Return the current solver parameters as RosenbrockSolverParameters or BackwardEulerSolverParameters, depending on the solver type.

State Functions#

set_concentrations!(state::State, concentrations::Dict{String})

Set species concentrations. For a single grid cell, values may be scalars. For multiple grid cells, provide a Vector of length number_of_grid_cells.

# Single grid cell
set_concentrations!(state, Dict{String,Any}("A" => 1.0, "B" => 0.0))

# Multiple grid cells
set_concentrations!(state, Dict{String,Any}("A" => [1.0, 2.0, 3.0]))
get_concentrations(state: :State) → Dict{String, Vector{Float64}}#

Return species concentrations for all grid cells. Keys are species names; each value is a vector of length number_of_grid_cells.

set_conditions!(state::State; temperatures, pressures, air_densities)

Set environmental conditions. All keyword arguments accept scalars (single grid cell) or vectors (multiple grid cells). If air_densities is omitted, it is computed from the Ideal Gas Law using the provided temperature and pressure.

# Single cell, auto air_density
set_conditions!(state, temperatures = 298.0, pressures = 101325.0)

# Multiple cells
set_conditions!(state,
    temperatures = [298.0, 310.0, 280.0],
    pressures    = [101325.0, 95000.0, 105000.0])
get_conditions(state: :State) → Dict{String, Vector{Float64}}#

Return environmental conditions for all grid cells. Keys are "temperature" (K), "pressure" (Pa), and "air_density" (mol m⁻³).

set_user_defined_rate_parameters!(state::State, params::Dict{String})

Set user-defined rate parameters (e.g. photolysis rates or emission fluxes). Values may be scalars or vectors matching the number of grid cells.

set_user_defined_rate_parameters!(state,
    Dict{String,Any}("USER.reaction 1" => 0.001, "USER.reaction 2" => 0.002))
get_user_defined_rate_parameters(state: :State) → Dict{String, Vector{Float64}}#

Return user-defined rate parameters for all grid cells.

get_species_ordering(state: :State) → Dict{String, Int}#

Return the mapping of species names to their 0-based indices in the internal concentration array.

get_user_defined_rate_parameters_ordering(state: :State) → Dict{String, Int}#

Return the mapping of user-defined rate parameter names to their 0-based indices.

Multi-Grid-Cell Example#

using Musica

micm = MICM(config_path = "configs/v0/analytical")
state = create_state(micm, number_of_grid_cells = 3)

set_concentrations!(state, Dict{String,Any}(
    "A" => [1.0, 2.0, 3.0],
    "B" => [0.0, 0.0, 0.0],
))

set_conditions!(state,
    temperatures = [298.0, 310.0, 280.0],
    pressures    = [101325.0, 95000.0, 105000.0],
)

set_user_defined_rate_parameters!(state, Dict{String,Any}(
    "USER.reaction 1" => [0.001, 0.002, 0.003],
    "USER.reaction 2" => [0.004, 0.005, 0.006],
))

result = solve!(micm, state, 60.0)
@assert result.state == Converged

concs = get_concentrations(state)
println("Final [A] per cell: ", concs["A"])

Testing#

To run the Julia test suite:

cd julia
julia --project=. test/runtests.jl

Additional Resources#