Skip to content

Commit

Permalink
TUV-x photolysis rate constant and heating rate ordering access (#235)
Browse files Browse the repository at this point in the history
* return mappings for photo and heating rates from TUV-x

* return TUV-x rate ordering from Fortran wrapper

* address reviewr comments
  • Loading branch information
mattldawson authored Nov 1, 2024
1 parent dbbdb22 commit ff9f8de
Show file tree
Hide file tree
Showing 17 changed files with 463 additions and 24 deletions.
82 changes: 75 additions & 7 deletions fortran/test/unit/tuvx.F90
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,13 @@ program test_tuvx_connection

integer, parameter :: number_of_layers = 3
integer, parameter :: number_of_wavelengths = 5
integer, parameter :: number_of_reactions = 2
integer, parameter :: number_of_reactions = 3
integer, parameter :: number_of_heating_rates = 2

real(dk), parameter :: expected_photolysis_rate_constants(4,2) = reshape( [ &
real(dk), parameter :: expected_photolysis_rate_constants(4,3) = 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] )
2.49575956372508e-27_dk, 4.58686176250519e-20_dk, 2.22679622672858e-13_dk, 2.29392676897831e-07_dk, &
1.78278752667774e-27_dk, 3.28516384208994e-20_dk, 1.69678305465474e-13_dk, 1.97484189784941e-07_dk ], [4,3] )
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] )
Expand Down Expand Up @@ -65,7 +67,7 @@ end subroutine test_tuvx
!> Test TUV-x with a fixed configuration
subroutine test_tuvx_fixed( )

use musica_util, only : assert, error_t
use musica_util, only : assert, error_t, mappings_t
use musica_tuvx_grid, only : grid_t
use musica_tuvx_grid_map, only : grid_map_t
use musica_tuvx_profile, only : profile_t
Expand All @@ -79,10 +81,11 @@ subroutine test_tuvx_fixed( )
type(profile_map_t), pointer :: profiles
type(radiator_map_t), pointer :: radiators
type(tuvx_t), pointer :: tuvx
type(mappings_t), pointer :: photo_mappings, heating_mappings
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)
real(dk) :: heating_rates(number_of_layers + 1, number_of_heating_rates)
integer :: i, j
real(dk) :: a, b

Expand All @@ -101,11 +104,43 @@ subroutine test_tuvx_fixed( )
a = photo_rate_consts(j,i)
b = expected_photolysis_rate_constants(j,i)
ASSERT_NEAR( a, b, 1.0e-5_dk )
end do
end do
do i = 1, number_of_heating_rates
do j = 1, number_of_layers + 1
a = heating_rates(j,i)
b = expected_heating_rates(j,i)
ASSERT_NEAR( a, b, 1.0e-5_dk )
end do
end do
photo_mappings => tuvx%get_photolysis_rate_constants_ordering( error )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%size(), number_of_reactions )
ASSERT_EQ( photo_mappings%name( 1 ), "jfoo" )
ASSERT_EQ( photo_mappings%index( 1 ), 1 )
ASSERT_EQ( photo_mappings%index( "jfoo", error ), 1 )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%name( 2 ), "jbar" )
ASSERT_EQ( photo_mappings%index( 2 ), 2 )
ASSERT_EQ( photo_mappings%index( "jbar", error ), 2 )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%name( 3 ), "jbaz" )
ASSERT_EQ( photo_mappings%index( 3 ), 3 )
ASSERT_EQ( photo_mappings%index( "jbaz", error ), 3 )
ASSERT( error%is_success() )
heating_mappings => tuvx%get_heating_rates_ordering( error )
ASSERT( error%is_success() )
ASSERT_EQ( heating_mappings%size(), number_of_heating_rates )
ASSERT_EQ( heating_mappings%name( 1 ), "jfoo" )
ASSERT_EQ( heating_mappings%index( 1 ), 1 )
ASSERT_EQ( heating_mappings%index( "jfoo", error ), 1 )
ASSERT( error%is_success() )
ASSERT_EQ( heating_mappings%name( 2 ), "jbar" )
ASSERT_EQ( heating_mappings%index( 2 ), 2 )
ASSERT_EQ( heating_mappings%index( "jbar", error ), 2 )
ASSERT( error%is_success() )
deallocate( photo_mappings )
deallocate( heating_mappings )
deallocate( tuvx )
deallocate( radiators )
deallocate( profiles )
Expand All @@ -118,7 +153,7 @@ 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_util, only : assert, error_t, mappings_t
use musica_tuvx_grid, only : grid_t
use musica_tuvx_grid_map, only : grid_map_t
use musica_tuvx_profile, only : profile_t
Expand All @@ -134,10 +169,11 @@ subroutine test_tuvx_data_from_host( )
type(profile_map_t), pointer :: profiles
type(radiator_map_t), pointer :: radiators
type(tuvx_t), pointer :: tuvx
type(mappings_t), pointer :: photo_mappings, heating_mappings
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)
real(dk) :: heating_rates(number_of_layers + 1, number_of_heating_rates)
integer :: i, j
real(dk) :: a, b
real(dk), allocatable :: temp_vals(:)
Expand Down Expand Up @@ -219,6 +255,10 @@ subroutine test_tuvx_data_from_host( )
a = photo_rate_consts(j,i)
b = expected_photolysis_rate_constants(j,i)
ASSERT_NEAR( a, b, 1.0e-5_dk )
end do
end do
do i = 1, number_of_heating_rates
do j = 1, number_of_layers + 1
a = heating_rates(j,i)
b = expected_heating_rates(j,i)
ASSERT_NEAR( a, b, 1.0e-5_dk )
Expand Down Expand Up @@ -273,6 +313,34 @@ subroutine test_tuvx_data_from_host( )
ASSERT_EQ( temp_vals(1), 287.5_dk )
ASSERT_EQ( temp_vals(2), 267.5_dk )
ASSERT_EQ( temp_vals(3), 257.5_dk )
photo_mappings => tuvx%get_photolysis_rate_constants_ordering( error )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%size(), number_of_reactions )
ASSERT_EQ( photo_mappings%name( 1 ), "jfoo" )
ASSERT_EQ( photo_mappings%index( 1 ), 1 )
ASSERT_EQ( photo_mappings%index( "jfoo", error ), 1 )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%name( 2 ), "jbar" )
ASSERT_EQ( photo_mappings%index( 2 ), 2 )
ASSERT_EQ( photo_mappings%index( "jbar", error ), 2 )
ASSERT( error%is_success() )
ASSERT_EQ( photo_mappings%name( 3 ), "jbaz" )
ASSERT_EQ( photo_mappings%index( 3 ), 3 )
ASSERT_EQ( photo_mappings%index( "jbaz", error ), 3 )
ASSERT( error%is_success() )
heating_mappings => tuvx%get_heating_rates_ordering( error )
ASSERT( error%is_success() )
ASSERT_EQ( heating_mappings%size(), number_of_heating_rates )
ASSERT_EQ( heating_mappings%name( 1 ), "jfoo" )
ASSERT_EQ( heating_mappings%index( 1 ), 1 )
ASSERT_EQ( heating_mappings%index( "jfoo", error ), 1 )
ASSERT( error%is_success() )
ASSERT_EQ( heating_mappings%name( 2 ), "jbar" )
ASSERT_EQ( heating_mappings%index( 2 ), 2 )
ASSERT_EQ( heating_mappings%index( "jbar", error ), 2 )
ASSERT( error%is_success() )
deallocate( photo_mappings )
deallocate( heating_mappings )
deallocate( temp_vals )
deallocate( tuvx )
deallocate( heights )
Expand Down
66 changes: 66 additions & 0 deletions fortran/tuvx/tuvx.F90
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,24 @@ function get_radiator_map_c(tuvx, error) bind(C, name="GetRadiatorMap")
type(c_ptr) :: get_radiator_map_c
end function get_radiator_map_c

function get_photolysis_rate_constants_ordering_c(tuvx, error) &
bind(C, name="GetPhotolysisRateConstantsOrdering")
use musica_util, only: error_t_c, mappings_t_c
use iso_c_binding, only: c_ptr
type(c_ptr), value, intent(in) :: tuvx
type(error_t_c), intent(inout) :: error
type(mappings_t_c) :: get_photolysis_rate_constants_ordering_c
end function get_photolysis_rate_constants_ordering_c

function get_heating_rates_ordering_c(tuvx, error) &
bind(C, name="GetHeatingRatesOrdering")
use musica_util, only: error_t_c, mappings_t_c
use iso_c_binding, only: c_ptr
type(c_ptr), value, intent(in) :: tuvx
type(error_t_c), intent(inout) :: error
type(mappings_t_c) :: get_heating_rates_ordering_c
end function get_heating_rates_ordering_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
Expand All @@ -88,6 +106,10 @@ end subroutine run_tuvx_c
procedure :: get_profiles
! Create a radiator map
procedure :: get_radiators
! Get photolysis rate constant ordering
procedure :: get_photolysis_rate_constants_ordering
! Get heating rate ordering
procedure :: get_heating_rates_ordering
! Run the calculator
procedure :: run
! Deallocate the tuvx instance
Expand Down Expand Up @@ -212,6 +234,50 @@ end function get_radiators

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> Get the photolysis rate constant ordering
function get_photolysis_rate_constants_ordering(this, error) &
result(mappings)
use musica_util, only: error_t, error_t_c, mappings_t

! Arguments
class(tuvx_t), intent(inout) :: this
type(error_t), intent(inout) :: error

! Return value
type(mappings_t), pointer :: mappings

! Local variables
type(error_t_c) :: error_c

mappings => &
mappings_t(get_photolysis_rate_constants_ordering_c(this%ptr_, error_c))
error = error_t(error_c)

end function get_photolysis_rate_constants_ordering

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> Get the heating rate ordering
function get_heating_rates_ordering(this, error) result(mappings)
use musica_util, only: error_t, error_t_c, mappings_t

! Arguments
class(tuvx_t), intent(inout) :: this
type(error_t), intent(inout) :: error

! Return value
type(mappings_t), pointer :: mappings

! Local variables
type(error_t_c) :: error_c

mappings => mappings_t(get_heating_rates_ordering_c(this%ptr_, error_c))
error = error_t(error_c)

end function get_heating_rates_ordering

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> Run the calculator
subroutine run(this, solar_zenith_angle, earth_sun_distance, &
photolysis_rate_constants, heating_rates, error)
Expand Down
60 changes: 59 additions & 1 deletion include/musica/tuvx/tuvx.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,16 @@ namespace musica
/// @return a radiator map pointer
RadiatorMap *CreateRadiatorMap(Error *error);

/// @brief Returns the ordering of photolysis rate constants
/// @param error Error struct to indicate success or failure
/// @return Array of photolysis rate constant name-index pairs
Mappings GetPhotolysisRateConstantsOrdering(Error *error);

/// @brief Returns the ordering of heating rates
/// @param error Error struct to indicate success or failure
/// @return Array of heating rate name-index pairs
Mappings GetHeatingRatesOrdering(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]
Expand Down Expand Up @@ -75,11 +85,57 @@ namespace musica

// The external C API for TUVX
// callable by wrappers in other languages

/// @brief Creates a TUVX instance by passing a configuration file path and host-defined grids, profiles, and radiators
/// @param config_path Path to configuration file
/// @param grids Grid map from host application
/// @param profiles Profile map from host application
/// @param radiators Radiator map from host application
/// @param error Error struct to indicate success or failure
TUVX *CreateTuvx(const char *config_path, GridMap *grids, ProfileMap *profiles, RadiatorMap *radiators, Error *error);

/// @brief Deletes a TUVX instance
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
void DeleteTuvx(const TUVX *tuvx, Error *error);

/// @brief Returns the set of grids used by TUVX
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
/// @return Grid map
GridMap *GetGridMap(TUVX *tuvx, Error *error);

/// @brief Returns the set of profiles used by TUVX
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
/// @return Profile map
ProfileMap *GetProfileMap(TUVX *tuvx, Error *error);

/// @brief Returns the set of radiators used by TUVX
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
/// @return Radiator map
RadiatorMap *GetRadiatorMap(TUVX *tuvx, Error *error);

/// @brief Returns the ordering photolysis rate constants
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
/// @return Array of photolysis rate constant name-index pairs
Mappings GetPhotolysisRateConstantsOrdering(TUVX *tuvx, Error *error);

/// @brief Returns the ordering of heating rates
/// @param tuvx Pointer to TUVX instance
/// @param error Error struct to indicate success or failure
/// @return Array of heating rate name-index pairs
Mappings GetHeatingRatesOrdering(TUVX *tuvx, Error *error);

/// @brief Run the TUV-x photolysis calculator
/// @param tuvx Pointer to TUVX instance
/// @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 RunTuvx(
TUVX *tuvx,
const double solar_zenith_angle,
Expand All @@ -88,7 +144,7 @@ namespace musica
double *const heating_rates,
Error *const error);

// for use by musica interanlly. If tuvx ever gets rewritten in C++, these functions will
// for use by musica internally. 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(
Expand All @@ -103,6 +159,8 @@ namespace musica
void *InternalGetGridMap(void *tuvx, int *error_code);
void *InternalGetProfileMap(void *tuvx, int *error_code);
void *InternalGetRadiatorMap(void *tuvx, int *error_code);
Mappings InternalGetPhotolysisRateConstantsOrdering(void *tuvx, int *error_code);
Mappings InternalGetHeatingRatesOrdering(void *tuvx, int *error_code);
void InternalRunTuvx(
void *tuvx,
const int number_of_layers,
Expand Down
5 changes: 5 additions & 0 deletions include/musica/util.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,11 @@ namespace musica
/// @return The Configuration
Configuration LoadConfigurationFromFile(const char* filename, Error* error);

/// @brief Allocates an array of Mappings
/// @param size The size of the array
/// @return The array of Mappings
Mapping* AllocateMappingArray(const std::size_t size);

/// @brief Allocate a new Mappings struct
/// @param size The size of the Mappings
/// @return The Mappings
Expand Down
15 changes: 15 additions & 0 deletions src/test/data/tuvx/fixed/config.json
Original file line number Diff line number Diff line change
Expand Up @@ -143,6 +143,21 @@
"heating": {
"energy term": 550.0
}
},
{
"name": "jbaz",
"cross section": {
"type": "base",
"netcdf files": [
{
"file path": "test/data/tuvx/fixed/foo_cross_section.nc"
}
]
},
"quantum yield": {
"type": "base",
"constant value": 1.0
}
}
]
}
Expand Down
15 changes: 15 additions & 0 deletions src/test/data/tuvx/from_host/config.json
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,21 @@
"heating": {
"energy term": 550.0
}
},
{
"name": "jbaz",
"cross section": {
"type": "base",
"netcdf files": [
{
"file path": "test/data/tuvx/from_host/foo_cross_section.nc"
}
]
},
"quantum yield": {
"type": "base",
"constant value": 1.0
}
}
]
}
Expand Down
Loading

0 comments on commit ff9f8de

Please sign in to comment.