Add Julia/C++ bindings for MUSICA

.github/workflows/julia.yml

@@ -33,12 +33,28 @@ jobs:
         working-directory: julia
         run: julia --project=. -e 'using Pkg; Pkg.instantiate()'
 
+      - name: Get CxxWrap prefix path
+        id: cxxwrap
+        working-directory: julia
+        run: echo "prefix=$(julia --project=. -e 'using CxxWrap; print(CxxWrap.prefix_path())')" >> $GITHUB_OUTPUT
+
       - name: Configure MUSICA with Julia support
-        run: cmake -S . -B build -D CMAKE_BUILD_TYPE=Release -D MUSICA_ENABLE_JULIA=ON -D MUSICA_ENABLE_MICM=ON -D MUSICA_ENABLE_TUVX=OFF -D MUSICA_ENABLE_CARMA=OFF -D MUSICA_ENABLE_TESTS=OFF -D CMAKE_POLICY_VERSION_MINIMUM=3.5
+        run: cmake -S . -B build -D CMAKE_BUILD_TYPE=Release -D MUSICA_ENABLE_JULIA=ON -D MUSICA_ENABLE_MICM=ON -D MUSICA_ENABLE_TUVX=OFF -D MUSICA_ENABLE_CARMA=OFF -D MUSICA_ENABLE_TESTS=OFF -D CMAKE_POSITION_INDEPENDENT_CODE=ON -D CMAKE_PREFIX_PATH=${{ steps.cxxwrap.outputs.prefix }}
 
       - name: Build MUSICA with Julia support
         run: cmake --build build --verbose
 
+      - name: Set MUSICA_JULIA_LIB path
+        working-directory: julia
+        run: |
+            if [ -f deps/lib/libmusica_julia.so ]; then
+              echo "MUSICA_JULIA_LIB=$(pwd)/deps/lib/libmusica_julia.so" >> $GITHUB_ENV
+            elif [ -f deps/lib/libmusica_julia.dylib ]; then
+              echo "MUSICA_JULIA_LIB=$(pwd)/deps/lib/libmusica_julia.dylib" >> $GITHUB_ENV
+            else
+              echo "ERROR: libmusica_julia not found" && ls -R deps/lib/ && exit 1
+            fi
+
       - name: Precompile Julia package with built library
         working-directory: julia
         run: julia --project=. -e 'using Pkg; Pkg.precompile()'

docker/Dockerfile.julia

@@ -24,9 +24,10 @@ COPY --from=julia-base ${JULIA_PATH} ${JULIA_PATH}
 # Copy the musica code
 COPY . musica
 
-# Install Julia dependencies
+# Install Julia dependencies and get CxxWrap prefix path
 RUN cd musica/julia \
-    && julia --project=. -e 'using Pkg; Pkg.instantiate()'
+    && julia --project=. -e 'using Pkg; Pkg.instantiate()' \
+    && julia --project=. -e 'using CxxWrap; open("/tmp/cxxwrap_prefix", "w") do f; print(f, CxxWrap.prefix_path()); end'
 
 # Build
 RUN cd musica \
@@ -39,10 +40,13 @@ RUN cd musica \
              -D MUSICA_ENABLE_TUVX=OFF \
              -D MUSICA_ENABLE_CARMA=OFF \
              -D MUSICA_ENABLE_TESTS=OFF \
+             -D CMAKE_POSITION_INDEPENDENT_CODE=ON \
+             -D CMAKE_PREFIX_PATH=$(cat /tmp/cxxwrap_prefix) \
     && cmake --build build
 
-# Precompile and test Julia package
+# Set library path and precompile Julia package
+ENV MUSICA_JULIA_LIB=/musica/julia/deps/lib/libmusica_julia.so
 RUN cd musica/julia \
-    && julia --project=. -e 'using Pkg; Pkg.precompile()' 
+    && julia --project=. -e 'using Pkg; Pkg.precompile()'
 
 WORKDIR /musica

julia/CMakeLists.txt

@@ -7,11 +7,12 @@ add_library(musica_julia SHARED
 
 target_compile_features(musica_julia PUBLIC cxx_std_20)
 
-set_target_properties(musica PROPERTIES POSITION_INDEPENDENT_CODE ON)
+set_target_properties(musica_julia PROPERTIES POSITION_INDEPENDENT_CODE ON)
 
 target_link_libraries(musica_julia
     musica::musica
     JlCxx::cxxwrap_julia
+    JlCxx::cxxwrap_julia_stl
 )
 
 # Set output properties

julia/Project.toml

@@ -5,15 +5,14 @@ version = "0.14.4"
 
 [deps]
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
-Musica_jll = "f5e5459d-1a87-5129-97d9-bab22ca84fbb"
 
 [extras]
+Musica_jll = "f5e5459d-1a87-5129-97d9-bab22ca84fbb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test"]
 
 [compat]
 CxxWrap = "^0.16"
-Musica_jll = "0.14.4"
 julia = "1.10 - 1.11"

julia/src/Musica.jl

@@ -1,19 +1,79 @@
+# Copyright (C) 2023-2026 University Corporation for Atmospheric Research
+# SPDX-License-Identifier: Apache-2.0
+
 module Musica
 
 using CxxWrap
-using Musica_jll
 
-function __init__()
-    # Allow overriding the JLL library for raw CMake development builds
-    lib = get(ENV, "MUSICA_JULIA_LIB", nothing)
-    if lib !== nothing && isfile(lib)
-        @wrapmodule(() -> lib)
+# Determine library path at compile time (must be top-level for @wrapmodule)
+const _lib_path = let
+    env_lib = get(ENV, "MUSICA_JULIA_LIB", nothing)
+    if env_lib !== nothing && isfile(env_lib)
+        env_lib
     else
-        @wrapmodule(() -> Musica_jll.libmusica_julia)
+        try
+            @eval using Musica_jll
+            Musica_jll.libmusica_julia
+        catch
+            error("MUSICA_JULIA_LIB environment variable must be set to the path of libmusica_julia.so " *
+                  "(or install Musica_jll)")
+        end
     end
+end
+
+# @wrapmodule must be at module top level so types are defined during precompilation
+@wrapmodule(() -> _lib_path)
+
+function __init__()
     @initcxx
 end
 
+# Type aliases for CxxWrap pointer types (used in MICM and State structs)
+const MICMPtr = CxxWrap.CxxWrapCore.CxxPtr{CppMICM}
+const StatePtr = CxxWrap.CxxWrapCore.CxxPtr{CppState}
+
+# Include MICM submodule files (order matters for dependencies)
+include("micm/constants.jl")
+include("micm/solver.jl")
+include("micm/solver_result.jl")
+include("micm/conditions.jl")
+include("micm/solver_parameters.jl")
+include("micm/utils.jl")
+include("micm/state.jl")
+include("micm/micm.jl")
+
+# Version
 export get_version
 
-end
+# Constants
+export AVOGADRO, BOLTZMANN, GAS_CONSTANT
+
+# Solver types
+export SolverType
+export Rosenbrock, RosenbrockStandardOrder, BackwardEuler, BackwardEulerStandardOrder, CudaRosenbrock
+
+# Solver results
+export SolverState, SolverStats, SolverResult
+export NotYetCalled, Running, Converged, ConvergenceExceededMaxSteps
+export StepSizeTooSmall, RepeatedlySingularMatrix, NaNDetected, InfDetected
+export AcceptingUnconvergedIntegration
+
+# Conditions
+export Conditions
+
+# Solver parameters
+export RosenbrockSolverParameters, BackwardEulerSolverParameters
+
+# MICM and State types
+export MICM, State
+
+# Functions
+export create_state, solve!
+export set_concentrations!, get_concentrations
+export set_conditions!, get_conditions
+export set_user_defined_rate_parameters!, get_user_defined_rate_parameters
+export get_species_ordering, get_user_defined_rate_parameters_ordering
+export set_solver_parameters!, get_solver_parameters
+export solver_type
+
+end # module Musica

julia/src/micm/conditions.jl

@@ -0,0 +1,32 @@
+# Copyright (C) 2023-2026 University Corporation for Atmospheric Research
+# SPDX-License-Identifier: Apache-2.0
+
+"""
+    Conditions
+
+Environmental conditions for a grid cell.
+
+If `air_density` is not provided and both `temperature` and `pressure` are given,
+air density is calculated from the Ideal Gas Law.
+
+# Fields
+- `temperature::Float64`: Temperature in Kelvin
+- `pressure::Float64`: Pressure in Pascals
+- `air_density::Float64`: Air density in mol m⁻³
+"""
+mutable struct Conditions
+    temperature::Float64
+    pressure::Float64
+    air_density::Float64
+end
+
+function Conditions(; temperature::Real=0.0, pressure::Real=0.0, air_density::Union{Real, Nothing}=nothing)
+    if air_density === nothing
+        if temperature > 0.0 && pressure > 0.0
+            air_density = pressure / (GAS_CONSTANT * temperature)
+        else
+            air_density = 0.0
+        end
+    end
+    return Conditions(Float64(temperature), Float64(pressure), Float64(air_density))
+end

julia/src/micm/constants.jl

@@ -0,0 +1,6 @@
+# Copyright (C) 2023-2026 University Corporation for Atmospheric Research
+# SPDX-License-Identifier: Apache-2.0
+
+const AVOGADRO = 6.02214076e23       # mol^-1
+const BOLTZMANN = 1.380649e-23       # J K^-1
+const GAS_CONSTANT = AVOGADRO * BOLTZMANN  # J K^-1 mol^-1

julia/src/micm/micm.jl

@@ -0,0 +1,101 @@
+# Copyright (C) 2023-2026 University Corporation for Atmospheric Research
+# SPDX-License-Identifier: Apache-2.0
+
+"""
+    MICM
+
+Wrapper around the C++ MICM chemical kinetics solver.
+
+# Constructor
+
+    MICM(; config_path, solver_type=RosenbrockStandardOrder, solver_parameters=nothing)
+
+- `config_path::String`: Path to configuration file or directory
+- `solver_type::SolverType`: Type of solver to use
+- `solver_parameters`: Optional `RosenbrockSolverParameters` or `BackwardEulerSolverParameters`
+"""
+mutable struct MICM
+    _ptr::MICMPtr
+    _solver_type::SolverType
+    _vector_size::Int
+
+    function MICM(;
+        config_path::String,
+        solver_type::SolverType=RosenbrockStandardOrder,
+        solver_parameters::Union{RosenbrockSolverParameters, BackwardEulerSolverParameters, Nothing}=nothing,
+    )
+        vs = Int(cpp_get_vector_size(Int(solver_type)))
+        vs > 0 || error("Invalid vector size: $vs")
+
+        ptr = cpp_create_solver(config_path, Int(solver_type))
+
+        obj = new(ptr, solver_type, vs)
+        finalizer(obj) do m
+            cpp_delete_solver(m._ptr)
+        end
+
+        if solver_parameters !== nothing
+            set_solver_parameters!(obj, solver_parameters)
+        end
+        return obj
+    end
+end
+
+"""
+    solver_type(micm::MICM) -> SolverType
+
+Get the type of solver used.
+"""
+solver_type(micm::MICM) = micm._solver_type
+
+"""
+    create_state(micm::MICM; number_of_grid_cells=1) -> State
+
+Create a new state object for this solver.
+"""
+function create_state(micm::MICM; number_of_grid_cells::Int=1)
+    return State(micm._ptr, number_of_grid_cells, micm._vector_size, micm)
+end
+
+"""
+    solve!(micm::MICM, state::State, time_step::Real) -> SolverResult
+
+Solve the chemical system for the given state and time step (in seconds).
+"""
+function solve!(micm::MICM, state::State, time_step::Real)
+    cpp_result = cpp_micm_solve(micm._ptr, state._ptr, Float64(time_step))
+    return SolverResult(cpp_result)
+end
+
+"""
+    set_solver_parameters!(micm::MICM, params::RosenbrockSolverParameters)
+    set_solver_parameters!(micm::MICM, params::BackwardEulerSolverParameters)
+
+Set solver-specific parameters. Parameter type must match the solver type.
+"""
+function set_solver_parameters!(micm::MICM, params::RosenbrockSolverParameters)
+    cpp_params = to_cpp_rosenbrock(params)
+    cpp_set_rosenbrock_params(micm._ptr, cpp_params)
+end
+
+function set_solver_parameters!(micm::MICM, params::BackwardEulerSolverParameters)
+    cpp_params = to_cpp_backward_euler(params)
+    cpp_set_backward_euler_params(micm._ptr, cpp_params)
+end
+
+"""
+    get_solver_parameters(micm::MICM)
+
+Get the current solver parameters. Returns `RosenbrockSolverParameters` or
+`BackwardEulerSolverParameters` depending on the solver type.
+"""
+function get_solver_parameters(micm::MICM)
+    st = micm._solver_type
+    if st == Rosenbrock || st == RosenbrockStandardOrder || st == CudaRosenbrock
+        return from_cpp_rosenbrock(cpp_get_rosenbrock_params(micm._ptr))
+    elseif st == BackwardEuler || st == BackwardEulerStandardOrder
+        return from_cpp_backward_euler(cpp_get_backward_euler_params(micm._ptr))
+    else
+        error("Unknown solver type: $st")
+    end
+end

julia/src/micm/solver.jl

@@ -0,0 +1,17 @@
+# Copyright (C) 2023-2026 University Corporation for Atmospheric Research
+# SPDX-License-Identifier: Apache-2.0
+
+"""
+    SolverType
+
+Enum representing the type of MICM solver to use.
+
+Values match the C++ `musica::MICMSolver` enum.
+"""
+@enum SolverType begin
+    Rosenbrock = 1
+    RosenbrockStandardOrder = 2
+    BackwardEuler = 3
+    BackwardEulerStandardOrder = 4
+    CudaRosenbrock = 5
+end