add fortran MICM solve function that takes c pointers

NCAR · Nov 4, 2024 · 10fbd0d · 10fbd0d
1 parent ff12c14
commit 10fbd0d
Show file tree

Hide file tree

Showing 10 changed files with 575 additions and 126 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -30,7 +30,7 @@ option(MUSICA_BUILD_DOCS "Build the documentation" OFF)
 option(MUSICA_ENABLE_MICM "Enable MICM" ON)
 option(MUSICA_ENABLE_TUVX "Enable TUV-x" ON)
 
-set(MUSICA_SET_MICM_VECTOR_MATRIX_SIZE "1" CACHE STRING "Set MICM vector-ordered matrix dimension")
+set(MUSICA_SET_MICM_VECTOR_MATRIX_SIZE "4" CACHE STRING "Set MICM vector-ordered matrix dimension")
 
 cmake_dependent_option(
   MUSICA_ENABLE_PYTHON_LIBRARY "Adds pybind11, a lightweight header-only library that exposes C++ types in Python and vice versa" OFF "MUSICA_BUILD_C_CXX_INTERFACE" OFF)

diff --git a/fortran/micm.F90 b/fortran/micm.F90
@@ -135,10 +135,14 @@ end function get_user_defined_reaction_rates_ordering_c
   type :: micm_t
     type(mappings_t), pointer :: species_ordering => null()
     type(mappings_t), pointer :: user_defined_reaction_rates => null()
-    type(c_ptr), private         :: ptr = c_null_ptr
+    type(c_ptr), private      :: ptr = c_null_ptr
+    integer,     private      :: number_of_grid_cells = 0
+    integer,     private      :: solver_type = 0
   contains
     ! Solve the chemical system
-    procedure :: solve
+    procedure, private :: solve_arrays
+    procedure, private :: solve_c_ptrs
+    generic :: solve => solve_arrays, solve_c_ptrs
     ! Get species properties
     procedure :: get_species_property_string
     procedure :: get_species_property_double
@@ -208,6 +212,8 @@ function constructor(config_path, solver_type, num_grid_cells, error)  result( t
     end do
     c_config_path(n+1) = c_null_char
 
+    this%number_of_grid_cells = num_grid_cells
+    this%solver_type = solver_type
     this%ptr = create_micm_c( c_config_path, int(solver_type, kind=c_int), &
                               int(num_grid_cells, kind=c_int), error_c )
     error = error_t(error_c)
@@ -236,8 +242,12 @@ function constructor(config_path, solver_type, num_grid_cells, error)  result( t
 
   end function constructor
 
-  subroutine solve(this, time_step, temperature, pressure, air_density, concentrations, &
-              user_defined_reaction_rates, solver_state, solver_stats, error)
+  !> Solves the chemical system
+  !!
+  !! This function accepts fortran arrays and checks their sizes
+  !! against the number of grid cells and the species/rate parameter ordering.
+  subroutine solve_arrays(this, time_step, temperature, pressure, air_density, &
+      concentrations, user_defined_reaction_rates, solver_state, solver_stats, error)
     use iso_c_binding, only: c_loc
     use iso_fortran_env, only: real64
     use musica_util, only: string_t, string_t_c, error_t_c, error_t
@@ -246,8 +256,8 @@ subroutine solve(this, time_step, temperature, pressure, air_density, concentrat
     real(real64), target, intent(in)    :: temperature(:)
     real(real64), target, intent(in)    :: pressure(:)
     real(real64), target, intent(in)    :: air_density(:)
-    real(real64), target, intent(inout) :: concentrations(:)
-    real(real64), target, intent(in)    :: user_defined_reaction_rates(:)
+    real(real64), target, intent(inout) :: concentrations(:,:)
+    real(real64), target, intent(in)    :: user_defined_reaction_rates(:,:)
     type(string_t),       intent(out)   :: solver_state
     type(solver_stats_t), intent(out)   :: solver_stats
     type(error_t),        intent(out)   :: error
@@ -256,6 +266,54 @@ subroutine solve(this, time_step, temperature, pressure, air_density, concentrat
     type(solver_stats_t_c) :: solver_stats_c
     type(error_t_c)        :: error_c
 
+    if (size(temperature) .ne. this%number_of_grid_cells) then
+        error = error_t(1, "MICM_SOLVE", "Temperature array size does not match number of grid cells")
+        return
+    end if
+    if (size(pressure) .ne. this%number_of_grid_cells) then
+        error = error_t(1, "MICM_SOLVE", "Pressure array size does not match number of grid cells")
+        return
+    end if
+    if (size(air_density) .ne. this%number_of_grid_cells) then
+        error = error_t(1, "MICM_SOLVE", "Air density array size does not match number of grid cells")
+        return
+    end if
+    if (this%solver_type .eq. Rosenbrock .or. this%solver_type .eq. BackwardEuler) then
+        if (size(concentrations, 1) .ne. this%number_of_grid_cells) then
+            error = error_t(1, "MICM_SOLVE", "Concentrations array dimension 1 does not match number of grid cells")
+            return
+        end if
+        if (size(concentrations, 2) .ne. this%species_ordering%size()) then
+            error = error_t(1, "MICM_SOLVE", "Concentrations array dimension 2 does not match species ordering")
+            return
+        end if
+        if (size(user_defined_reaction_rates, 1) .ne. this%number_of_grid_cells) then
+            error = error_t(1, "MICM_SOLVE", "User defined reaction rates array dimension 1 does not match number of grid cells")
+            return
+        end if
+        if (size(user_defined_reaction_rates, 2) .ne. this%user_defined_reaction_rates%size()) then
+            error = error_t(1, "MICM_SOLVE", "User defined reaction rates array dimension 2 does not match user defined reaction rates ordering")
+            return
+        end if
+    else
+        if (size(concentrations, 1) .ne. this%species_ordering%size()) then
+            error = error_t(1, "MICM_SOLVE", "Concentrations array dimension 1 does not match species ordering")
+            return
+        end if
+        if (size(concentrations, 2) .ne. this%number_of_grid_cells) then
+            error = error_t(1, "MICM_SOLVE", "Concentrations array dimension 2 does not match number of grid cells")
+            return
+        end if
+        if (size(user_defined_reaction_rates, 1) .ne. this%user_defined_reaction_rates%size()) then
+            error = error_t(1, "MICM_SOLVE", "User defined reaction rates array dimension 1 does not match user defined reaction rates ordering")
+            return
+        end if
+        if (size(user_defined_reaction_rates, 2) .ne. this%number_of_grid_cells) then
+            error = error_t(1, "MICM_SOLVE", "User defined reaction rates array dimension 2 does not match number of grid cells")
+            return
+        end if
+    end if
+
     call micm_solve_c(this%ptr, real(time_step, kind=c_double), c_loc(temperature), &
                       c_loc(pressure), c_loc(air_density), c_loc(concentrations), &
                       c_loc(user_defined_reaction_rates), &
@@ -265,7 +323,40 @@ subroutine solve(this, time_step, temperature, pressure, air_density, concentrat
     solver_stats = solver_stats_t(solver_stats_c)
     error = error_t(error_c)
 
-  end subroutine solve
+  end subroutine solve_arrays
+
+  !> Solves the chemical system
+  !!
+  !! This function accepts c pointers and does not check their sizes.
+  !! The user is responsible for ensuring the sizes are correct.
+  subroutine solve_c_ptrs(this, time_step, temperature, pressure, air_density, &
+      concentrations, user_defined_reaction_rates, solver_state, solver_stats, error)
+    use iso_fortran_env, only: real64
+    use musica_util, only: string_t, string_t_c, error_t_c, error_t
+    class(micm_t),        intent(in)    :: this
+    real(real64),         intent(in)    :: time_step
+    type(c_ptr),          intent(in)    :: temperature
+    type(c_ptr),          intent(in)    :: pressure
+    type(c_ptr),          intent(in)    :: air_density
+    type(c_ptr),          intent(in)    :: concentrations
+    type(c_ptr),          intent(in)    :: user_defined_reaction_rates
+    type(string_t),       intent(out)   :: solver_state
+    type(solver_stats_t), intent(out)   :: solver_stats
+    type(error_t),        intent(out)   :: error
+
+    type(string_t_c)       :: solver_state_c
+    type(solver_stats_t_c) :: solver_stats_c
+    type(error_t_c)        :: error_c
+
+    call micm_solve_c(this%ptr, real(time_step, kind=c_double), temperature, pressure, &
+                      air_density, concentrations, user_defined_reaction_rates, &
+                      solver_state_c, solver_stats_c, error_c)
+
+    solver_state = string_t(solver_state_c)
+    solver_stats = solver_stats_t(solver_stats_c)
+    error = error_t(error_c)
+
+  end subroutine solve_c_ptrs
 
   !> Constructor for solver_stats_t object that takes ownership of solver_stats_t_c
   function solver_stats_t_constructor( c_solver_stats ) result( new_solver_stats )

diff --git a/fortran/test/fetch_content_integration/CMakeLists.txt b/fortran/test/fetch_content_integration/CMakeLists.txt
@@ -48,6 +48,8 @@ if (MUSICA_ENABLE_MICM)
       LINKER_LANGUAGE Fortran
   )
 
+  target_compile_definitions(test_micm_fortran_api PUBLIC MICM_VECTOR_MATRIX_SIZE=${MUSICA_SET_MICM_VECTOR_MATRIX_SIZE})
+
   add_test(
     NAME test_micm_fortran_api
     COMMAND $<TARGET_FILE:test_micm_fortran_api>