Add TUV-x run function (#223)

* pass maps to TUV-x internal constructor

* add TUV-x run function

* remove unneeded valgrind suppressions

* add tuvx run test with host data

* add TUV-x run function to fortran wrapper

* add TUV-x fortran test

* fix dependency bug

* address reviewer comments

fortran/test/unit/tuvx.F90

@@ -7,10 +7,29 @@
 !> Test module for the tuvx connection
 program test_tuvx_connection
   use musica_assert
+  use musica_util, only : dk => musica_dk
 
   implicit none
 
+  integer, parameter :: number_of_layers = 3
+  integer, parameter :: number_of_wavelengths = 5
+  integer, parameter :: number_of_reactions = 2
+
+  real(dk), parameter :: expected_photolysis_rate_constants(4,2) = reshape( [ &
+    8.91393763338872e-28_dk, 1.64258192104497e-20_dk, 8.48391527327371e-14_dk, 9.87420948924703e-08_dk, &
+    2.49575956372508e-27_dk, 4.58686176250519e-20_dk, 2.22679622672858e-13_dk, 2.29392676897831e-07_dk ], [4,2] )
+  real(dk), parameter :: expected_heating_rates(4,2) = reshape( [ &
+    1.12394047546984e-46_dk, 2.04518267143613e-39_dk, 7.44349752571804e-33_dk, 5.42628100199216e-28_dk, &
+    5.14970120496081e-46_dk, 9.37067648164478e-39_dk, 3.41659389501112e-32_dk, 5.46672356294259e-27_dk ], [4,2] )
+
+#define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ )
+#define ASSERT_EQ( a, b ) call assert( a == b, __FILE__, __LINE__ )
+#define ASSERT_NE( a, b ) call assert( a /= b, __FILE__, __LINE__ )
+#define ASSERT_NEAR( a, b, tol ) call assert( abs(a - b) < abs(a + b) * tol, __FILE__, __LINE__ )
+
   call test_tuvx()
+  call test_tuvx_fixed()
+  call test_tuvx_data_from_host()
 
 contains
 
@@ -41,6 +60,208 @@ subroutine test_tuvx( )
 
   end subroutine test_tuvx
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Test TUV-x with a fixed configuration
+  subroutine test_tuvx_fixed( )
+
+    use musica_util,              only : assert, error_t
+    use musica_tuvx_grid,         only : grid_t
+    use musica_tuvx_grid_map,     only : grid_map_t
+    use musica_tuvx_profile,      only : profile_t
+    use musica_tuvx_profile_map,  only : profile_map_t
+    use musica_tuvx_radiator_map, only : radiator_map_t
+    use musica_tuvx,              only : tuvx_t
+
+    character(len=*), parameter :: my_name = "TUV-x fixed configuration tests"
+    character(len=*), parameter :: config_path = "test/data/tuvx/fixed/config.json"
+    type(grid_map_t),     pointer :: grids
+    type(profile_map_t),  pointer :: profiles
+    type(radiator_map_t), pointer :: radiators
+    type(tuvx_t),         pointer :: tuvx
+    type(error_t) :: error
+
+    real(dk) :: photo_rate_consts(number_of_layers + 1, number_of_reactions)
+    real(dk) :: heating_rates(number_of_layers + 1, number_of_reactions)
+    integer :: i, j
+    real(dk) :: a, b
+
+    grids => grid_map_t( error )
+    ASSERT( error%is_success() )
+    profiles => profile_map_t( error )
+    ASSERT( error%is_success() )
+    radiators => radiator_map_t( error )
+    ASSERT( error%is_success() )
+    tuvx => tuvx_t( config_path, grids, profiles, radiators, error )
+    ASSERT( error%is_success() )
+    call tuvx%run( 0.1_dk, 1.1_dk, photo_rate_consts, heating_rates, error )
+    ASSERT( error%is_success() )
+    do i = 1, number_of_reactions
+      do j = 1, number_of_layers + 1
+        a = photo_rate_consts(j,i)
+        b = expected_photolysis_rate_constants(j,i)
+        ASSERT_NEAR( a, b, 1.0e-5_dk )
+        a = heating_rates(j,i)
+        b = expected_heating_rates(j,i)
+        ASSERT_NEAR( a, b, 1.0e-5_dk )
+      end do
+    end do
+    deallocate( tuvx )
+    deallocate( radiators )
+    deallocate( profiles )
+    deallocate( grids )
+
+  end subroutine test_tuvx_fixed
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Test TUV-x with host-supplied grids and profiles
+  subroutine test_tuvx_data_from_host( )
+
+    use musica_util,              only : assert, error_t
+    use musica_tuvx_grid,         only : grid_t
+    use musica_tuvx_grid_map,     only : grid_map_t
+    use musica_tuvx_profile,      only : profile_t
+    use musica_tuvx_profile_map,  only : profile_map_t
+    use musica_tuvx_radiator_map, only : radiator_map_t
+    use musica_tuvx,              only : tuvx_t
+
+    character(len=*), parameter :: my_name = "TUV-x from-host configuration tests"
+    character(len=*), parameter :: config_path = "test/data/tuvx/from_host/config.json"
+    type(grid_t),         pointer :: heights, wavelengths
+    type(grid_map_t),     pointer :: grids
+    type(profile_t),      pointer :: temperatures
+    type(profile_map_t),  pointer :: profiles
+    type(radiator_map_t), pointer :: radiators
+    type(tuvx_t),         pointer :: tuvx
+    type(error_t) :: error
+
+    real(dk) :: photo_rate_consts(number_of_layers + 1, number_of_reactions)
+    real(dk) :: heating_rates(number_of_layers + 1, number_of_reactions)
+    integer :: i, j
+    real(dk) :: a, b
+    real(dk), allocatable :: temp_vals(:)
+
+    heights => grid_t( "height", "km", 3, error )
+    ASSERT( error%is_success() )
+    call heights%set_edges( [ 0.0_dk, 1.0_dk, 2.0_dk, 3.0_dk ], error )
+    ASSERT( error%is_success() )
+    call heights%set_midpoints( [ 0.5_dk, 1.5_dk, 2.5_dk ], error )
+    wavelengths => grid_t( "wavelength", "nm", 5, error )
+    ASSERT( error%is_success() )
+    call wavelengths%set_edges( [ 300.0_dk, 400.0_dk, 500.0_dk, 600.0_dk, 700.0_dk, 800.0_dk ], error )
+    ASSERT( error%is_success() )
+    call wavelengths%set_midpoints( [ 350.0_dk, 450.0_dk, 550.0_dk, 650.0_dk, 750.0_dk ], error )
+    ASSERT( error%is_success() )
+    grids => grid_map_t( error )
+    ASSERT( error%is_success() )
+    call grids%add( heights, error )
+    ASSERT( error%is_success() )
+    call grids%add( wavelengths, error )
+    ASSERT( error%is_success() )
+    temperatures => profile_t( "temperature", "K", heights, error )
+    ASSERT( error%is_success() )
+    profiles => profile_map_t( error )
+    ASSERT( error%is_success() )
+    call profiles%add( temperatures, error )
+    ASSERT( error%is_success() )
+    radiators => radiator_map_t( error )
+    ASSERT( error%is_success() )
+    tuvx => tuvx_t( config_path, grids, profiles, radiators, error )
+    ASSERT( error%is_success() )
+    deallocate( heights )
+    deallocate( wavelengths )
+    deallocate( temperatures )
+    deallocate( grids )
+    deallocate( profiles )
+    deallocate( radiators )
+    grids => tuvx%get_grids( error )
+    ASSERT( error%is_success() )
+    profiles => tuvx%get_profiles( error )
+    ASSERT( error%is_success() )
+    radiators => tuvx%get_radiators( error )
+    ASSERT( error%is_success() )
+    heights => grids%get( "height", "km", error )
+    ASSERT( error%is_success() )
+    wavelengths => grids%get( "wavelength", "nm", error )
+    ASSERT( error%is_success() )
+    temperatures => profiles%get( "temperature", "K", error )
+    ASSERT( error%is_success() )
+    call temperatures%set_edge_values( [ 300.0_dk, 275.0_dk, 260.0_dk, 255.0_dk ], error )
+    ASSERT( error%is_success() )
+    call temperatures%set_midpoint_values( [ 287.5_dk, 267.5_dk, 257.5_dk ], error )
+    ASSERT( error%is_success() )
+    call tuvx%run( 0.1_dk, 1.1_dk, photo_rate_consts, heating_rates, error )
+    ASSERT( error%is_success() )
+    do i = 1, number_of_reactions
+      do j = 1, number_of_layers + 1
+        a = photo_rate_consts(j,i)
+        b = expected_photolysis_rate_constants(j,i)
+        ASSERT_NEAR( a, b, 1.0e-5_dk )
+        a = heating_rates(j,i)
+        b = expected_heating_rates(j,i)
+        ASSERT_NEAR( a, b, 1.0e-5_dk )
+      end do
+    end do
+    allocate( temp_vals(4) )
+    call heights%get_edges( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 0.0_dk )
+    ASSERT_EQ( temp_vals(2), 1.0_dk )
+    ASSERT_EQ( temp_vals(3), 2.0_dk )
+    ASSERT_EQ( temp_vals(4), 3.0_dk )
+    deallocate( temp_vals )
+    allocate( temp_vals(3) )
+    call heights%get_midpoints( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 0.5_dk )
+    ASSERT_EQ( temp_vals(2), 1.5_dk )
+    ASSERT_EQ( temp_vals(3), 2.5_dk )
+    deallocate( temp_vals )
+    allocate( temp_vals(6) )
+    call wavelengths%get_edges( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 300.0_dk )
+    ASSERT_EQ( temp_vals(2), 400.0_dk )
+    ASSERT_EQ( temp_vals(3), 500.0_dk )
+    ASSERT_EQ( temp_vals(4), 600.0_dk )
+    ASSERT_EQ( temp_vals(5), 700.0_dk )
+    ASSERT_EQ( temp_vals(6), 800.0_dk )
+    deallocate( temp_vals )
+    allocate( temp_vals(5) )
+    call wavelengths%get_midpoints( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 350.0_dk )
+    ASSERT_EQ( temp_vals(2), 450.0_dk )
+    ASSERT_EQ( temp_vals(3), 550.0_dk )
+    ASSERT_EQ( temp_vals(4), 650.0_dk )
+    ASSERT_EQ( temp_vals(5), 750.0_dk )
+    deallocate( temp_vals )
+    allocate( temp_vals(4) )
+    call temperatures%get_edge_values( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 300.0_dk )
+    ASSERT_EQ( temp_vals(2), 275.0_dk )
+    ASSERT_EQ( temp_vals(3), 260.0_dk )
+    ASSERT_EQ( temp_vals(4), 255.0_dk )
+    deallocate( temp_vals )
+    allocate( temp_vals(3) )
+    call temperatures%get_midpoint_values( temp_vals, error )
+    ASSERT( error%is_success() )
+    ASSERT_EQ( temp_vals(1), 287.5_dk )
+    ASSERT_EQ( temp_vals(2), 267.5_dk )
+    ASSERT_EQ( temp_vals(3), 257.5_dk )
+    deallocate( temp_vals )
+    deallocate( tuvx )
+    deallocate( heights )
+    deallocate( wavelengths )
+    deallocate( temperatures )
+    deallocate( radiators )
+    deallocate( profiles )
+    deallocate( grids ) 
+
+  end subroutine test_tuvx_data_from_host
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 end program test_tuvx_connection

fortran/tuvx/tuvx.F90

@@ -63,6 +63,18 @@ function get_radiator_map_c(tuvx, error) bind(C, name="GetRadiatorMap")
       type(error_t_c), intent(inout) :: error
       type(c_ptr)                    :: get_radiator_map_c
     end function get_radiator_map_c
+
+    subroutine run_tuvx_c(tuvx, solar_zenith_angle, earth_sun_distance, &
+        photolysis_rate_constants, heating_rates, error) bind(C, name="RunTuvx")
+      use musica_util, only: error_t_c
+      use iso_c_binding, only: c_ptr, c_double
+      type(c_ptr), value,         intent(in)    :: tuvx
+      real(kind=c_double), value, intent(in)    :: solar_zenith_angle
+      real(kind=c_double), value, intent(in)    :: earth_sun_distance
+      type(c_ptr), value,         intent(in)    :: photolysis_rate_constants
+      type(c_ptr), value,         intent(in)    :: heating_rates
+      type(error_t_c),            intent(inout) :: error
+    end subroutine run_tuvx_c
   end interface
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -76,6 +88,8 @@ end function get_radiator_map_c
     procedure :: get_profiles
     ! Create a radiator map
     procedure :: get_radiators
+    ! Run the calculator
+    procedure :: run
     ! Deallocate the tuvx instance
     final :: finalize
   end type tuvx_t
@@ -196,6 +210,36 @@ function get_radiators(this, error) result(radiator_map)
 
   end function get_radiators
 
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Run the calculator
+  subroutine run(this, solar_zenith_angle, earth_sun_distance, &
+      photolysis_rate_constants, heating_rates, error)
+    use iso_c_binding, only: c_double, c_ptr, c_loc
+    use musica_util, only: error_t, error_t_c, dk => musica_dk
+
+    ! Arguments
+    class(tuvx_t),         intent(inout) :: this
+    real(kind=dk),         intent(in)    :: solar_zenith_angle             ! radians
+    real(kind=dk),         intent(in)    :: earth_sun_distance             ! AU
+    real(kind=dk), target, intent(inout) :: photolysis_rate_constants(:,:) ! s-1 (layer, reaction)
+    real(kind=dk), target, intent(inout) :: heating_rates(:,:)             ! K s-1 (layer, reaction)
+    type(error_t),         intent(inout) :: error
+
+    ! Local variables
+    type(error_t_c) :: error_c
+    type(c_ptr)     :: photo_rate_c, heating_c
+
+    photo_rate_c = c_loc(photolysis_rate_constants(1,1))
+    heating_c    = c_loc(heating_rates(1,1))
+    call run_tuvx_c(this%ptr_, &
+                    real(solar_zenith_angle, kind=c_double), &
+                    real(earth_sun_distance, kind=c_double), &
+                    photo_rate_c, heating_c, error_c)
+    error = error_t(error_c)
+
+  end subroutine run
+
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !> Deallocate the tuvx instance

include/musica/tuvx/grid_map.hpp

@@ -44,6 +44,8 @@ namespace musica
    private:
     void *grid_map_;
     bool owns_grid_map_;
+
+    friend class TUVX;
   };
 
 #ifdef __cplusplus

include/musica/tuvx/profile_map.hpp

@@ -45,6 +45,8 @@ namespace musica
    private:
     void *profile_map_;
     bool owns_profile_map_;
+
+    friend class TUVX;
   };
 
 #ifdef __cplusplus

include/musica/tuvx/radiator_map.hpp

@@ -43,6 +43,8 @@ namespace musica
    private:
     void *radiator_map_;
     bool owns_radiator_map_;
+
+    friend class TUVX;
   };
 
 #ifdef __cplusplus

include/musica/tuvx/tuvx.hpp

@@ -48,10 +48,20 @@ namespace musica
     /// @return a radiator map pointer
     RadiatorMap *CreateRadiatorMap(Error *error);
 
+    /// @brief Run the TUV-x photolysis calculator
+    /// @param solar_zenith_angle Solar zenith angle [radians]
+    /// @param earth_sun_distance Earth-Sun distance [AU]
+    /// @param photolysis_rate_constants Photolysis rate constant for each layer and reaction [s^-1]
+    /// @param heating_rates Heating rates for each layer and reaction [K/s]
+    /// @param error Error struct to indicate success or failure
+    void Run(const double solar_zenith_angle, const double earth_sun_distance,
+             double * const photolysis_rate_constants, double * const heating_rates, Error * const error);
+
     ~TUVX();
 
    private:
     void *tuvx_;
+    int number_of_layers_;
   };
 
 #ifdef __cplusplus
@@ -66,15 +76,19 @@ namespace musica
     GridMap *GetGridMap(TUVX *tuvx, Error *error);
     ProfileMap *GetProfileMap(TUVX *tuvx, Error *error);
     RadiatorMap *GetRadiatorMap(TUVX *tuvx, Error *error);
+    void RunTuvx(TUVX *tuvx, const double solar_zenith_angle, const double earth_sun_distance,
+                 double * const photolysis_rate_constants, double * const heating_rates, Error * const error);
 
     // for use by musica interanlly. If tuvx ever gets rewritten in C++, these functions will
     // go away but the C API will remain the same and downstream projects (like CAM-SIMA) will
     // not need to change
-    void *InternalCreateTuvx(const char *config_path, std::size_t config_path_length, int *error_code);
+    void *InternalCreateTuvx(const char *config_path, std::size_t config_path_length, void *grid_map, void *profile_map, void *radiator_map, int *number_of_layers, int *error_code);
     void InternalDeleteTuvx(void *tuvx, int *error_code);
     void *InternalGetGridMap(void *tuvx, int *error_code);
     void *InternalGetProfileMap(void *tuvx, int *error_code);
     void *InternalGetRadiatorMap(void *tuvx, int *error_code);
+    void InternalRunTuvx(void *tuvx, const int number_of_layers, const double solar_zenith_angle, const double earth_sun_distance,
+                         double *photolysis_rate_constants, double *heating_rates, int *error_code);
 
 #ifdef __cplusplus
   }