Skip to content

Commit e014f6f

Browse files
Copilotmattldawsonboulderdaze
authored
Add Tutorial Chapter 3: Multiple grid cells MICM box model (#509)
* Initial plan * Create Chapter 3 tutorial and Fortran example for multiple grid cells MICM box model Co-authored-by: mattldawson <[email protected]> * Fix multi-grid cell array indexing and improve documentation accuracy Co-authored-by: mattldawson <[email protected]> * address the review comments * address the review comments * fix a bug * add missing cmake scripts for a test * fix bugs in test * remove analysis * address the review comment --------- Co-authored-by: copilot-swe-agent[bot] <[email protected]> Co-authored-by: mattldawson <[email protected]> Co-authored-by: Jiwon Gim <[email protected]>
1 parent c7a0aee commit e014f6f

File tree

5 files changed

+304
-2
lines changed

5 files changed

+304
-2
lines changed

‎docs/source/tutorials/chapter3.rst‎

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
.. _chapter3:
2+
3+
Chapter 3
4+
=========
5+
6+
MICM Box Model with Multiple Grid Cells in Fortran
7+
--------------------------------------------------
8+
9+
In this third MUSICA Fortran example,
10+
we will extend the MICM box model from Chapter 2 to handle multiple grid cells
11+
in a single call to the MICM solver.
12+
This demonstrates how MUSICA can efficiently solve chemical mechanisms
13+
for multiple independent air masses simultaneously.
14+
15+
Each grid cell represents an independent set of well-mixed air masses,
16+
allowing us to model different atmospheric conditions (temperature, pressure, concentrations)
17+
and observe how the same chemical mechanism evolves differently under various conditions.
18+
19+
We will use the same analytical configuration files from Chapter 2
20+
(``config.json``, ``species.json``, and ``reactions.json``)
21+
saved in the subdirectory ``configs/v0/analytical``.
22+
23+
For reference, these files define:
24+
25+
- A system of six chemical species: A, B, C, D, E, and F
26+
- Four reactions including two Arrhenius-type reactions and two user-defined rate reactions
27+
- The same mechanism will be solved simultaneously across multiple grid cells
28+
29+
To create a multiple grid cells box model, save the following Fortran code to a file named ``micm_multiple_grid_cells.F90``:
30+
31+
.. literalinclude:: ../../../fortran/test/fetch_content_integration/test_micm_multiple_grid_cells.F90
32+
:language: f90
33+
34+
Key differences from the single grid cell example in Chapter 2:
35+
36+
**Grid Cell Configuration:**
37+
38+
- ``num_grid_cells`` is set to 3 instead of 1
39+
- Each grid cell can have different environmental conditions (temperature, pressure)
40+
- Each grid cell can start with different initial concentrations
41+
42+
**Setting Conditions and Initial Values:**
43+
44+
The ``state`` object now manages data for multiple grid cells using strides:
45+
46+
- ``state%species_strides%variable`` provides the stride for accessing species data across cells
47+
- ``state%rate_parameters_strides%variable`` provides the stride for rate parameter data
48+
- Concentrations are stored in a flat array with stride-based indexing
49+
50+
.. code-block:: fortran
51+
52+
! Different temperatures for each cell
53+
do cell_id = 1, num_grid_cells
54+
state%conditions(cell_id)%temperature = 273.0 + (cell_id - 1) * 10.0
55+
state%conditions(cell_id)%pressure = 1.0e5
56+
state%conditions(cell_id)%air_density = state%conditions(cell_id)%pressure / &
57+
(GAS_CONSTANT * state%conditions(cell_id)%temperature)
58+
end do
59+
60+
**Accessing Multi-Cell Data:**
61+
62+
To access concentration data for a specific species in a specific grid cell:
63+
64+
.. code-block:: fortran
65+
66+
cell_stride = state%species_strides%grid_cell
67+
var_stride = state%species_strides%variable
68+
cell_concentration = state%concentrations(1 + (cell_id - 1) * cell_stride + (species_index - 1) * var_stride)
69+
70+
Fortran's 1-based indexing requires subtracting 1 from the 0-based ``cell_id``
71+
and ``species_index`` to adjust offsets.
72+
73+
**Simultaneous Solving:**
74+
75+
A call to ``micm%solve`` with a state that represents multiple grid cells is only
76+
required once per time step; all grid cells are solved simultaneously in a multi-grid cell micm state.
77+
78+
.. code-block:: fortran
79+
80+
call micm%solve(time_step, state, solver_state, solver_stats, error)
81+
82+
This is much more efficient than calling the solver separately for each grid cell,
83+
especially when dealing with large numbers of atmospheric grid cells in climate or weather models.
84+
85+
To build and run the example, follow the same build instructions as in Chapter 1 and 2.
86+
The compilation command would be:
87+
88+
.. code-block:: bash
89+
90+
gfortran -o micm_multiple_grid_cells micm_multiple_grid_cells.F90 -I<MUSICA_DIR>/include -L<MUSICA_DIR>/lib64 -lmusica-fortran -lmusica -lstdc++ -lyaml-cpp
91+
92+
Assuming you name the executable ``micm_multiple_grid_cells``, you can run the program as follows:
93+
94+
.. code-block:: bash
95+
96+
$ ./micm_multiple_grid_cells
97+
Creating MICM solver with 3 grid cells...
98+
Creating State for multiple grid cells...
99+
Species in the mechanism:
100+
Species Name:A, Index: 1
101+
Species Name:B, Index: 2
102+
Species Name:C, Index: 5
103+
Species Name:D, Index: 3
104+
Species Name:E, Index: 4
105+
Species Name:F, Index: 6
106+
107+
Initial concentrations by grid cell:
108+
Grid Cell 1 (T= 273.0K):
109+
1.000 1.000 1.000 1.000 1.000 1.000
110+
Grid Cell 2 (T= 283.0K):
111+
2.000 2.000 2.000 2.000 2.000 2.000
112+
Grid Cell 3 (T= 293.0K):
113+
0.500 0.500 0.500 0.500 0.500 0.500
114+
115+
Solving for all grid cells simultaneously...
116+
117+
Final concentrations by grid cell:
118+
Grid Cell 1 (T= 273.0K):
119+
0.382 1.468 0.670 1.116 1.150 1.214
120+
Grid Cell 2 (T= 283.0K):
121+
0.764 2.936 1.340 2.232 2.300 2.428
122+
Grid Cell 3 (T= 293.0K):
123+
0.191 0.734 0.335 0.558 0.575 0.607
124+
125+
Solver completed successfully for all 3 grid cells!
126+
127+
The results differ in each grid cell due to the varying conditions.
128+

‎docs/source/tutorials/tutorials.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Fortran
1111
1. :ref:`installing MUSICA <chapter0>`
1212
2. :ref:`first Fortran MUSICA program <chapter1>`
1313
3. :ref:`box model example <chapter2>`
14+
4. :ref:`multiple grid cells box model <chapter3>`
1415

1516
.. toctree::
1617
:hidden:

‎fortran/test/fetch_content_integration/CMakeLists.txt‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ if (MUSICA_ENABLE_MICM)
4040
add_executable(test_micm_fortran_api test_micm_api.F90)
4141
add_executable(test_get_micm_version test_get_micm_version.F90)
4242
add_executable(test_micm_box_model test_micm_box_model.F90)
43+
add_executable(test_micm_multiple_grid_cells test_micm_multiple_grid_cells.F90)
4344

4445
target_link_libraries(test_micm_fortran_api
4546
PRIVATE
@@ -90,6 +91,23 @@ if (MUSICA_ENABLE_MICM)
9091
COMMAND $<TARGET_FILE:test_micm_box_model>
9192
WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
9293
)
94+
95+
target_link_libraries(test_micm_multiple_grid_cells
96+
PRIVATE
97+
musica::musica-fortran
98+
)
99+
100+
set_target_properties(test_micm_multiple_grid_cells
101+
PROPERTIES
102+
LINKER_LANGUAGE Fortran
103+
)
104+
105+
add_test(
106+
NAME test_micm_multiple_grid_cells
107+
COMMAND $<TARGET_FILE:test_micm_multiple_grid_cells>
108+
WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
109+
)
110+
93111
endif()
94112

95113
# API Test

‎fortran/test/fetch_content_integration/test_micm_box_model.F90‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,11 +4,17 @@ program test_micm_box_model
44
use, intrinsic :: ieee_arithmetic
55

66
use iso_fortran_env, only : real64
7-
use musica_util, only: error_t, string_t, mapping_t
7+
use musica_util, only: assert, error_t, string_t, mapping_t
88
use musica_micm, only: micm_t, solver_stats_t
99
use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
1010
use musica_state, only: conditions_t, state_t
1111

12+
#define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ )
13+
#define ASSERT_EQ( a, b ) call assert( a == b, __FILE__, __LINE__ )
14+
#define ASSERT_NE( a, b ) call assert( a /= b, __FILE__, __LINE__ )
15+
#define ASSERT_NEAR( a, b, tol ) call assert( abs(a - b) < abs(a + b) * tol, __FILE__, __LINE__ )
16+
#define ASSERT_GT( a, b ) call assert( a > b, __FILE__, __LINE__ )
17+
1218
implicit none
1319

1420
call box_model_arrays()
@@ -36,10 +42,12 @@ subroutine box_model_arrays()
3642
num_grid_cells = 1
3743

3844
write(*,*) "Creating MICM solver..."
39-
micm => micm_t(config_path, solver_type, error)
45+
micm => micm_t(config_path, solver_type, error)
46+
ASSERT( error%is_success() )
4047

4148
write(*,*) "Creating State..."
4249
state => micm%get_state(num_grid_cells, error)
50+
ASSERT( error%is_success() )
4351

4452
time_step = 200
4553

@@ -64,6 +72,8 @@ subroutine box_model_arrays()
6472

6573
write(*,*) "Solving starts..."
6674
call micm%solve(time_step, state, solver_state, solver_stats, error)
75+
ASSERT( error%is_success() )
76+
6777
write(*,*) "After solving, concentrations", state%concentrations
6878
deallocate( micm )
6979
deallocate( state )
@@ -92,9 +102,11 @@ subroutine box_model_c_ptrs()
92102

93103
write(*,*) "Creating MICM solver..."
94104
micm => micm_t(config_path, solver_type, error)
105+
ASSERT( error%is_success() )
95106

96107
write(*,*) "Creating State..."
97108
state => micm%get_state(num_grid_cells, error)
109+
ASSERT( error%is_success() )
98110

99111
state%conditions(1)%temperature = 273.0
100112
state%conditions(1)%pressure = 1.0e5
@@ -117,6 +129,7 @@ subroutine box_model_c_ptrs()
117129

118130
write(*,*) "Solving starts..."
119131
call micm%solve(time_step, state, solver_state, solver_stats, error)
132+
ASSERT( error%is_success() )
120133
write(*,*) "After solving, concentrations", state%concentrations
121134

122135
deallocate( micm )

‎fortran/test/fetch_content_integration/test_micm_multiple_grid_cells.F90‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
program test_micm_multiple_grid_cells
2+
3+
use, intrinsic :: iso_c_binding
4+
use, intrinsic :: ieee_arithmetic
5+
6+
use iso_fortran_env, only : real64
7+
use musica_util, only: assert, error_t, string_t, mapping_t
8+
use musica_micm, only: micm_t, solver_stats_t
9+
use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
10+
use musica_state, only: conditions_t, state_t
11+
12+
#define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ )
13+
#define ASSERT_EQ( a, b ) call assert( a == b, __FILE__, __LINE__ )
14+
#define ASSERT_NE( a, b ) call assert( a /= b, __FILE__, __LINE__ )
15+
#define ASSERT_NEAR( a, b, tol ) call assert( abs(a - b) < abs(a + b) * tol, __FILE__, __LINE__ )
16+
#define ASSERT_GT( a, b ) call assert( a > b, __FILE__, __LINE__ )
17+
18+
implicit none
19+
20+
call multiple_grid_cells_example()
21+
22+
contains
23+
24+
!> Runs a multiple grid cells box model using the MICM solver
25+
subroutine multiple_grid_cells_example()
26+
27+
character(len=256) :: config_path
28+
integer :: solver_type
29+
integer :: num_grid_cells
30+
real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1
31+
real(real64) :: time_step
32+
type(string_t) :: solver_state
33+
type(solver_stats_t) :: solver_stats
34+
type(error_t) :: error
35+
type(micm_t), pointer :: micm
36+
type(state_t), pointer :: state
37+
integer :: i, j, cell_id
38+
39+
config_path = "configs/v0/analytical"
40+
solver_type = RosenbrockStandardOrder
41+
num_grid_cells = 3 ! Use 3 grid cells to demonstrate multiple cells
42+
43+
write(*,*) "Creating MICM solver with", num_grid_cells, "grid cells..."
44+
micm => micm_t(config_path, solver_type, error)
45+
ASSERT( error%is_success() )
46+
47+
write(*,*) "Creating State for multiple grid cells..."
48+
state => micm%get_state(num_grid_cells, error)
49+
ASSERT( error%is_success() )
50+
51+
time_step = 200
52+
53+
! Set up conditions for each grid cell
54+
do cell_id = 1, num_grid_cells
55+
state%conditions(cell_id)%temperature = 273.0 + (cell_id - 1) * 10.0 ! Different temperatures
56+
state%conditions(cell_id)%pressure = 1.0e5
57+
state%conditions(cell_id)%air_density = state%conditions(cell_id)%pressure / &
58+
(GAS_CONSTANT * state%conditions(cell_id)%temperature)
59+
end do
60+
61+
! Set initial concentrations for each grid cell
62+
associate( cell_stride => state%species_strides%grid_cell, &
63+
var_stride => state%species_strides%variable )
64+
do cell_id = 1, num_grid_cells
65+
! Set different initial concentrations for each cell
66+
do i = 1, 6 ! Assuming 6 species
67+
! Cell 1: all species start at 1.0
68+
! Cell 2: all species start at 2.0
69+
! Cell 3: all species start at 0.5
70+
if (cell_id == 1) then
71+
state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 1.0
72+
else if (cell_id == 2) then
73+
state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 2.0
74+
else
75+
state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 0.5
76+
end if
77+
end do
78+
end do
79+
end associate
80+
81+
! Set rate parameters for each grid cell
82+
associate( cell_stride => state%rate_parameters_strides%grid_cell, &
83+
var_stride => state%rate_parameters_strides%variable )
84+
do cell_id = 1, num_grid_cells
85+
! Set the same rate parameters for all cells
86+
state%rate_parameters(1 + (cell_id - 1) * cell_stride + 0 * var_stride) = 0.001
87+
state%rate_parameters(1 + (cell_id - 1) * cell_stride + 1 * var_stride) = 0.002
88+
end do
89+
end associate
90+
91+
! Print species information
92+
write(*,*) "Species in the mechanism:"
93+
do i = 1, state%species_ordering%size()
94+
write(*,*) "Species Name:", trim(state%species_ordering%name(i)), &
95+
", Index:", state%species_ordering%index(i)
96+
end do
97+
98+
! Print initial concentrations for each grid cell
99+
write(*,*) ""
100+
write(*,*) "Initial concentrations by grid cell:"
101+
associate( cell_stride => state%species_strides%grid_cell, &
102+
var_stride => state%species_strides%variable )
103+
do cell_id = 1, num_grid_cells
104+
write(*,'(A,I0,A,F6.1,A)') "Grid Cell ", cell_id, " (T=", &
105+
state%conditions(cell_id)%temperature, "K):"
106+
do i = 1, 6
107+
write(*,'(A,F8.3)', advance='no') " ", &
108+
state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride)
109+
end do
110+
write(*,*) ""
111+
end do
112+
end associate
113+
114+
write(*,*) ""
115+
write(*,*) "Solving for all grid cells simultaneously..."
116+
call micm%solve(time_step, state, solver_state, solver_stats, error)
117+
ASSERT( error%is_success() )
118+
119+
! Print final concentrations for each grid cell
120+
write(*,*) ""
121+
write(*,*) "Final concentrations by grid cell:"
122+
associate( cell_stride => state%species_strides%grid_cell, &
123+
var_stride => state%species_strides%variable )
124+
do cell_id = 1, num_grid_cells
125+
write(*,'(A,I0,A,F6.1,A)') "Grid Cell ", cell_id, " (T=", &
126+
state%conditions(cell_id)%temperature, "K):"
127+
do i = 1, 6
128+
write(*,'(A,F8.3)', advance='no') " ", &
129+
state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride)
130+
end do
131+
write(*,*) ""
132+
end do
133+
end associate
134+
135+
write(*,*) ""
136+
write(*,*) "Solver completed successfully for all", num_grid_cells, "grid cells!"
137+
138+
deallocate( micm )
139+
deallocate( state )
140+
end subroutine multiple_grid_cells_example
141+
142+
end program test_micm_multiple_grid_cells

0 commit comments

Comments
 (0)