Merge pull request #145 from NCAR/83-solver-result

Return MICM solver state details from MicmSolve()

docker/Dockerfile

@@ -1,5 +1,7 @@
 FROM fedora:35
 
+ARG BUILD_TYPE=Release
+
 RUN dnf -y update \
     && dnf -y install \
         cmake \
@@ -23,6 +25,7 @@ COPY . musica
 RUN cd musica \
     && cmake -S . \
              -B build \
+             -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \
     && cd build \
     && make install -j 8 
 

docker/Dockerfile.fortran-gcc.integration

@@ -1,6 +1,7 @@
 FROM fedora:35
 
 ARG MUSICA_GIT_TAG=main
+ARG BUILD_TYPE=Release
 
 RUN dnf -y update \
   && dnf -y install \
@@ -37,20 +38,18 @@ COPY . musica
 RUN cd musica \
     && cmake -S . \
              -B build \
-             -D CMAKE_BUILD_TYPE=Release \
+             -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \
     && cd build \
     && make install -j
 
 # Build and install MUSICA
 RUN cd musica/fortran/test/fetch_content_integration \
     && mkdir build && cd build \
     && cmake .. \
-        -D CMAKE_BUILD_TYPE=Release \
+        -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \
         -D MUSICA_GIT_TAG=${MUSICA_GIT_TAG} \
         -D MUSICA_ENABLE_MICM=ON \
         -D MUSICA_ENABLE_TUVX=OFF \
     && make -j
 
-WORKDIR musica/fortran/test/fetch_content_integration/build
-RUN cp -r /musica/build/_deps/tuvx-src/examples/ .
-RUN cp -r /musica/build/_deps/tuvx-src/data/ .
+WORKDIR musica/fortran/test/fetch_content_integration/build

docker/Dockerfile.fortran-nvhpc

@@ -68,6 +68,4 @@ RUN cd musica/fortran/test/fetch_content_integration \
         -D MUSICA_GIT_TAG=${MUSICA_GIT_TAG} \
     && make -j
 
-WORKDIR musica/fortran/test/fetch_content_integration/build
-RUN cp -r /musica/build/_deps/tuvx-src/examples/ .
-RUN cp -r /musica/build/_deps/tuvx-src/data/ .
+WORKDIR musica/fortran/test/fetch_content_integration/build

fortran/micm.F90

@@ -3,21 +3,29 @@
 !
 module musica_micm
 #define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ )
-   use iso_c_binding, only: c_ptr, c_char, c_int, c_bool, c_double, c_null_char, &
+   use iso_c_binding, only: c_ptr, c_char, c_int, c_int64_t, c_bool, c_double, c_null_char, &
                             c_size_t, c_f_pointer, c_funptr, c_null_ptr, c_associated
-   use musica_util, only: assert, mapping_t
+   use iso_fortran_env, only: int64
+   use musica_util, only: assert, mapping_t, string_t, string_t_c
    implicit none
 
-   public :: get_micm_version, micm_t
+   public :: micm_t, solver_stats_t, get_micm_version
    private
 
+   !> Wrapper for c solver stats
+   type, bind(c) :: solver_stats_t_c
+      integer(c_int64_t) :: function_calls_ = 0_c_int64_t
+      integer(c_int64_t) :: jacobian_updates_ = 0_c_int64_t
+      integer(c_int64_t) :: number_of_steps_ = 0_c_int64_t
+      integer(c_int64_t) :: accepted_ = 0_c_int64_t
+      integer(c_int64_t) :: rejected_ = 0_c_int64_t
+      integer(c_int64_t) :: decompositions_ = 0_c_int64_t
+      integer(c_int64_t) :: solves_ = 0_c_int64_t
+      integer(c_int64_t) :: singular_ = 0_c_int64_t
+      real(c_double)     :: final_time_ = 0._c_double
+   end type solver_stats_t_c
 
    interface
-      function get_micm_version_c() bind(C, name="get_micm_version")
-         use musica_util, only: string_t_c
-         type(string_t_c) :: get_micm_version_c
-      end function get_micm_version_c
-
       function create_micm_c(config_path, error) bind(C, name="CreateMicm")
          use musica_util, only: error_t_c
          import c_ptr, c_int, c_char
@@ -34,9 +42,10 @@ subroutine delete_micm_c(micm, error) bind(C, name="DeleteMicm")
       end subroutine delete_micm_c
 
       subroutine micm_solve_c(micm, time_step, temperature, pressure, air_density, num_concentrations, concentrations, &
-                              num_user_defined_reaction_rates, user_defined_reaction_rates, error) bind(C, name="MicmSolve")
-         use musica_util, only: error_t_c
-         import c_ptr, c_double, c_int
+                  num_user_defined_reaction_rates, user_defined_reaction_rates, solver_state, solver_stats, error) &
+                  bind(C, name="MicmSolve")
+         use musica_util, only: string_t_c, error_t_c
+         import c_ptr, c_double, c_int, solver_stats_t_c
          type(c_ptr), value, intent(in)         :: micm
          real(kind=c_double), value, intent(in) :: time_step
          real(kind=c_double), value, intent(in) :: temperature
@@ -46,9 +55,16 @@ subroutine micm_solve_c(micm, time_step, temperature, pressure, air_density, num
          real(kind=c_double), intent(inout)     :: concentrations(num_concentrations)
          integer(kind=c_int), value, intent(in) :: num_user_defined_reaction_rates
          real(kind=c_double), intent(inout)     :: user_defined_reaction_rates(num_user_defined_reaction_rates)
+         type(string_t_c), intent(out)          :: solver_state
+         type(solver_stats_t_c), intent(out)    :: solver_stats
          type(error_t_c), intent(inout)         :: error
       end subroutine micm_solve_c
 
+      function get_micm_version_c() bind(C, name="MicmVersion")
+         use musica_util, only: string_t_c
+         type(string_t_c) :: get_micm_version_c
+      end function get_micm_version_c
+
       function get_species_property_string_c(micm, species_name, property_name, error) bind(c, name="GetSpeciesPropertyString")
          use musica_util, only: error_t_c, string_t_c
          import c_ptr, c_char
@@ -129,6 +145,33 @@ end subroutine delete_mappings_c
       procedure constructor
    end interface micm_t
 
+   !> Solver stats type
+   type :: solver_stats_t
+      integer(int64) :: function_calls_
+      integer(int64) :: jacobian_updates_
+      integer(int64) :: number_of_steps_
+      integer(int64) :: accepted_
+      integer(int64) :: rejected_
+      integer(int64) :: decompositions_
+      integer(int64) :: solves_
+      integer(int64) :: singular_
+      real           :: final_time_
+   contains
+      procedure :: function_calls => solver_stats_t_function_calls
+      procedure :: jacobian_updates => solver_stats_t_jacobian_updates
+      procedure :: number_of_steps => solver_stats_t_number_of_steps
+      procedure :: accepted => solver_stats_t_accepted
+      procedure :: rejected => solver_stats_t_rejected
+      procedure :: decompositions => solver_stats_t_decompositions
+      procedure :: solves => solver_stats_t_solves
+      procedure :: singular => solver_stats_t_singular
+      procedure :: final_time => solver_stats_t_final_time
+   end type solver_stats_t
+
+   interface solver_stats_t
+      procedure solver_stats_t_constructor
+   end interface solver_stats_t
+
 contains
 
    function get_micm_version() result(value)
@@ -190,26 +233,142 @@ function constructor(config_path, error)  result( this )
    end function constructor
 
    subroutine solve(this, time_step, temperature, pressure, air_density, num_concentrations, concentrations, &
-                     num_user_defined_reaction_rates, user_defined_reaction_rates, error)
-      use musica_util, only: error_t_c, error_t
-      class(micm_t)                  :: this
-      real(c_double),  intent(in)    :: time_step
-      real(c_double),  intent(in)    :: temperature
-      real(c_double),  intent(in)    :: pressure
-      real(c_double),  intent(in)    :: air_density
-      integer(c_int),  intent(in)    :: num_concentrations
-      real(c_double),  intent(inout) :: concentrations(*)
-      integer(c_int),  intent(in)    :: num_user_defined_reaction_rates
-      real(c_double),  intent(inout) :: user_defined_reaction_rates(*)
-      type(error_t),   intent(out)   :: error
-
-      type(error_t_c)                :: error_c
-      call micm_solve_c(this%ptr, time_step, temperature, pressure, air_density, num_concentrations, &
-                        concentrations, num_user_defined_reaction_rates, user_defined_reaction_rates, &
-                        error_c)
+               num_user_defined_reaction_rates, user_defined_reaction_rates, solver_state, solver_stats, error)
+      use musica_util, only: string_t, string_t_c, error_t_c, error_t
+      class(micm_t)                       :: this
+      real(c_double),       intent(in)    :: time_step
+      real(c_double),       intent(in)    :: temperature
+      real(c_double),       intent(in)    :: pressure
+      real(c_double),       intent(in)    :: air_density
+      integer(c_int),       intent(in)    :: num_concentrations
+      real(c_double),       intent(inout) :: concentrations(*)
+      integer(c_int),       intent(in)    :: num_user_defined_reaction_rates
+      real(c_double),       intent(inout) :: 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, time_step, temperature, pressure, air_density, num_concentrations, concentrations, &
+            num_user_defined_reaction_rates, 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
 
+   !> 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 )
+      use iso_fortran_env, only: int64
+      use musica_util, only: string_t
+      type(solver_stats_t_c), intent(inout) :: c_solver_stats
+      type(solver_stats_t)                  :: new_solver_stats
+
+      new_solver_stats%function_calls_ = c_solver_stats%function_calls_
+      new_solver_stats%jacobian_updates_ = c_solver_stats%jacobian_updates_
+      new_solver_stats%number_of_steps_ = c_solver_stats%number_of_steps_
+      new_solver_stats%accepted_ = c_solver_stats%accepted_
+      new_solver_stats%rejected_ = c_solver_stats%rejected_
+      new_solver_stats%decompositions_ = c_solver_stats%decompositions_
+      new_solver_stats%solves_ = c_solver_stats%solves_
+      new_solver_stats%singular_ = c_solver_stats%singular_
+      new_solver_stats%final_time_ = real( c_solver_stats%final_time_ )
+
+   end function solver_stats_t_constructor
+
+   !> Get the number of forcing function calls
+   function solver_stats_t_function_calls( this ) result( function_calls )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: function_calls
+
+      function_calls = this%function_calls_
+
+   end function solver_stats_t_function_calls
+
+   !> Get the number of jacobian function calls
+   function solver_stats_t_jacobian_updates( this ) result( jacobian_updates )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: jacobian_updates
+
+      jacobian_updates = this%jacobian_updates_
+
+   end function solver_stats_t_jacobian_updates
+
+   !> Get the total number of internal time steps taken
+   function solver_stats_t_number_of_steps( this ) result( number_of_steps )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: number_of_steps
+
+      number_of_steps = this%number_of_steps_
+
+   end function solver_stats_t_number_of_steps
+
+   !> Get the number of accepted integrations
+   function solver_stats_t_accepted( this ) result( accepted )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: accepted
+
+      accepted = this%accepted_
+
+   end function solver_stats_t_accepted
+
+   !> Get the number of rejected integrations
+   function solver_stats_t_rejected( this ) result( rejected )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: rejected
+
+      rejected = this%rejected_
+
+   end function solver_stats_t_rejected
+
+   !> Get the number of LU decompositions
+   function solver_stats_t_decompositions( this ) result( decompositions )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: decompositions
+
+      decompositions = this%decompositions_
+
+   end function solver_stats_t_decompositions
+
+   !> Get the number of linear solves
+   function solver_stats_t_solves( this ) result( solves )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: solves
+
+      solves = this%solves_
+
+   end function solver_stats_t_solves
+
+   !> Get the number of times a singular matrix is detected
+   function solver_stats_t_singular( this ) result( singular )
+      use iso_fortran_env, only: int64
+      class(solver_stats_t), intent(in) :: this
+      integer(int64)                    :: singular
+
+      singular = this%function_calls_
+
+   end function solver_stats_t_singular
+
+   !> Get the final time the solver iterated to
+   function solver_stats_t_final_time( this ) result( final_time )
+      class(solver_stats_t), intent(in) :: this
+      real                              :: final_time
+
+      final_time = this%final_time_
+
+   end function solver_stats_t_final_time
+
    function get_species_property_string(this, species_name, property_name, error) result(value)
       use musica_util, only: error_t_c, error_t, string_t, string_t_c, to_c_string
       class(micm_t), intent(inout) :: this
@@ -276,4 +435,4 @@ subroutine finalize(this)
       ASSERT(error%is_success())
    end subroutine finalize
 
-end module musica_micm
+end module musica_micm

fortran/test/fetch_content_integration/CMakeLists.txt

@@ -89,4 +89,14 @@ if (MUSICA_ENABLE_TUVX)
     COMMAND $<TARGET_FILE:test_tuvx_fortran_api>
     WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
   )
+
+  add_custom_target(
+    copy_tuvx_examples_dir ALL ${CMAKE_COMMAND} -E copy_directory
+    ${CMAKE_CURRENT_SOURCE_DIR}/../../../build/_deps/tuvx-src/examples ${CMAKE_BINARY_DIR}/examples
+  )
+
+  add_custom_target(
+    copy_tuvx_data_dir ALL ${CMAKE_COMMAND} -E copy_directory
+    ${CMAKE_CURRENT_SOURCE_DIR}/../../../build/_deps/tuvx-src/data ${CMAKE_BINARY_DIR}/data
+  )
 endif()

fortran/test/fetch_content_integration/test_micm_api.F90

@@ -5,7 +5,7 @@ program test_micm_api
 
   use, intrinsic :: iso_c_binding
   use, intrinsic :: ieee_arithmetic
-  use musica_micm, only: micm_t, get_micm_version
+  use musica_micm, only: micm_t, solver_stats_t, get_micm_version
   use musica_util, only: assert, error_t, mapping_t, string_t
 
 #include "micm/util/error.hpp"
@@ -31,15 +31,17 @@ subroutine test_api()
     integer(c_int)                :: num_concentrations, num_user_defined_reaction_rates
     real(c_double), dimension(5)  :: concentrations 
     real(c_double), dimension(3)  :: user_defined_reaction_rates 
-    integer                       :: i
     character(len=256)            :: config_path
     character(len=:), allocatable :: string_value
     real(c_double)                :: double_value
     integer(c_int)                :: int_value
     logical(c_bool)               :: bool_value
+    type(string_t)                :: solver_state
+    type(solver_stats_t)          :: solver_stats
     type(error_t)                 :: error
     real(c_double), parameter     :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1
-
+    integer                       :: i
+
     time_step = 200
     temperature = 272.5
     pressure = 101253.4
@@ -71,11 +73,21 @@ subroutine test_api()
     write(*,*) "[test micm fort api] Initial concentrations", concentrations
 
     write(*,*) "[test micm fort api] Solving starts..."
-    call micm%solve(time_step, temperature, pressure, air_density, num_concentrations, concentrations, &
-                    num_user_defined_reaction_rates, user_defined_reaction_rates, error)
+    call micm%solve(time_step, temperature, pressure,  air_density, num_concentrations, concentrations, &
+        num_user_defined_reaction_rates, user_defined_reaction_rates, solver_state, solver_stats, error)
     ASSERT( error%is_success() )
 
-    write(*,*) "[test micm fort api] After solving, concentrations", concentrations
+    write(*,*) "[test micm fort api] After solving, concentrations: ", concentrations
+    write(*,*) "[test micm fort api] Solver state: ", solver_state%get_char_array()
+    write(*,*) "[test micm fort api] Function calls: ", solver_stats%function_calls()
+    write(*,*) "[test micm fort api] Jacobian updates: ", solver_stats%jacobian_updates()
+    write(*,*) "[test micm fort api] Number of steps: ", solver_stats%number_of_steps()
+    write(*,*) "[test micm fort api] Accepted: ", solver_stats%accepted()
+    write(*,*) "[test micm fort api] Rejected: ", solver_stats%rejected()
+    write(*,*) "[test micm fort api] Decompositions: ", solver_stats%decompositions()
+    write(*,*) "[test micm fort api] Solves: ", solver_stats%solves()
+    write(*,*) "[test micm fort api] Singular: ", solver_stats%singular()
+    write(*,*) "[test micm fort api] Final time: ", solver_stats%final_time()
 
     string_value = micm%get_species_property_string( "O3", "__long name", error )
     ASSERT( error%is_success() )