ABBA

ABBA Box/Column Model (MUSICA MICM)

Simple MUSICA MICM Testing

MUSICA MICM Fortran API
Simple verifiable solutions
Number of levels can be 1 (box model) or more (column model)
User defined MICM reactions so concentration units can be arbitrary
Run stand alone and running within CheMPAS

System

$$ AB \leftrightarrow A + B $$

Expected equilibrium

With rate constants $k_f = 2.0$ for $AB \rightarrow A + B$ and $k_r = 1.0$ for $A + B \rightarrow AB$, starting from $[AB] = 1$, $[A] = 0$, $[B] = 0$ in a closed box:

$x$ is the amount of AB converted:

$$ \begin{aligned} \\ [AB] &= 1 - x \\ [A] &= x \\ [B] &= x \end{aligned} $$

Equilibrium $k_f [AB] = k_r [A][B]$

$$ 2 (1 - x) = x^2 \quad\Rightarrow\quad x^2 + 2x - 2 = 0 $$

$$ x = -1 + \sqrt{3} \approx 0.732 $$

So at equilibrium:

$[A] ~ 0.732$
$[B] ~ 0.732$
$[AB] ~ 0.268$

These values match the defaults if USER_DEFINED reactions follow simple mass-action kinetics.

Single-grid-cell box model using the MUSICA MICM Fortran API and the ABBA user-defined mechanism (abba.yaml). Defaults: 1 s time step, 4 s total integration, initial AB=1, A=0, B=0. Prints a formatted table of time, concentrations, and solver step counts per call plus solver statistics.

Files

abba_box.F90 – driver program (shown in sections below).
abba.yaml – three-species mechanism with two USER_DEFINED reactions (AB → A+B and reverse).
Makefile – builds abba_box.x via pkg-config musica-fortran.

Build

make

Ensure MUSICA is loaded or PKG_CONFIG_PATH points to musica-fortran.pc.

Run

./abba_box.x [dt_seconds] [total_time_seconds]

Examples:

./abba_box.x # dt=1.0, total_time=4.0
./abba_box.x 0.5 8.0 # dt=0.5, total_time=8.0

Program Walkthrough (complete code, sectioned)

Mechanism (`abba.yaml`)

version: "1.0.0"
name: ABBA
species:
  - name: A
    __absolute tolerance: 1e-4
    __long name: atom_A
    __do advect: true
  - name: B
    __absolute tolerance: 1e-4
    __long name: atom_B
    __do advect: true
  - name: AB
    __absolute tolerance: 1e-4
    __long name: molecule_AB
    __do advect: true

Mechanism (phases and reactions)

phases:
  - name: gas
    species:
      - A
      - B
      - AB

reactions:
  - type: USER_DEFINED
    gas phase: gas
    reactants:
      - species name: AB
        coefficient: 1
    products:
      - species name: A
        coefficient: 1
      - species name: B
        coefficient: 1
    scaling factor: 2.0
    name: forward_AB_to_A_B

  - type: USER_DEFINED
    gas phase: gas
    reactants:
      - species name: A
        coefficient: 1
      - species name: B
        coefficient: 1
    products:
      - species name: AB
        coefficient: 1
    scaling factor: 1.0
    name: reverse_A_B_to_AB

Program Walkthrough (complete code, sectioned)

Program header and constants

program abba_box
  ! Single-grid-cell ABBA box driver using the MUSICA MICM Fortran interface.
  use iso_fortran_env, only : real64, int64
  use musica_util,       only : error_t, string_t
  use musica_state,      only : state_t
  use musica_micm,       only : micm_t, solver_stats_t, RosenbrockStandardOrder

  implicit none

  integer,          parameter :: NUM_LEVELS = 1
  real(real64),     parameter :: GAS_CONSTANT = 8.314_real64
  real(real64),     parameter :: DEFAULT_TEMPERATURE = 300.0_real64
  real(real64),     parameter :: DEFAULT_PRESSURE = 101325._real64
  real(real64),     parameter :: DEFAULT_TOTAL_TIME = 4.0_real64
  character(len=*), parameter :: CONFIG_FILE = "abba.yaml"

State, arrays, and runtime options

  type(error_t)           :: error
  type(micm_t), pointer   :: micm => null()
  type(state_t), pointer  :: state => null()
  type(string_t)          :: solver_state
  type(solver_stats_t)    :: solver_stats
  real(real64), allocatable :: temperature(:), pressure(:), air_density(:)
  real(real64), allocatable :: init_ab(:), init_a(:), init_b(:)
  real(real64) :: time_step, elapsed, total_time
  integer :: num_steps, step, ios
  integer(int64) :: step_solver_steps
  integer :: ab_index, a_index, b_index, cell_stride, var_stride
  character(len=256) :: arg

Allocation and initial conditions

  allocate(temperature(NUM_LEVELS), pressure(NUM_LEVELS), &
           air_density(NUM_LEVELS), init_ab(NUM_LEVELS), init_a(NUM_LEVELS), &
           init_b(NUM_LEVELS))

  temperature = DEFAULT_TEMPERATURE
  pressure = DEFAULT_PRESSURE
  air_density = pressure / (GAS_CONSTANT * temperature)
  init_ab = 1.0_real64
  init_a = 0.0_real64
  init_b = 0.0_real64

CLI overrides and run setup

  time_step = 1.0_real64
  total_time = DEFAULT_TOTAL_TIME
  num_steps = 0

  if (command_argument_count() >= 1) then
    ! Optional override of time step (seconds).
    call get_command_argument(1, arg)
    read(arg, *, iostat=ios) time_step
    if (ios /= 0 .or. time_step <= 0._real64) time_step = 1.0_real64
  end if
  if (command_argument_count() >= 2) then
    ! Optional override of total integration time (seconds).
    call get_command_argument(2, arg)
    read(arg, *, iostat=ios) total_time
    if (ios /= 0 .or. total_time <= 0._real64) total_time = DEFAULT_TOTAL_TIME
  end if
  num_steps = max(1, int(total_time / time_step))

  write(*,'(a,f8.3,a,i0,a,f8.3)') "Starting ABBA box with dt=", time_step, &
       " s, steps=", num_steps, ", total T=", total_time

MICM setup and initial state

  micm => micm_t(CONFIG_FILE, RosenbrockStandardOrder, error)
  call ensure_success("creating MICM solver", error)

  state => micm%get_state(NUM_LEVELS, error)
  call ensure_success("allocating MICM state", error)

  call load_conditions(state, temperature, pressure, air_density)
  call set_species_profile(state, "AB", init_ab, error)
  call ensure_success("setting AB initial conditions", error)
  call set_species_profile(state, "A", init_a, error)
  call ensure_success("setting A initial conditions", error)
  call set_species_profile(state, "B", init_b, error)
  call ensure_success("setting B initial conditions", error)
  call assign_rate_parameters(state)

Species indexing and initial table row

  cell_stride = state%species_strides%grid_cell
  var_stride  = state%species_strides%variable
  ab_index = state%species_ordering%index("AB", error)
  call ensure_success("retrieving AB index", error)
  a_index = state%species_ordering%index("A", error)
  call ensure_success("retrieving A index", error)
  b_index = state%species_ordering%index("B", error)
  call ensure_success("retrieving B index", error)

  call print_table_header()
  elapsed = 0.0_real64
  call print_table_row(state, elapsed, cell_stride, var_stride, ab_index, &
       a_index, b_index, 0_int64)

Time-stepping loop

  do step = 1, num_steps
    call micm%solve(time_step, state, solver_state, solver_stats, error)
    call ensure_success("MICM solve step", error)
    step_solver_steps = solver_stats%number_of_steps()
    elapsed = elapsed + time_step
    call print_table_row(state, elapsed, cell_stride, var_stride, ab_index, &
         a_index, b_index, step_solver_steps)
  end do
  write(*,'(a)') "Chemistry solver status: " // &
       trim(solver_state%get_char_array())
  write(*,'(a,i0)') "Total internal steps accepted: ", solver_stats%accepted()
  call print_solver_stats(solver_stats)

  if (associated(state)) deallocate(state)
  if (associated(micm)) deallocate(micm)

Error handling helper

contains

  subroutine ensure_success(context, err)
    ! Abort with a readable message if an error is present.
    character(len=*), intent(in) :: context
    type(error_t),    intent(in) :: err
    character(len=:), allocatable :: category, message
    if (err%is_success()) return
    category = err%category()
    message = err%message()
    write(*,'(a)') "[error] " // trim(context) // ": " // trim(category) &
         // " - " // trim(message)
    stop 1
  end subroutine ensure_success

Condition loading and species setters

  subroutine load_conditions(state, t, p, rho)
    ! Copy simple meteorology into the MICM state structure.
    type(state_t), intent(inout) :: state
    real(real64), intent(in)     :: t(:), p(:), rho(:)
    integer :: cell
    do cell = 1, size(t)
      state%conditions(cell)%temperature = t(cell)
      state%conditions(cell)%pressure = p(cell)
      state%conditions(cell)%air_density = rho(cell)
    end do
  end subroutine load_conditions

  subroutine set_species_profile(state, species_name, values, err)
    ! Write a 1-D profile into the contiguous concentrations array.
    use iso_fortran_env, only : real64
    type(state_t), intent(inout) :: state
    character(len=*), intent(in) :: species_name
    real(real64), intent(in)     :: values(:)
    type(error_t), intent(inout) :: err
    integer :: species_index, cell_stride, var_stride, cell, idx

    species_index = state%species_ordering%index(species_name, err)
    if (.not. err%is_success()) return

    cell_stride = state%species_strides%grid_cell
    var_stride  = state%species_strides%variable

    do cell = 1, size(values)
      idx = 1 + (cell - 1) * cell_stride + (species_index - 1) * var_stride
      state%concentrations(idx) = values(cell)
    end do
  end subroutine set_species_profile

Rate parameters (placeholder)

  subroutine assign_rate_parameters(state)
    ! Provide a default value for any rate parameters (none for ABBA by default).
    type(state_t), intent(inout) :: state
    integer :: cell_stride, var_stride, rp_index, cell, idx
    if (state%number_of_rate_parameters == 0) return

    cell_stride = state%rate_parameters_strides%grid_cell
    var_stride  = state%rate_parameters_strides%variable

    do rp_index = 1, state%rate_parameters_ordering%size()
      do cell = 1, NUM_LEVELS
        idx = 1 + (cell - 1) * cell_stride + (rp_index - 1) * var_stride
        state%rate_parameters(idx) = 1.0_real64
      end do
    end do
  end subroutine assign_rate_parameters

Tabular output helpers

  subroutine print_table_header()
    write(*,'(/,1x,a12,3x,a12,3x,a12,3x,a12,3x,a12)') "time (s)", "AB", "A", &
         "B", "steps"
    write(*,'(a)')    " ----------------------------------------------------------------"
  end subroutine print_table_header

  subroutine print_table_row(state, time_s, cell_stride, var_stride, &
       ab_idx, a_idx, b_idx, n_steps)
    ! Emit a single formatted row of time and concentrations.
    type(state_t), intent(in) :: state
    real(real64), intent(in)  :: time_s
    integer,      intent(in)  :: cell_stride, var_stride
    integer,      intent(in)  :: ab_idx, a_idx, b_idx
    integer(int64), intent(in) :: n_steps
    integer :: cell, idx_ab, idx_a, idx_b
    cell = 1
    idx_ab = 1 + (cell - 1) * cell_stride + (ab_idx - 1) * var_stride
    idx_a  = 1 + (cell - 1) * cell_stride + (a_idx  - 1) * var_stride
    idx_b  = 1 + (cell - 1) * cell_stride + (b_idx  - 1) * var_stride
    write(*,'(1x,f12.2,3x,f12.4,3x,f12.4,3x,f12.4,3x,i12)') time_s, &
         state%concentrations(idx_ab), state%concentrations(idx_a), &
         state%concentrations(idx_b), n_steps
  end subroutine print_table_row

Solver statistics dump

  subroutine print_solver_stats(stats)
    type(solver_stats_t), intent(in) :: stats
  write(*,'(/,a)') "Solver statistics:"
  write(*,'(a,i0)') "  function calls: ", stats%function_calls()
  write(*,'(a,i0)') "  jacobian updates: ", stats%jacobian_updates()
  write(*,'(a,i0)') "  number of steps: ", stats%number_of_steps()
  write(*,'(a,i0)') "  accepted: ", stats%accepted()
  write(*,'(a,i0)') "  rejected: ", stats%rejected()
  write(*,'(a,i0)') "  LU decompositions: ", stats%decompositions()
  write(*,'(a,i0)') "  linear solves: ", stats%solves()
  write(*,'(a,f10.3)') "  final time (s): ", stats%final_time()
  end subroutine print_solver_stats

end program abba_box

Notes

The per-step steps column shows the internal solver step count for that micm%solve call (not cumulative).
For cumulative stats, use the printed solver statistics after the loop.

Example output (dt=1.0 s, total_time=4.0 s)

Starting ABBA box with dt=   1.000 s, steps=4, total T=   4.000

     time (s)             AB              A              B          steps
 ----------------------------------------------------------------
         0.00         1.0000         0.0000         0.0000              0
         1.00         0.2965         0.7035         0.7035             14
         2.00         0.2688         0.7312         0.7312             10
         3.00         0.2680         0.7320         0.7320              9
         4.00         0.2679         0.7321         0.7321              9
Chemistry solver status: Converged
Total internal steps accepted: 9

Solver statistics:
  function calls: 18
  jacobian updates: 9
  number of steps: 9
  accepted: 9
  rejected: 0
  LU decompositions: 9
  linear solves: 27
  final time (s):      1.000

ABBA

ABBA Box/Column Model (MUSICA MICM)

Simple MUSICA MICM Testing

System

Expected equilibrium

Files

Build

Run

Program Walkthrough (complete code, sectioned)

Mechanism (abba.yaml)

Mechanism (phases and reactions)

Program Walkthrough (complete code, sectioned)

Program header and constants

State, arrays, and runtime options

Allocation and initial conditions

CLI overrides and run setup

MICM setup and initial state

Species indexing and initial table row

Time-stepping loop

Error handling helper

Condition loading and species setters

Rate parameters (placeholder)

Tabular output helpers

Solver statistics dump

Notes

Example output (dt=1.0 s, total_time=4.0 s)

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally

Mechanism (`abba.yaml`)