Skip to content

MUSICA-MICM chemistry solver #1425

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
dwfncar wants to merge 1 commit into
base: develop
Choose a base branch
from
+199 −20
Draft

MUSICA-MICM chemistry solver #1425

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 18 additions & 0 deletions src/core_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1737,6 +1737,15 @@
description="Graupel mixing ratio"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>

<var name="qAB" array_group="passive" units="kg kg^{-1}"
description="Molecular AB mixing ratio"/>

<var name="qA" array_group="passive" units="kg kg^{-1}"
description="Atomic A mixing ratio"/>

<var name="qB" array_group="passive" units="kg kg^{-1}"
description="Atomic B mixing ratio"/>

<var name="ni" array_group="number" units="nb kg^{-1}"
description="Cloud ice number concentration"
packages="bl_mynn_in;mp_thompson_in;mp_thompson_aers_in"/>
Expand Down Expand Up @@ -2103,6 +2112,15 @@
description="Tendency of graupel mass per unit volume divided by d(zeta)/dz"
packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>

<var name="tend_qAB" name_in_code="qAB" array_group="passive" units="kg m^{-3} s^{-1}"
description="Tendency of molecular AB mass per unit volume divided by d(zeta)/dz"/>

<var name="tend_qA" name_in_code="qA" array_group="passive" units="kg m^{-3} s^{-1}"
description="Tendency of atomic A mass per unit volume divided by d(zeta)/dz"/>

<var name="tend_qB" name_in_code="qB" array_group="passive" units="kg m^{-3} s^{-1}"
description="Tendency of atomic B mass per unit volume divided by d(zeta)/dz"/>

<var name="tend_ni" name_in_code="ni" array_group="number" units="nb m^{-3} s^{-1}"
description="Tendency of cloud ice number concentration multiplied by dry air density divided by d(zeta)/dz"
packages="bl_mynn_in;mp_thompson_in;mp_thompson_aers_in"/>
Expand Down
21 changes: 19 additions & 2 deletions src/core_atmosphere/chemistry/mpas_atm_chemistry.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@
!-------------------------------------------------------------------------
module mpas_atm_chemistry

use mpas_kind_types, only : StrKIND, RKIND

implicit none

private
Expand Down Expand Up @@ -86,16 +88,31 @@ end subroutine chemistry_init
!> This routine steps the chemistry packages, updating their state
!> based on the current simulation time and conditions.
!------------------------------------------------------------------------
subroutine chemistry_step()
subroutine chemistry_step(dt)

#ifdef MPAS_USE_MUSICA
use iso_fortran_env, only: real64
use mpas_musica, only: musica_step
#endif
use mpas_log, only : mpas_log_write
use mpas_derived_types, only: MPAS_LOG_CRIT

real (kind=RKIND), intent(in) :: dt

#ifdef MPAS_USE_MUSICA
real(real64) :: time_step
integer :: error_code
character(len=:), allocatable :: error_message

time_step = dt

call mpas_log_write('Stepping chemistry packages...')
! call musica_step()

call musica_step(time_step, error_code, error_message)

if (error_code /= 0) then
call mpas_log_write(error_message, messageType=MPAS_LOG_CRIT)
end if
#endif

end subroutine chemistry_step
Expand Down
169 changes: 152 additions & 17 deletions src/core_atmosphere/chemistry/musica/mpas_musica.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,21 +44,30 @@ module mpas_musica
!> setting up necessary parameters and data structures for the simulation.
!> For now, we will load fixed configurations for MICM and TUV-x.
!> Later, we will gradually replace the fixed configuration elements
!> with runtime updates using the MUSICA API
!> with runtime updates using the MUSICA API.
!------------------------------------------------------------------------
subroutine musica_init(filename_of_micm_configuration, &
number_of_grid_cells, &
error_code, error_message)

use iso_fortran_env, only: real64

use musica_micm, only : RosenbrockStandardOrder
use musica_util, only : error_t, string_t

use mpas_log, only : mpas_log_write

character(len=*), intent(in) :: filename_of_micm_configuration
integer, intent(in) :: number_of_grid_cells
integer, intent(out) :: error_code
character(len=:), allocatable, intent(out) :: error_message
character(len=*), intent(in) :: filename_of_micm_configuration
integer, intent(in) :: number_of_grid_cells
integer, intent(out) :: error_code
character(len=:), allocatable, intent(out) :: error_message

! TEMPORARY: ABBA specific initial concentrations
real(real64), allocatable :: init_ab(:)
real(real64), allocatable :: init_a(:)
real(real64), allocatable :: init_b(:)

character(len=*), parameter :: INITIAL_CONCENTRATION_PROPERTY = "__initial concentration"

type(error_t) :: error
type(string_t) :: micm_version
Expand Down Expand Up @@ -99,25 +108,48 @@ end subroutine musica_init
! subroutine musica_step
!
!> \brief Steps the MUSICA chemistry package
!> \author CheMPAS-A Developers
!> \date 20 August 2025
!> \author David Fillmore
!> \date 17 November 2025
!> \details
!> This subroutine steps the MUSICA chemistry package, updating its state
!> based on the current simulation time and conditions.
!> It first calls the TUV-x package to compute photolysis rates,
!> then calls the MICM package to solve the ODEs for chemical species
!> concentrations.
!> This subroutine steps the MICM ODE solver.
!------------------------------------------------------------------------
subroutine musica_step()
subroutine musica_step(time_step, error_code, error_message)

use iso_fortran_env, only: real64

use musica_micm, only : solver_stats_t
use musica_util, only : error_t, string_t

use mpas_log, only : mpas_log_write

real(real64), intent(in) :: time_step
integer, intent(out) :: error_code
character(len=:), allocatable, intent(out) :: error_message

type(string_t) :: solver_state
type(solver_stats_t) :: solver_stats
type(error_t) :: error

if (.not. musica_is_initialized) return

call mpas_log_write('Stepping MUSICA chemistry package...')
call mpas_log_write('[MUSICA] Stepping MICM solver...')

call micm%solve(time_step, state, solver_state, solver_stats, error)
if (has_error_occurred(error, error_message, error_code)) return

if (solver_state%get_char_array() /= 'Converged') then
call mpas_log_write('[MUSICA Warning] MICM solver failure')
end if

call log_solver_stats(solver_stats)

! Here we would typically call the TUV-x and MICM packages to perform
! the necessary computations, but for now we will just log the step.
! TEMPORARY: ABBA specific logging
call mpas_log_write('[MUSICA] MICM AB concentrations')
call log_concentrations(state, state%number_of_grid_cells, 'AB', error)
call mpas_log_write('[MUSICA] MICM A concentrations')
call log_concentrations(state, state%number_of_grid_cells, 'A', error)
call mpas_log_write('[MUSICA] MICM B concentrations')
call log_concentrations(state, state%number_of_grid_cells, 'B', error)

end subroutine musica_step

Expand All @@ -143,6 +175,110 @@ subroutine musica_finalize()

end subroutine musica_finalize

subroutine assign_rate_parameters(state, number_of_grid_cells)
! Provide a default value for any rate parameters (none for ABBA by default).
use iso_fortran_env, only : real64

type(state_t), intent(inout) :: state
integer, intent(in) :: number_of_grid_cells

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, number_of_grid_cells
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

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
use musica_util, only : error_t, string_t

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

subroutine log_concentrations(state, number_of_grid_cells, species_name, err)
use iso_fortran_env, only : real64
use mpas_log, only : mpas_log_write
use mpas_kind_types, only : RKIND
use musica_util, only : error_t, string_t

type(state_t), intent(in) :: state
integer, intent(in) :: number_of_grid_cells
character(len=*), intent(in) :: species_name
type(error_t), intent(inout) :: err

integer :: species_index, cell_stride, var_stride, cell, idx
real (kind=RKIND) :: concentration

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, number_of_grid_cells
idx = 1 + (cell - 1) * cell_stride + (species_index - 1) * var_stride
concentration = state%concentrations(idx)
call mpas_log_write(' MICM concentration: $r', realArgs=[concentration])
end do
end subroutine log_concentrations

subroutine log_solver_stats(solver_stats)
use mpas_log, only : mpas_log_write
use mpas_kind_types, only : RKIND
use musica_micm, only : solver_stats_t

type(solver_stats_t), intent(in) :: solver_stats

integer :: function_calls, jacobian_updates, number_of_steps
integer :: accepted, rejected, decompositions, solves

real (kind=RKIND) :: final_time

function_calls = solver_stats%function_calls()
jacobian_updates = solver_stats%jacobian_updates()
number_of_steps = solver_stats%number_of_steps()
accepted = solver_stats%accepted()
rejected = solver_stats%rejected()
decompositions = solver_stats%decompositions()
solves = solver_stats%solves()
final_time = solver_stats%final_time()

call mpas_log_write('[MUSICA] MICM Solver statistics ...')
call mpas_log_write(' MICM function calls: $i', intArgs=[function_calls])
call mpas_log_write(' MICM jacobian updates: $i', intArgs=[jacobian_updates])
call mpas_log_write(' MICM number of steps: $i', intArgs=[number_of_steps])
call mpas_log_write(' MICM accepted: $i', intArgs=[accepted])
call mpas_log_write(' MICM rejected: $i', intArgs=[rejected])
call mpas_log_write(' MICM LU decompositions: $i', intArgs=[decompositions])
call mpas_log_write(' MICM linear solves: $i', intArgs=[solves])
call mpas_log_write(' MICM final time (s): $r', realArgs=[final_time])
end subroutine log_solver_stats

!-------------------------------------------------------------------------
! function has_error_occurred
!
Expand Down Expand Up @@ -175,7 +311,6 @@ logical function has_error_occurred(error, error_message, error_code)
error_message = '[MUSICA Error]: ' // error%category( ) // '[' // &
trim( adjustl( error_code_str ) ) // ']: ' // error%message( )
has_error_occurred = .true.

end function has_error_occurred

end module mpas_musica
2 changes: 1 addition & 1 deletion src/core_atmosphere/mpas_atm_core.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1038,7 +1038,7 @@ subroutine atm_do_timestep(domain, dt, itimestep)
#endif

#ifdef DO_CHEMISTRY
call chemistry_step()
call chemistry_step(dt)
#endif

call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group)
Expand Down
9 changes: 9 additions & 0 deletions src/core_init_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1163,6 +1163,15 @@
<var name="qr" array_group="moist" units="kg kg^{-1}"
description="Rain water mixing ratio"/>

<var name="qAB" array_group="passive" units="kg kg^{-1}"
description="Molecular AB mixing ratio"/>

<var name="qA" array_group="passive" units="kg kg^{-1}"
description="Atomic A mixing ratio"/>

<var name="qB" array_group="passive" units="kg kg^{-1}"
description="Atomic B mixing ratio"/>

<var name="nifa" array_group="number" units="nb kg^{-1}"
description="Gocart ice-friendly aerosol number concentration"
packages="mp_thompson_aers_in"/>
Expand Down