From 1639f0b54946c2dfc971961ed822dbcc38585b00 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 3 Oct 2024 09:15:56 -0600 Subject: [PATCH 01/13] Auto-format code using Clang-Format (#231) Co-authored-by: GitHub Actions --- include/musica/micm.hpp | 2 +- src/micm/micm.cpp | 15 +++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/include/musica/micm.hpp b/include/musica/micm.hpp index 6d754984..bb97cc7a 100644 --- a/include/musica/micm.hpp +++ b/include/musica/micm.hpp @@ -17,12 +17,12 @@ #include #include -#include #include #include #include #include #include +#include #include #ifndef MICM_VECTOR_MATRIX_SIZE diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index ceefe22f..ff037d7f 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -328,7 +328,8 @@ namespace musica solver_config.ReadAndParse(std::filesystem::path(config_path)); solver_parameters_ = std::make_unique(solver_config.GetSolverParams()); - auto solver = std::make_unique(micm::CpuSolverBuilder( + auto solver = + std::make_unique(micm::CpuSolverBuilder( micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) .SetSystem(solver_parameters_->system_) .SetReactions(solver_parameters_->processes_) @@ -336,7 +337,8 @@ namespace musica .SetIgnoreUnusedSpecies(true) .Build()); auto state = solver->GetState(); - rosenbrock_standard_ = std::pair, StandardState>(std::move(solver), std::move(state)); + rosenbrock_standard_ = + std::pair, StandardState>(std::move(solver), std::move(state)); DeleteError(error); *error = NoError(); @@ -401,8 +403,9 @@ namespace musica .SetIgnoreUnusedSpecies(true) .Build()); auto state = solver->GetState(); - - backward_euler_standard_ = std::pair, StandardState>(std::move(solver), std::move(state)); + + backward_euler_standard_ = + std::pair, StandardState>(std::move(solver), std::move(state)); *error = NoError(); } @@ -426,8 +429,8 @@ namespace musica { try { - auto& solver = solver_state_pair.first; - auto& state = solver_state_pair.second; + auto &solver = solver_state_pair.first; + auto &state = solver_state_pair.second; const std::size_t num_species = state.variables_.NumColumns(); const std::size_t num_custom_rate_parameters = state.custom_rate_parameters_.NumColumns(); From dbbdb22f5f2807e27c2695db85291951ba178634 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Wed, 9 Oct 2024 14:19:43 -0600 Subject: [PATCH 02/13] fix bug in getting number of grid sections (#233) * fix bug in getting number of grid sections * address reviewer comments --- fortran/test/unit/tuvx.F90 | 26 ++++++++++++++++++++++++-- src/tuvx/grid.cpp | 2 +- src/tuvx/interface_grid.F90 | 12 ++++++------ 3 files changed, 31 insertions(+), 9 deletions(-) diff --git a/fortran/test/unit/tuvx.F90 b/fortran/test/unit/tuvx.F90 index ce477652..b7238889 100644 --- a/fortran/test/unit/tuvx.F90 +++ b/fortran/test/unit/tuvx.F90 @@ -144,9 +144,9 @@ subroutine test_tuvx_data_from_host( ) 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 ) + call heights%set_edges( [ 0.0_dk, 1.2_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 ) + call heights%set_midpoints( [ 0.5_dk, 1.3_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 ) @@ -157,6 +157,27 @@ subroutine test_tuvx_data_from_host( ) ASSERT( error%is_success() ) call grids%add( heights, error ) ASSERT( error%is_success() ) + ASSERT_EQ( heights%number_of_sections( error ), 3 ) + ASSERT( error%is_success() ) + 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.2_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.3_dk ) + ASSERT_EQ( temp_vals(3), 2.5_dk ) + deallocate( temp_vals ) + 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 ) + ASSERT( error%is_success() ) call grids%add( wavelengths, error ) ASSERT( error%is_success() ) temperatures => profile_t( "temperature", "K", heights, error ) @@ -204,6 +225,7 @@ subroutine test_tuvx_data_from_host( ) end do end do allocate( temp_vals(4) ) + ASSERT_EQ( heights%number_of_sections( error ), 3 ) call heights%get_edges( temp_vals, error ) ASSERT( error%is_success() ) ASSERT_EQ( temp_vals(1), 0.0_dk ) diff --git a/src/tuvx/grid.cpp b/src/tuvx/grid.cpp index 071b7b0e..1e39baa6 100644 --- a/src/tuvx/grid.cpp +++ b/src/tuvx/grid.cpp @@ -96,7 +96,7 @@ namespace musica std::size_t Grid::GetNumSections(Error *error) { int error_code = 0; - std::size_t n_sections = InternalGetNumSections(grid_, &error_code); + std::size_t n_sections = InternalGetNumSections(updater_, &error_code); if (error_code != 0) { *error = Error{ 1, CreateString(MUSICA_ERROR_CATEGORY), CreateString("Failed to get number of sections") }; diff --git a/src/tuvx/interface_grid.F90 b/src/tuvx/interface_grid.F90 index 41fb7333..70faad54 100644 --- a/src/tuvx/interface_grid.F90 +++ b/src/tuvx/interface_grid.F90 @@ -117,23 +117,23 @@ end subroutine internal_delete_grid_updater !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - function internal_get_num_sections(grid, error_code) & + function internal_get_num_sections(updater, error_code) & bind(C, name="InternalGetNumSections") result(num_sections) use iso_c_binding, only: c_ptr, c_f_pointer, c_int, c_size_t - use tuvx_grid_from_host, only: grid_from_host_t + use tuvx_grid_from_host, only: grid_updater_t ! arguments - type(c_ptr), value, intent(in) :: grid + type(c_ptr), value, intent(in) :: updater integer(kind=c_int), intent(out) :: error_code ! output integer(kind=c_size_t) :: num_sections ! variables - type(grid_from_host_t), pointer :: f_grid + type(grid_updater_t), pointer :: f_updater - call c_f_pointer(grid, f_grid) - num_sections = f_grid%size( ) + call c_f_pointer(updater, f_updater) + num_sections = f_updater%grid_%size( ) end function internal_get_num_sections From ff9f8dee693cc165223ec2a6f4e0bb2f4c62409b Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 1 Nov 2024 12:57:50 -0700 Subject: [PATCH 03/13] TUV-x photolysis rate constant and heating rate ordering access (#235) * return mappings for photo and heating rates from TUV-x * return TUV-x rate ordering from Fortran wrapper * address reviewr comments --- fortran/test/unit/tuvx.F90 | 82 +++++++++++++++++-- fortran/tuvx/tuvx.F90 | 66 +++++++++++++++ include/musica/tuvx/tuvx.hpp | 60 +++++++++++++- include/musica/util.hpp | 5 ++ src/test/data/tuvx/fixed/config.json | 15 ++++ src/test/data/tuvx/from_host/config.json | 15 ++++ src/test/unit/tuvx/tuvx_run_from_config.cpp | 65 +++++++++++++-- src/tuvx/CMakeLists.txt | 2 +- src/tuvx/interface.F90 | 81 +++++++++++++++++- src/tuvx/interface_grid.F90 | 1 - src/tuvx/interface_grid_map.F90 | 1 - src/tuvx/interface_profile.F90 | 1 - src/tuvx/interface_profile_map.F90 | 1 - src/tuvx/interface_radiator_map.F90 | 1 - .../{tuvx_util.F90 => interface_util.F90} | 50 ++++++++++- src/tuvx/tuvx.cpp | 36 ++++++++ src/util.cpp | 5 ++ 17 files changed, 463 insertions(+), 24 deletions(-) rename src/tuvx/{tuvx_util.F90 => interface_util.F90} (86%) diff --git a/fortran/test/unit/tuvx.F90 b/fortran/test/unit/tuvx.F90 index b7238889..2192ff3b 100644 --- a/fortran/test/unit/tuvx.F90 +++ b/fortran/test/unit/tuvx.F90 @@ -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] ) @@ -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 @@ -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 @@ -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 ) @@ -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 @@ -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(:) @@ -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 ) @@ -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 ) diff --git a/fortran/tuvx/tuvx.F90 b/fortran/tuvx/tuvx.F90 index fe9d970a..d5c2ca2f 100644 --- a/fortran/tuvx/tuvx.F90 +++ b/fortran/tuvx/tuvx.F90 @@ -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 @@ -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 @@ -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) diff --git a/include/musica/tuvx/tuvx.hpp b/include/musica/tuvx/tuvx.hpp index 7fe798b2..fad4c10b 100644 --- a/include/musica/tuvx/tuvx.hpp +++ b/include/musica/tuvx/tuvx.hpp @@ -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] @@ -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, @@ -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( @@ -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, diff --git a/include/musica/util.hpp b/include/musica/util.hpp index a791e8a2..ce968479 100644 --- a/include/musica/util.hpp +++ b/include/musica/util.hpp @@ -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 diff --git a/src/test/data/tuvx/fixed/config.json b/src/test/data/tuvx/fixed/config.json index be54d22b..8ecbc9e2 100644 --- a/src/test/data/tuvx/fixed/config.json +++ b/src/test/data/tuvx/fixed/config.json @@ -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 + } } ] } diff --git a/src/test/data/tuvx/from_host/config.json b/src/test/data/tuvx/from_host/config.json index 4fc9674f..2191eec1 100644 --- a/src/test/data/tuvx/from_host/config.json +++ b/src/test/data/tuvx/from_host/config.json @@ -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 + } } ] } diff --git a/src/test/unit/tuvx/tuvx_run_from_config.cpp b/src/test/unit/tuvx/tuvx_run_from_config.cpp index b01f4405..8cb440fd 100644 --- a/src/test/unit/tuvx/tuvx_run_from_config.cpp +++ b/src/test/unit/tuvx/tuvx_run_from_config.cpp @@ -6,9 +6,10 @@ using namespace musica; // Expected values for photolysis rate constants and heating rates // were determined by running the stand-alone TUV-x model with the fixed configuration. -const double expected_photolysis_rate_constants[2][4] = { +const double expected_photolysis_rate_constants[3][4] = { { 8.91393763338872e-28, 1.64258192104497e-20, 8.48391527327371e-14, 9.87420948924703e-08 }, - { 2.49575956372508e-27, 4.58686176250519e-20, 2.22679622672858e-13, 2.29392676897831e-07 } + { 2.49575956372508e-27, 4.58686176250519e-20, 2.22679622672858e-13, 2.29392676897831e-07 }, + { 1.78278752667774e-27, 3.28516384208994e-20, 1.69678305465474e-13, 1.97484189784941e-07} }; const double expected_heating_rates[2][4] = { { 1.12394047546984e-46, 2.04518267143613e-39, 7.44349752571804e-33, 5.42628100199216e-28 }, @@ -30,6 +31,7 @@ class TuvxRunTest : public ::testing::Test int number_of_layers; int number_of_wavelengths; int number_of_reactions; + int number_of_heating_rates; double* photolysis_rate_constants; double* heating_rates; @@ -46,6 +48,7 @@ class TuvxRunTest : public ::testing::Test number_of_layers = 0; number_of_wavelengths = 0; number_of_reactions = 0; + number_of_heating_rates = 0; photolysis_rate_constants = nullptr; heating_rates = nullptr; } @@ -73,9 +76,10 @@ class TuvxRunTest : public ::testing::Test ASSERT_TRUE(IsSuccess(error)); number_of_layers = 3; number_of_wavelengths = 5; - number_of_reactions = 2; + number_of_reactions = 3; + number_of_heating_rates = 2; photolysis_rate_constants = new double[(number_of_layers + 1) * number_of_reactions]; - heating_rates = new double[(number_of_layers + 1) * number_of_reactions]; + heating_rates = new double[(number_of_layers + 1) * number_of_heating_rates]; DeleteError(&error); } @@ -99,9 +103,10 @@ class TuvxRunTest : public ::testing::Test ASSERT_TRUE(IsSuccess(error)); number_of_layers = 3; number_of_wavelengths = 5; - number_of_reactions = 2; + number_of_reactions = 3; + number_of_heating_rates = 2; photolysis_rate_constants = new double[(number_of_layers + 1) * number_of_reactions]; - heating_rates = new double[(number_of_layers + 1) * number_of_reactions]; + heating_rates = new double[(number_of_layers + 1) * number_of_heating_rates]; DeleteError(&error); } @@ -148,12 +153,36 @@ TEST_F(TuvxRunTest, CreateTuvxInstanceWithJsonConfig) photolysis_rate_constants[i * (number_of_layers + 1) + j], expected_photolysis_rate_constants[i][j], expected_photolysis_rate_constants[i][j] * 1.0e-5); + } + } + for (int i = 0; i < number_of_heating_rates; i++) + { + for (int j = 0; j < number_of_layers + 1; j++) + { EXPECT_NEAR( heating_rates[i * (number_of_layers + 1) + j], expected_heating_rates[i][j], expected_heating_rates[i][j] * 1.0e-5); } } + Mappings photo_rate_labels = GetPhotolysisRateConstantsOrdering(tuvx, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(photo_rate_labels.size_, 3); + ASSERT_STREQ(photo_rate_labels.mappings_[0].name_.value_, "jfoo"); + ASSERT_EQ(photo_rate_labels.mappings_[0].index_, 0); + ASSERT_STREQ(photo_rate_labels.mappings_[1].name_.value_, "jbar"); + ASSERT_EQ(photo_rate_labels.mappings_[1].index_, 1); + ASSERT_STREQ(photo_rate_labels.mappings_[2].name_.value_, "jbaz"); + ASSERT_EQ(photo_rate_labels.mappings_[2].index_, 2); + Mappings heating_rate_labels = GetHeatingRatesOrdering(tuvx, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(heating_rate_labels.size_, 2); + ASSERT_STREQ(heating_rate_labels.mappings_[0].name_.value_, "jfoo"); + ASSERT_EQ(heating_rate_labels.mappings_[0].index_, 0); + ASSERT_STREQ(heating_rate_labels.mappings_[1].name_.value_, "jbar"); + ASSERT_EQ(heating_rate_labels.mappings_[1].index_, 1); + DeleteMappings(&photo_rate_labels); + DeleteMappings(&heating_rate_labels); DeleteError(&error); } @@ -225,12 +254,36 @@ TEST_F(TuvxRunTest, CreateTuvxInstanceWithJsonConfigAndHostData) photolysis_rate_constants[i * (number_of_layers + 1) + j], expected_photolysis_rate_constants[i][j], expected_photolysis_rate_constants[i][j] * 1.0e-5); + } + } + for (int i = 0; i < number_of_heating_rates; i++) + { + for (int j = 0; j < number_of_layers + 1; j++) + { EXPECT_NEAR( heating_rates[i * (number_of_layers + 1) + j], expected_heating_rates[i][j], expected_heating_rates[i][j] * 1.0e-5); } } + Mappings photo_rate_labels = GetPhotolysisRateConstantsOrdering(tuvx, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(photo_rate_labels.size_, 3); + ASSERT_STREQ(photo_rate_labels.mappings_[0].name_.value_, "jfoo"); + ASSERT_EQ(photo_rate_labels.mappings_[0].index_, 0); + ASSERT_STREQ(photo_rate_labels.mappings_[1].name_.value_, "jbar"); + ASSERT_EQ(photo_rate_labels.mappings_[1].index_, 1); + ASSERT_STREQ(photo_rate_labels.mappings_[2].name_.value_, "jbaz"); + ASSERT_EQ(photo_rate_labels.mappings_[2].index_, 2); + Mappings heating_rate_labels = GetHeatingRatesOrdering(tuvx, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(heating_rate_labels.size_, 2); + ASSERT_STREQ(heating_rate_labels.mappings_[0].name_.value_, "jfoo"); + ASSERT_EQ(heating_rate_labels.mappings_[0].index_, 0); + ASSERT_STREQ(heating_rate_labels.mappings_[1].name_.value_, "jbar"); + ASSERT_EQ(heating_rate_labels.mappings_[1].index_, 1); + DeleteMappings(&photo_rate_labels); + DeleteMappings(&heating_rate_labels); GetGridEdges(heights, height_edges, 4, &error); ASSERT_TRUE(IsSuccess(error)); ASSERT_EQ(height_edges[0], 0.0); diff --git a/src/tuvx/CMakeLists.txt b/src/tuvx/CMakeLists.txt index 46fa14bc..bee91711 100644 --- a/src/tuvx/CMakeLists.txt +++ b/src/tuvx/CMakeLists.txt @@ -7,6 +7,7 @@ target_sources(musica interface_profile_map.F90 interface_radiator.F90 interface_radiator_map.F90 + interface_util.F90 grid.cpp grid_map.cpp profile.cpp @@ -14,5 +15,4 @@ target_sources(musica radiator.cpp radiator_map.cpp tuvx.cpp - tuvx_util.F90 ) \ No newline at end of file diff --git a/src/tuvx/interface.F90 b/src/tuvx/interface.F90 index ac393318..66d39b7d 100644 --- a/src/tuvx/interface.F90 +++ b/src/tuvx/interface.F90 @@ -9,7 +9,6 @@ module tuvx_interface use tuvx_grid_warehouse, only : grid_warehouse_t use tuvx_profile_warehouse, only : profile_warehouse_t use tuvx_radiator_warehouse, only : radiator_warehouse_t -use musica_tuvx_util, only : to_f_string, string_t_c use musica_string, only : string_t implicit none @@ -158,6 +157,86 @@ function internal_get_radiator_map(tuvx, error_code) result(radiator_map_ptr) & end function internal_get_radiator_map +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function internal_get_photolysis_rate_constants_ordering(tuvx, error_code) & + result(photolysis_rate_constant_ordering) & + bind(C, name="InternalGetPhotolysisRateConstantsOrdering") + use iso_c_binding, only: c_ptr, c_f_pointer, c_int + use tuvx_interface_util, only: create_string_t_c, mappings_t_c, & + mapping_t_c, allocate_mappings_c + + ! arguments + type(c_ptr), value, intent(in) :: tuvx + integer(kind=c_int), intent(out) :: error_code + + ! result + type(mappings_t_c) :: photolysis_rate_constant_ordering + + ! variables + type(core_t), pointer :: core + type(string_t), allocatable :: labels(:) + type(mapping_t_c), pointer :: mappings(:) + type(c_ptr) :: mappings_ptr + integer :: i, n_labels + + error_code = 0 + call c_f_pointer(tuvx, core) + + labels = core%photolysis_reaction_labels() + n_labels = size(labels) + mappings_ptr = allocate_mappings_c(int(n_labels, kind=c_size_t)) + call c_f_pointer(mappings_ptr, mappings, [ n_labels ]) + do i = 1, n_labels + mappings(i)%name_ = create_string_t_c(labels(i)%val_) + mappings(i)%index_ = int(i-1, kind=c_size_t) + end do + + photolysis_rate_constant_ordering%mappings_ = c_loc(mappings) + photolysis_rate_constant_ordering%size_ = n_labels + + end function internal_get_photolysis_rate_constants_ordering + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function internal_get_heating_rates_ordering(tuvx, error_code) & + result(heating_rates_ordering) & + bind(C, name="InternalGetHeatingRatesOrdering") + use iso_c_binding, only: c_ptr, c_f_pointer, c_int + use tuvx_interface_util, only: create_string_t_c, mappings_t_c, & + mapping_t_c, allocate_mappings_c + + ! arguments + type(c_ptr), value, intent(in) :: tuvx + integer(kind=c_int), intent(out) :: error_code + + ! result + type(mappings_t_c) :: heating_rates_ordering + + ! variables + type(core_t), pointer :: core + type(string_t), allocatable :: labels(:) + type(mapping_t_c), pointer :: mappings(:) + type(c_ptr) :: mappings_ptr + integer :: i, n_labels + + error_code = 0 + call c_f_pointer(tuvx, core) + + labels = core%heating_rate_labels() + n_labels = size(labels) + mappings_ptr = allocate_mappings_c(int(n_labels, kind=c_size_t)) + call c_f_pointer(mappings_ptr, mappings, [ n_labels ]) + do i = 1, n_labels + mappings(i)%name_ = create_string_t_c(labels(i)%val_) + mappings(i)%index_ = int(i-1, kind=c_size_t) + end do + + heating_rates_ordering%mappings_ = c_loc(mappings) + heating_rates_ordering%size_ = n_labels + + end function internal_get_heating_rates_ordering + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine internal_run_tuvx(tuvx, number_of_layers, solar_zenith_angle, & diff --git a/src/tuvx/interface_grid.F90 b/src/tuvx/interface_grid.F90 index 70faad54..3328e75b 100644 --- a/src/tuvx/interface_grid.F90 +++ b/src/tuvx/interface_grid.F90 @@ -17,7 +17,6 @@ function internal_create_grid(grid_name, grid_name_length, units, & units_length, num_sections, error_code) & bind(C, name="InternalCreateGrid") result(grid) use iso_c_binding, only: c_ptr, c_f_pointer, c_char, c_loc, c_size_t, c_int - use musica_tuvx_util, only: to_f_string use musica_string, only: string_t use tuvx_grid_from_host, only: grid_from_host_t diff --git a/src/tuvx/interface_grid_map.F90 b/src/tuvx/interface_grid_map.F90 index c0eaf250..12c306a6 100644 --- a/src/tuvx/interface_grid_map.F90 +++ b/src/tuvx/interface_grid_map.F90 @@ -6,7 +6,6 @@ module tuvx_interface_grid_map use iso_c_binding, only : c_ptr, c_loc, c_int, c_size_t, c_char use tuvx_grid_warehouse, only : grid_warehouse_t use tuvx_grid, only : grid_t - use musica_tuvx_util, only : to_f_string, string_t_c use musica_string, only : string_t implicit none diff --git a/src/tuvx/interface_profile.F90 b/src/tuvx/interface_profile.F90 index eaccc102..83f3c2ac 100644 --- a/src/tuvx/interface_profile.F90 +++ b/src/tuvx/interface_profile.F90 @@ -17,7 +17,6 @@ function internal_create_profile(profile_name, profile_name_length, units, & units_length, grid_updater_c, error_code) & bind(C, name="InternalCreateProfile") result(profile) use iso_c_binding, only: c_ptr, c_f_pointer, c_char, c_loc, c_size_t, c_int - use musica_tuvx_util, only: to_f_string use musica_string, only: string_t use tuvx_grid_from_host, only: grid_updater_t use tuvx_profile_from_host, only: profile_from_host_t diff --git a/src/tuvx/interface_profile_map.F90 b/src/tuvx/interface_profile_map.F90 index 47e9cba7..6ba67a48 100644 --- a/src/tuvx/interface_profile_map.F90 +++ b/src/tuvx/interface_profile_map.F90 @@ -6,7 +6,6 @@ module tuvx_interface_profile_map use iso_c_binding, only : c_ptr, c_loc, c_int, c_size_t, c_char use tuvx_profile, only : profile_t use tuvx_profile_warehouse, only : profile_warehouse_t - use musica_tuvx_util, only : to_f_string, string_t_c use musica_string, only : string_t implicit none diff --git a/src/tuvx/interface_radiator_map.F90 b/src/tuvx/interface_radiator_map.F90 index 2ddea28e..5b3a25b0 100644 --- a/src/tuvx/interface_radiator_map.F90 +++ b/src/tuvx/interface_radiator_map.F90 @@ -6,7 +6,6 @@ module tuvx_interface_radiator_map use iso_c_binding, only : c_ptr, c_loc, c_int, c_size_t, c_char use tuvx_radiator, only : radiator_t use tuvx_radiator_warehouse, only : radiator_warehouse_t - use musica_tuvx_util, only : to_f_string, string_t_c use musica_string, only : string_t implicit none diff --git a/src/tuvx/tuvx_util.F90 b/src/tuvx/interface_util.F90 similarity index 86% rename from src/tuvx/tuvx_util.F90 rename to src/tuvx/interface_util.F90 index ad6b6ebb..a22f3ff3 100644 --- a/src/tuvx/tuvx_util.F90 +++ b/src/tuvx/interface_util.F90 @@ -1,6 +1,6 @@ ! Copyright (C) 2023-2024 National Center for Atmospheric Research ! SPDX-License-Identifier: Apache-2.0 -module musica_tuvx_util +module tuvx_interface_util use iso_c_binding, only: c_char, c_int, c_ptr, c_size_t, & c_null_ptr, c_f_pointer @@ -9,7 +9,8 @@ module musica_tuvx_util private public :: string_t_c, string_t, error_t_c, error_t, mapping_t_c, mapping_t, & - to_c_string, to_f_string, assert, delete_string_c + mappings_t_c, to_c_string, to_f_string, assert, delete_string_c, & + create_string_t_c, delete_string_t_c, allocate_mappings_c !> Wrapper for a c string type, bind(c) :: string_t_c @@ -81,12 +82,22 @@ module musica_tuvx_util module procedure mapping_constructor end interface mapping_t + type, bind(c) :: mappings_t_c + type(c_ptr) :: mappings_ = c_null_ptr + integer(c_size_t) :: size_ = 0_c_size_t + end type mappings_t_c + interface function create_string_c( string ) bind(c, name="CreateString") import :: string_t_c, c_char character(kind=c_char, len=1), intent(in) :: string(*) type(string_t_c) :: create_string_c end function create_string_c + function allocate_mappings_c( size ) bind(c, name="AllocateMappingArray") + import :: c_size_t, c_ptr + integer(c_size_t), value, intent(in) :: size + type(c_ptr) :: allocate_mappings_c + end function allocate_mappings_c subroutine delete_string_c( string ) bind(c, name="DeleteString") import :: string_t_c type(string_t_c), intent(inout) :: string @@ -312,6 +323,39 @@ function to_c_string( f_string ) result( c_string ) end function to_c_string +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Creates a new string_t_c from a fortran character array + function create_string_t_c( f_string ) result( c_string ) + + character(len=*), intent(in) :: f_string + type(string_t_c) :: c_string + + c_string = create_string_c( to_c_string( f_string ) ) + + end function create_string_t_c + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Delete a string_t_c object + subroutine delete_string_t_c( string ) & + bind(c, name="InternalDeleteString") + use iso_c_binding, only : c_char, c_null_ptr, c_size_t, & + c_f_pointer, c_associated + + type(string_t_c), intent(inout) :: string + + character(kind=c_char, len=1), pointer :: c_string_ptr(:) + + if ( c_associated( string%ptr_ ) ) then + call c_f_pointer( string%ptr_, c_string_ptr, [ string%size_ + 1 ] ) + deallocate( c_string_ptr ) + string%ptr_ = c_null_ptr + string%size_ = 0_c_size_t + end if + + end subroutine delete_string_t_c + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Convert a c string to a fortran character array @@ -361,4 +405,4 @@ end subroutine assert !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end module musica_tuvx_util \ No newline at end of file +end module tuvx_interface_util \ No newline at end of file diff --git a/src/tuvx/tuvx.cpp b/src/tuvx/tuvx.cpp index 3a2b3a9b..594ef880 100644 --- a/src/tuvx/tuvx.cpp +++ b/src/tuvx/tuvx.cpp @@ -66,6 +66,18 @@ namespace musica return tuvx->CreateRadiatorMap(error); } + Mappings GetPhotolysisRateConstantsOrdering(TUVX *tuvx, Error *error) + { + DeleteError(error); + return tuvx->GetPhotolysisRateConstantsOrdering(error); + } + + Mappings GetHeatingRatesOrdering(TUVX *tuvx, Error *error) + { + DeleteError(error); + return tuvx->GetHeatingRatesOrdering(error); + } + void RunTuvx( TUVX *const tuvx, const double solar_zenith_angle, @@ -172,6 +184,30 @@ namespace musica return radiator_map; } + Mappings TUVX::GetPhotolysisRateConstantsOrdering(Error *error) + { + *error = NoError(); + int error_code = 0; + Mappings mappings = InternalGetPhotolysisRateConstantsOrdering(tuvx_, &error_code); + if (error_code != 0) + { + *error = Error{ 1, CreateString(MUSICA_ERROR_CATEGORY), CreateString("Failed to get photolysis rate constants ordering") }; + } + return mappings; + } + + Mappings TUVX::GetHeatingRatesOrdering(Error *error) + { + *error = NoError(); + int error_code = 0; + Mappings mappings = InternalGetHeatingRatesOrdering(tuvx_, &error_code); + if (error_code != 0) + { + *error = Error{ 1, CreateString(MUSICA_ERROR_CATEGORY), CreateString("Failed to get heating rates ordering") }; + } + return mappings; + } + void TUVX::Run( const double solar_zenith_angle, const double earth_sun_distance, diff --git a/src/util.cpp b/src/util.cpp index abf70657..5ea9036f 100644 --- a/src/util.cpp +++ b/src/util.cpp @@ -139,6 +139,11 @@ namespace musica return mapping; } + Mapping* AllocateMappingArray(const std::size_t size) + { + return new Mapping[size]; + } + Mappings CreateMappings(std::size_t size) { Mappings mappings; From 186bdfeda84565ebeae0d63dd20579e27711997f Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Mon, 4 Nov 2024 09:20:29 -0600 Subject: [PATCH 04/13] Auto-format code using Clang-Format (#236) Co-authored-by: GitHub Actions --- src/test/unit/tuvx/tuvx_run_from_config.cpp | 2 +- src/tuvx/tuvx.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/test/unit/tuvx/tuvx_run_from_config.cpp b/src/test/unit/tuvx/tuvx_run_from_config.cpp index 8cb440fd..4f6cb07a 100644 --- a/src/test/unit/tuvx/tuvx_run_from_config.cpp +++ b/src/test/unit/tuvx/tuvx_run_from_config.cpp @@ -9,7 +9,7 @@ using namespace musica; const double expected_photolysis_rate_constants[3][4] = { { 8.91393763338872e-28, 1.64258192104497e-20, 8.48391527327371e-14, 9.87420948924703e-08 }, { 2.49575956372508e-27, 4.58686176250519e-20, 2.22679622672858e-13, 2.29392676897831e-07 }, - { 1.78278752667774e-27, 3.28516384208994e-20, 1.69678305465474e-13, 1.97484189784941e-07} + { 1.78278752667774e-27, 3.28516384208994e-20, 1.69678305465474e-13, 1.97484189784941e-07 } }; const double expected_heating_rates[2][4] = { { 1.12394047546984e-46, 2.04518267143613e-39, 7.44349752571804e-33, 5.42628100199216e-28 }, diff --git a/src/tuvx/tuvx.cpp b/src/tuvx/tuvx.cpp index 594ef880..0e5e4ab1 100644 --- a/src/tuvx/tuvx.cpp +++ b/src/tuvx/tuvx.cpp @@ -191,7 +191,8 @@ namespace musica Mappings mappings = InternalGetPhotolysisRateConstantsOrdering(tuvx_, &error_code); if (error_code != 0) { - *error = Error{ 1, CreateString(MUSICA_ERROR_CATEGORY), CreateString("Failed to get photolysis rate constants ordering") }; + *error = + Error{ 1, CreateString(MUSICA_ERROR_CATEGORY), CreateString("Failed to get photolysis rate constants ordering") }; } return mappings; } From c5c12a868e19c8dc888758009bd627c107286e46 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 4 Nov 2024 11:59:29 -0800 Subject: [PATCH 05/13] initialize config pointer to nullptr (#237) --- include/musica/util.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/musica/util.hpp b/include/musica/util.hpp index ce968479..33f664a5 100644 --- a/include/musica/util.hpp +++ b/include/musica/util.hpp @@ -41,7 +41,7 @@ namespace musica /// @brief A set of configuration data struct Configuration { - Yaml* data_; + Yaml* data_ = nullptr; }; /// @brief A struct to represent a mapping between a string and an index From 5b2e2f5b4afb5f7cbdc7112407d31b447bd891dd Mon Sep 17 00:00:00 2001 From: Jiwon Gim <55209567+boulderdaze@users.noreply.github.com> Date: Mon, 4 Nov 2024 14:11:38 -0700 Subject: [PATCH 06/13] Disable TUV-x tests when fetching content for TUV-x. (#239) * add tuvx fetch content options * fix a bug --- cmake/dependencies.cmake | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index b5c039bc..6a3e3b65 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -94,6 +94,9 @@ if (MUSICA_ENABLE_TUVX AND MUSICA_BUILD_C_CXX_INTERFACE) GIT_PROGRESS NOT ${FETCHCONTENT_QUIET} ) + set(TUVX_ENABLE_TESTS OFF) + set(TUVX_ENABLE_REGRESSION_TESTS OFF) + FetchContent_MakeAvailable(tuvx) endif() From 5817b5a4bd80f617856a6e5f1b1eefa062880f2b Mon Sep 17 00:00:00 2001 From: David Fillmore <1524012+dwfncar@users.noreply.github.com> Date: Tue, 5 Nov 2024 17:31:43 +0000 Subject: [PATCH 07/13] 195 backward euler python (#234) * Updated python wrapper to include Backward Euler options. * Added TestAnalyticalBackwardEuler. * Added TestAnalyticalStandardBackwardEuler. * Added multiple grid cell backward Euler tests. * Added places argument with default of 5. * Fixed wrapper. --------- Co-authored-by: David Fillmore --- python/test/test_analytical.py | 68 ++++++++++++++++++++++++++-------- python/wrapper.cpp | 22 ++++++++++- 2 files changed, 72 insertions(+), 18 deletions(-) diff --git a/python/test/test_analytical.py b/python/test/test_analytical.py index 07865bbc..fd344143 100644 --- a/python/test/test_analytical.py +++ b/python/test/test_analytical.py @@ -4,7 +4,7 @@ import random -def TestSingleGridCell(self, solver): +def TestSingleGridCell(self, solver, places=5): num_grid_cells = 1 time_step = 200.0 temperature = 272.5 @@ -78,17 +78,17 @@ def TestSingleGridCell(self, solver): (1.0 + (k1 * np.exp(-k2 * curr_time) - k2 * np.exp(-k1 * curr_time)) / (k2 - k1)) self.assertAlmostEqual( - ordered_concentrations[species_ordering["A"]], A_conc, places=5) + ordered_concentrations[species_ordering["A"]], A_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[species_ordering["B"]], B_conc, places=5) + ordered_concentrations[species_ordering["B"]], B_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[species_ordering["C"]], C_conc, places=5) + ordered_concentrations[species_ordering["C"]], C_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[species_ordering["D"]], D_conc, places=5) + ordered_concentrations[species_ordering["D"]], D_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[species_ordering["E"]], E_conc, places=5) + ordered_concentrations[species_ordering["E"]], E_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[species_ordering["F"]], F_conc, places=5) + ordered_concentrations[species_ordering["F"]], F_conc, places=places) curr_time += time_step @@ -102,7 +102,7 @@ def test_simulation(self): TestSingleGridCell(self, solver) -class TestAnalyicalStandardRosenbrock(unittest.TestCase): +class TestAnalyticalStandardRosenbrock(unittest.TestCase): def test_simulation(self): solver = musica.create_solver( "configs/analytical", @@ -111,7 +111,25 @@ def test_simulation(self): TestSingleGridCell(self, solver) -def TestMultipleGridCell(self, solver, num_grid_cells): +class TestAnalyticalBackwardEuler(unittest.TestCase): + def test_simulation(self): + solver = musica.create_solver( + "configs/analytical", + musica.micmsolver.backward_euler, + 1) + TestSingleGridCell(self, solver, places=2) + + +class TestAnalyticalStandardBackwardEuler(unittest.TestCase): + def test_simulation(self): + solver = musica.create_solver( + "configs/analytical", + musica.micmsolver.backward_euler_standard_order, + 1) + TestSingleGridCell(self, solver, places=2) + + +def TestMultipleGridCell(self, solver, num_grid_cells, places=5): time_step = 200.0 temperatures = [] pressures = [] @@ -221,17 +239,17 @@ def TestMultipleGridCell(self, solver, num_grid_cells): k1[i] * np.exp(-k2[i] * curr_time) - k2[i] * np.exp(-k1[i] * curr_time)) / (k2[i] - k1[i])) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["A"]], A_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["A"]], A_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["B"]], B_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["B"]], B_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["C"]], C_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["C"]], C_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["D"]], D_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["D"]], D_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["E"]], E_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["E"]], E_conc, places=places) self.assertAlmostEqual( - ordered_concentrations[i * len(concentrations.keys()) + species_ordering["F"]], F_conc, places=5) + ordered_concentrations[i * len(concentrations.keys()) + species_ordering["F"]], F_conc, places=places) curr_time += time_step @@ -245,7 +263,7 @@ def test_simulation(self): TestMultipleGridCell(self, solver, 3) -class TestAnalyicalStandardRosenbrockMultipleGridCells(unittest.TestCase): +class TestAnalyticalStandardRosenbrockMultipleGridCells(unittest.TestCase): def test_simulation(self): solver = musica.create_solver( "configs/analytical", @@ -254,5 +272,23 @@ def test_simulation(self): TestMultipleGridCell(self, solver, 3) +class TestAnalyticalBackwardEulerMultipleGridCells(unittest.TestCase): + def test_simulation(self): + solver = musica.create_solver( + "configs/analytical", + musica.micmsolver.backward_euler, + 3) + TestMultipleGridCell(self, solver, 3, places=2) + + +class TestAnalyticalStandardBackwardEulerMultipleGridCells(unittest.TestCase): + def test_simulation(self): + solver = musica.create_solver( + "configs/analytical", + musica.micmsolver.backward_euler_standard_order, + 3) + TestMultipleGridCell(self, solver, 3, places=2) + + if __name__ == '__main__': unittest.main() diff --git a/python/wrapper.cpp b/python/wrapper.cpp index d85ebeca..5e3f032f 100644 --- a/python/wrapper.cpp +++ b/python/wrapper.cpp @@ -14,7 +14,9 @@ PYBIND11_MODULE(musica, m) py::enum_(m, "micmsolver") .value("rosenbrock", musica::MICMSolver::Rosenbrock) - .value("rosenbrock_standard_order", musica::MICMSolver::RosenbrockStandardOrder); + .value("rosenbrock_standard_order", musica::MICMSolver::RosenbrockStandardOrder) + .value("backward_euler", musica::MICMSolver::BackwardEuler) + .value("backward_euler_standard_order", musica::MICMSolver::BackwardEulerStandardOrder); m.def( "create_solver", @@ -162,6 +164,14 @@ PYBIND11_MODULE(musica, m) { map = micm->rosenbrock_standard_.second.variable_map_; } + else if (micm->solver_type_ == musica::MICMSolver::BackwardEuler) + { + map = micm->backward_euler_.second.variable_map_; + } + else if (micm->solver_type_ == musica::MICMSolver::BackwardEulerStandardOrder) + { + map = micm->backward_euler_standard_.second.variable_map_; + } return map; }, @@ -182,8 +192,16 @@ PYBIND11_MODULE(musica, m) { map = micm->rosenbrock_standard_.second.custom_rate_parameter_map_; } + else if (micm->solver_type_ == musica::MICMSolver::BackwardEuler) + { + map = micm->backward_euler_.second.custom_rate_parameter_map_; + } + else if (micm->solver_type_ == musica::MICMSolver::BackwardEulerStandardOrder) + { + map = micm->backward_euler_standard_.second.custom_rate_parameter_map_; + } return map; }, "Return map of reaction rates"); -} \ No newline at end of file +} From f7a69ef440aabaa74eb1786830250d50f1782b6e Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 5 Nov 2024 10:31:14 -0800 Subject: [PATCH 08/13] Use fortran types in API functions (#238) * use fortran types in API functions * add fortran MICM solve function that takes c pointers * update python wrapper and address reviewer comments * fix fortran API test * address reviewer comments --- CMakeLists.txt | 2 +- fortran/micm.F90 | 166 +++++++-- .../fetch_content_integration/CMakeLists.txt | 2 + .../test_micm_api.F90 | 326 +++++++++++++----- .../test_micm_box_model.F90 | 86 ++++- fortran/test/unit/util.F90 | 6 + fortran/util.F90 | 29 +- include/musica/micm.hpp | 5 +- python/test/test_analytical.py | 153 ++++++-- python/test/test_chapman.py | 2 +- src/micm/micm.cpp | 32 +- src/test/unit/micm/micm_c_api.cpp | 201 +++++++++-- 12 files changed, 779 insertions(+), 231 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 45b1c62d..a0fb629a 100644 --- 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 index 22b41039..92383920 100644 --- a/fortran/micm.F90 +++ b/fortran/micm.F90 @@ -10,7 +10,7 @@ module musica_micm implicit none public :: micm_t, solver_stats_t, get_micm_version - public :: Rosenbrock, RosenbrockStandardOrder, BackwardEuler, BackwardEulerStandardOrder + public :: UndefinedSolver, Rosenbrock, RosenbrockStandardOrder, BackwardEuler, BackwardEulerStandardOrder private !> Wrapper for c solver stats @@ -28,6 +28,7 @@ module musica_micm ! We could use Fortran 2023 enum type feature if Fortran 2023 is supported ! https://fortran-lang.discourse.group/t/enumerator-type-in-bind-c-derived-type-best-practice/5947/2 enum, bind(c) + enumerator :: UndefinedSolver = 0 enumerator :: Rosenbrock = 1 enumerator :: RosenbrockStandardOrder = 2 enumerator :: BackwardEuler = 3 @@ -135,10 +136,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 = UndefinedSolver 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 @@ -191,9 +196,11 @@ function constructor(config_path, solver_type, num_grid_cells, error) result( t use musica_util, only: error_t_c, error_t, copy_mappings type(micm_t), pointer :: this character(len=*), intent(in) :: config_path - integer(c_int), intent(in) :: solver_type - integer(c_int), intent(in) :: num_grid_cells + integer, intent(in) :: solver_type + integer, intent(in) :: num_grid_cells type(error_t), intent(inout) :: error + + ! local variables character(len=1, kind=c_char) :: c_config_path(len_trim(config_path)+1) integer :: n, i type(error_t_c) :: error_c @@ -206,7 +213,10 @@ function constructor(config_path, solver_type, num_grid_cells, error) result( t end do c_config_path(n+1) = c_null_char - this%ptr = create_micm_c(c_config_path, solver_type, num_grid_cells, error_c) + 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) if (.not. error%is_success()) then deallocate(this) @@ -233,41 +243,121 @@ 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 + class(micm_t), intent(in) :: this + real(real64), intent(in) :: time_step + 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(:,:) + 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 + + 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), & + 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_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) :: this - real(c_double), intent(in) :: time_step - real(c_double), target, intent(in) :: temperature(:) - real(c_double), target, intent(in) :: pressure(:) - real(c_double), target, intent(in) :: air_density(:) - real(c_double), target, intent(inout) :: concentrations(:) - real(c_double), 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 + 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 - type(c_ptr) :: temperature_c, pressure_c, air_density_c, concentrations_c, & - user_defined_reaction_rates_c - - temperature_c = c_loc(temperature) - pressure_c = c_loc(pressure) - air_density_c = c_loc(air_density) - concentrations_c = c_loc(concentrations) - user_defined_reaction_rates_c = c_loc(user_defined_reaction_rates) - call micm_solve_c(this%ptr, time_step, temperature_c, pressure_c, air_density_c, & - concentrations_c, user_defined_reaction_rates_c, solver_state_c, & - solver_stats_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 + 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 ) @@ -359,8 +449,9 @@ end function solver_stats_t_solves !> Get the final time the solver iterated to function solver_stats_t_final_time( this ) result( final_time ) + use iso_fortran_env, only: real64 class(solver_stats_t), intent(in) :: this - real :: final_time + real(real64) :: final_time final_time = this%final_time_ @@ -382,15 +473,16 @@ function get_species_property_string(this, species_name, property_name, error) r end function get_species_property_string function get_species_property_double(this, species_name, property_name, error) result(value) + use iso_fortran_env, only: real64 use musica_util, only: error_t_c, error_t, to_c_string class(micm_t) :: this character(len=*), intent(in) :: species_name, property_name type(error_t), intent(inout) :: error - real(c_double) :: value + real(real64) :: value type(error_t_c) :: error_c - value = get_species_property_double_c(this%ptr, & - to_c_string(species_name), to_c_string(property_name), error_c) + value = real( get_species_property_double_c( this%ptr, to_c_string(species_name), & + to_c_string(property_name), error_c ), kind=real64 ) error = error_t(error_c) end function get_species_property_double @@ -399,11 +491,11 @@ function get_species_property_int(this, species_name, property_name, error) resu class(micm_t) :: this character(len=*), intent(in) :: species_name, property_name type(error_t), intent(inout) :: error - integer(c_int) :: value + integer :: value type(error_t_c) :: error_c - value = get_species_property_int_c(this%ptr, & - to_c_string(species_name), to_c_string(property_name), error_c) + value = int( get_species_property_int_c(this%ptr, & + to_c_string(species_name), to_c_string(property_name), error_c) ) error = error_t(error_c) end function get_species_property_int diff --git a/fortran/test/fetch_content_integration/CMakeLists.txt b/fortran/test/fetch_content_integration/CMakeLists.txt index d0707f7b..f6fda9b2 100644 --- 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 $ diff --git a/fortran/test/fetch_content_integration/test_micm_api.F90 b/fortran/test/fetch_content_integration/test_micm_api.F90 index fb0b634d..7168f853 100644 --- a/fortran/test/fetch_content_integration/test_micm_api.F90 +++ b/fortran/test/fetch_content_integration/test_micm_api.F90 @@ -5,6 +5,7 @@ program test_micm_api use, intrinsic :: iso_c_binding use, intrinsic :: ieee_arithmetic + use iso_fortran_env, only: real64 use musica_micm, only: micm_t, solver_stats_t, get_micm_version use musica_micm, only: Rosenbrock, RosenbrockStandardOrder, BackwardEuler, BackwardEulerStandardOrder use musica_util, only: assert, error_t, mapping_t, string_t, find_mapping_index @@ -16,14 +17,18 @@ program test_micm_api #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__ ) +#ifndef MICM_VECTOR_MATRIX_SIZE + #define MICM_VECTOR_MATRIX_SIZE 4 +#endif + implicit none type :: ArrheniusReaction - real(c_double) :: A_ = 1.0 - real(c_double) :: B_ = 0.0 - real(c_double) :: C_ = 0.0 - real(c_double) :: D_ = 300.0 - real(c_double) :: E_ = 0.0 + real(real64) :: A_ = 1.0 + real(real64) :: B_ = 0.0 + real(real64) :: C_ = 0.0 + real(real64) :: D_ = 300.0 + real(real64) :: E_ = 0.0 end type ArrheniusReaction call test_api() @@ -36,9 +41,9 @@ program test_micm_api function calculate_arrhenius( reaction, temperature, pressure ) result( rate ) type(ArrheniusReaction), intent(in) :: reaction - real(c_double), intent(in) :: temperature - real(c_double), intent(in) :: pressure - real(c_double) :: rate + real(real64), intent(in) :: temperature + real(real64), intent(in) :: pressure + real(real64) :: rate rate = reaction%A_ * exp( reaction%C_ / temperature ) & * (temperature / reaction%D_) ** reaction%B_ & * (1.0 + reaction%E_ * pressure) @@ -46,31 +51,31 @@ end function calculate_arrhenius subroutine test_api() - type(string_t) :: micm_version - type(micm_t), pointer :: micm - real(c_double) :: time_step - real(c_double), target :: temperature(1) - real(c_double), target :: pressure(1) - real(c_double), target :: air_density(1) - real(c_double), target, dimension(4) :: concentrations - real(c_double), target, dimension(3) :: user_defined_reaction_rates - character(len=256) :: config_path - integer(c_int) :: solver_type - integer(c_int) :: num_grid_cells - 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 - integer :: O2_index, O_index, O1D_index, O3_index - integer :: jO2_index, jO3a_index, jO3b_index + type(string_t) :: micm_version + type(micm_t), pointer :: micm + real(real64) :: time_step + real(real64) :: temperature(1) + real(real64) :: pressure(1) + real(real64) :: air_density(1) + real(real64), dimension(4,1) :: concentrations + real(real64), dimension(3,1) :: user_defined_reaction_rates + character(len=256) :: config_path + integer :: solver_type + integer :: num_grid_cells + character(len=:), allocatable :: string_value + real(real64) :: 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(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 + integer :: i + integer :: O2_index, O_index, O1D_index, O3_index + integer :: jO2_index, jO3a_index, jO3b_index config_path = "configs/chapman" - solver_type = Rosenbrock + solver_type = RosenbrockStandardOrder num_grid_cells = 1 time_step = 200 @@ -98,13 +103,13 @@ subroutine test_api() pressure(1) = 101253.4 air_density(:) = pressure(:) / ( GAS_CONSTANT * temperature(:) ) - concentrations(O2_index) = 0.75 - concentrations(O_index) = 0.0 - concentrations(O1D_index) = 0.0 - concentrations(O3_index) = 0.0000081 - user_defined_reaction_rates(jO2_index) = 2.7e-19 - user_defined_reaction_rates(jO3a_index) = 1.13e-9 - user_defined_reaction_rates(jO3b_index) = 5.8e-8 + concentrations(O2_index,1) = 0.75 + concentrations(O_index,1) = 0.0 + concentrations(O1D_index,1) = 0.0 + concentrations(O3_index,1) = 0.0000081 + user_defined_reaction_rates(jO2_index,1) = 2.7e-19 + user_defined_reaction_rates(jO3a_index,1) = 1.13e-9 + user_defined_reaction_rates(jO3b_index,1) = 5.8e-8 micm_version = get_micm_version() print *, "[test micm fort api] MICM version ", micm_version%get_char_array() @@ -143,10 +148,10 @@ subroutine test_api() deallocate( string_value ) double_value = micm%get_species_property_double( "O3", "molecular weight [kg mol-1]", error ) ASSERT( error%is_success() ) - ASSERT_EQ( double_value, 0.048_c_double ) + ASSERT_EQ( double_value, 0.048_real64 ) int_value = micm%get_species_property_int( "O3", "__atoms", error ) ASSERT( error%is_success() ) - ASSERT_EQ( int_value, 3_c_int ) + ASSERT_EQ( int_value, 3 ) bool_value = micm%get_species_property_bool( "O3", "__do advect", error ) ASSERT( error%is_success() ) ASSERT( logical( bool_value ) ) @@ -168,31 +173,160 @@ subroutine test_api() end subroutine test_api - subroutine test_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accuracy) + subroutine test_vector_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accuracy) + + type(micm_t), pointer, intent(inout) :: micm + integer, intent(in) :: NUM_GRID_CELLS + real(real64), intent(in) :: time_step + real, intent(in) :: test_accuracy + + integer, parameter :: NUM_SPECIES = 6 + integer, parameter :: NUM_USER_DEFINED_REACTION_RATES = 2 + real(real64), target :: temperature(NUM_GRID_CELLS) + real(real64), target :: pressure(NUM_GRID_CELLS) + real(real64), target :: air_density(NUM_GRID_CELLS) + real(real64), target :: concentrations(NUM_GRID_CELLS,NUM_SPECIES) + real(real64), target :: concentrations_c_ptrs(NUM_GRID_CELLS,NUM_SPECIES) + real(real64), target :: initial_concentrations(NUM_GRID_CELLS,NUM_SPECIES) + real(real64), target :: user_defined_reaction_rates(NUM_GRID_CELLS,NUM_USER_DEFINED_REACTION_RATES) + type(string_t) :: solver_state + type(solver_stats_t) :: solver_stats + integer :: solver_type + type(error_t) :: error + real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 + integer :: A_index, B_index, C_index, D_index, E_index, F_index + integer :: R1_index, R2_index + real(real64) :: initial_A, initial_C, initial_D, initial_F + real(real64) :: k1, k2, k3, k4 + real(real64) :: A, B, C, D, E, F + integer :: i_cell + real :: temp + type(ArrheniusReaction) :: r1, r2 + + A_index = micm%species_ordering%index( "A", error ) + ASSERT( error%is_success() ) + B_index = micm%species_ordering%index( "B", error ) + ASSERT( error%is_success() ) + C_index = micm%species_ordering%index( "C", error ) + ASSERT( error%is_success() ) + D_index = micm%species_ordering%index( "D", error ) + ASSERT( error%is_success() ) + E_index = micm%species_ordering%index( "E", error ) + ASSERT( error%is_success() ) + F_index = micm%species_ordering%index( "F", error ) + ASSERT( error%is_success() ) + + R1_index = micm%user_defined_reaction_rates%index( "USER.reaction 1", error ) + ASSERT( error%is_success() ) + R2_index = micm%user_defined_reaction_rates%index( "USER.reaction 2", error ) + ASSERT( error%is_success() ) + + do i_cell = 1, NUM_GRID_CELLS + call random_number( temp ) + temperature(i_cell) = 265.0 + temp * 20.0 + call random_number( temp ) + pressure(i_cell) = 100753.3 + temp * 1000.0 + air_density(i_cell) = pressure(i_cell) / ( GAS_CONSTANT * temperature(i_cell) ) + call random_number( temp ) + concentrations(i_cell,A_index) = 0.7 + temp * 0.1 + concentrations(i_cell,B_index) = 0.0 + call random_number( temp ) + concentrations(i_cell,C_index) = 0.35 + temp * 0.1 + call random_number( temp ) + concentrations(i_cell,D_index) = 0.75 + temp * 0.1 + concentrations(i_cell,E_index) = 0.0 + call random_number( temp ) + concentrations(i_cell,F_index) = 0.05 + temp * 0.1 + call random_number( temp ) + user_defined_reaction_rates(i_cell,R1_index) = 0.0005 + temp * 0.0001 + call random_number( temp ) + user_defined_reaction_rates(i_cell,R2_index) = 0.0015 + temp * 0.0001 + end do + initial_concentrations(:,:) = concentrations(:,:) + concentrations_c_ptrs(:,:) = concentrations(:,:) + + ! solve by passing fortran arrays + call micm%solve(time_step, temperature, pressure, air_density, concentrations, & + user_defined_reaction_rates, solver_state, solver_stats, error) + ASSERT( error%is_success() ) + ASSERT_EQ(solver_state%get_char_array(), "Converged") + + ! solve by passing C pointers + call micm%solve(time_step, c_loc(temperature), c_loc(pressure), c_loc(air_density), & + c_loc(concentrations_c_ptrs), c_loc(user_defined_reaction_rates), & + solver_state, solver_stats, error) + ASSERT( error%is_success() ) + ASSERT_EQ(solver_state%get_char_array(), "Converged") + + ! check concentrations + do i_cell = 1, NUM_GRID_CELLS + ASSERT_EQ(concentrations(i_cell,A_index), concentrations_c_ptrs(i_cell,A_index)) + ASSERT_EQ(concentrations(i_cell,B_index), concentrations_c_ptrs(i_cell,B_index)) + ASSERT_EQ(concentrations(i_cell,C_index), concentrations_c_ptrs(i_cell,C_index)) + ASSERT_EQ(concentrations(i_cell,D_index), concentrations_c_ptrs(i_cell,D_index)) + ASSERT_EQ(concentrations(i_cell,E_index), concentrations_c_ptrs(i_cell,E_index)) + ASSERT_EQ(concentrations(i_cell,F_index), concentrations_c_ptrs(i_cell,F_index)) + end do + + r1%A_ = 0.004 + r1%C_ = 50.0 + r2%A_ = 0.012 + r2%B_ = -2.0 + r2%C_ = 75.0 + r2%D_ = 50.0 + r2%E_ = 1.0e-6 + + do i_cell = 1, NUM_GRID_CELLS + initial_A = initial_concentrations(i_cell,A_index) + initial_C = initial_concentrations(i_cell,C_index) + initial_D = initial_concentrations(i_cell,D_index) + initial_F = initial_concentrations(i_cell,F_index) + k1 = user_defined_reaction_rates(i_cell,R1_index) + k2 = user_defined_reaction_rates(i_cell,R2_index) + k3 = calculate_arrhenius( r1, temperature(i_cell), pressure(i_cell) ) + k4 = calculate_arrhenius( r2, temperature(i_cell), pressure(i_cell) ) + A = initial_A * exp( -k3 * time_step ) + B = initial_A * (k3 / (k4 - k3)) * (exp(-k3 * time_step) - exp(-k4 * time_step)) + C = initial_C + initial_A * (1.0 + (k3 * exp(-k4 * time_step) - k4 * exp(-k3 * time_step)) / (k4 - k3)) + D = initial_D * exp( -k1 * time_step ) + E = initial_D * (k1 / (k2 - k1)) * (exp(-k1 * time_step) - exp(-k2 * time_step)) + F = initial_F + initial_D * (1.0 + (k1 * exp(-k2 * time_step) - k2 * exp(-k1 * time_step)) / (k2 - k1)) + ASSERT_NEAR(concentrations(i_cell,A_index), A, test_accuracy) + ASSERT_NEAR(concentrations(i_cell,B_index), B, test_accuracy) + ASSERT_NEAR(concentrations(i_cell,C_index), C, test_accuracy) + ASSERT_NEAR(concentrations(i_cell,D_index), D, test_accuracy) + ASSERT_NEAR(concentrations(i_cell,E_index), E, test_accuracy) + ASSERT_NEAR(concentrations(i_cell,F_index), F, test_accuracy) + end do + + end subroutine test_vector_multiple_grid_cells + + subroutine test_standard_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accuracy) type(micm_t), pointer, intent(inout) :: micm integer, intent(in) :: NUM_GRID_CELLS - real(c_double), intent(in) :: time_step + real(real64), intent(in) :: time_step real, intent(in) :: test_accuracy integer, parameter :: NUM_SPECIES = 6 integer, parameter :: NUM_USER_DEFINED_REACTION_RATES = 2 - real(c_double), target :: temperature(NUM_GRID_CELLS) - real(c_double), target :: pressure(NUM_GRID_CELLS) - real(c_double), target :: air_density(NUM_GRID_CELLS) - real(c_double), target :: concentrations(NUM_GRID_CELLS * NUM_SPECIES) - real(c_double), target :: initial_concentrations(NUM_GRID_CELLS * NUM_SPECIES) - real(c_double), target :: user_defined_reaction_rates(NUM_GRID_CELLS * NUM_USER_DEFINED_REACTION_RATES) + real(real64), target :: temperature(NUM_GRID_CELLS) + real(real64), target :: pressure(NUM_GRID_CELLS) + real(real64), target :: air_density(NUM_GRID_CELLS) + real(real64), target :: concentrations(NUM_SPECIES, NUM_GRID_CELLS) + real(real64), target :: concentrations_c_ptrs(NUM_SPECIES, NUM_GRID_CELLS) + real(real64), target :: initial_concentrations(NUM_SPECIES, NUM_GRID_CELLS) + real(real64), target :: user_defined_reaction_rates(NUM_USER_DEFINED_REACTION_RATES, NUM_GRID_CELLS) type(string_t) :: solver_state type(solver_stats_t) :: solver_stats integer(c_int) :: solver_type type(error_t) :: error - real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1 + real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 integer :: A_index, B_index, C_index, D_index, E_index, F_index integer :: R1_index, R2_index - real(c_double) :: initial_A, initial_C, initial_D, initial_F - real(c_double) :: k1, k2, k3, k4 - real(c_double) :: A, B, C, D, E, F + real(real64) :: initial_A, initial_C, initial_D, initial_F + real(real64) :: k1, k2, k3, k4 + real(real64) :: A, B, C, D, E, F integer :: i_cell real :: temp type(ArrheniusReaction) :: r1, r2 @@ -222,28 +356,46 @@ subroutine test_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accura pressure(i_cell) = 100753.3 + temp * 1000.0 air_density(i_cell) = pressure(i_cell) / ( GAS_CONSTANT * temperature(i_cell) ) call random_number( temp ) - concentrations((i_cell-1)*NUM_SPECIES+A_index) = 0.7 + temp * 0.1 - concentrations((i_cell-1)*NUM_SPECIES+B_index) = 0.0 + concentrations(A_index, i_cell) = 0.7 + temp * 0.1 + concentrations(B_index, i_cell) = 0.0 call random_number( temp ) - concentrations((i_cell-1)*NUM_SPECIES+C_index) = 0.35 + temp * 0.1 + concentrations(C_index, i_cell) = 0.35 + temp * 0.1 call random_number( temp ) - concentrations((i_cell-1)*NUM_SPECIES+D_index) = 0.75 + temp * 0.1 - concentrations((i_cell-1)*NUM_SPECIES+E_index) = 0.0 + concentrations(D_index, i_cell) = 0.75 + temp * 0.1 + concentrations(E_index, i_cell) = 0.0 call random_number( temp ) - concentrations((i_cell-1)*NUM_SPECIES+F_index) = 0.05 + temp * 0.1 + concentrations(F_index, i_cell) = 0.05 + temp * 0.1 call random_number( temp ) - user_defined_reaction_rates((i_cell-1)*NUM_USER_DEFINED_REACTION_RATES+R1_index) = 0.0005 + temp * 0.0001 + user_defined_reaction_rates(R1_index, i_cell) = 0.0005 + temp * 0.0001 call random_number( temp ) - user_defined_reaction_rates((i_cell-1)*NUM_USER_DEFINED_REACTION_RATES+R2_index) = 0.0015 + temp * 0.0001 + user_defined_reaction_rates(R2_index, i_cell) = 0.0015 + temp * 0.0001 end do - initial_concentrations(:) = concentrations(:) + initial_concentrations(:,:) = concentrations(:,:) + concentrations_c_ptrs(:,:) = concentrations(:,:) + ! solve by passing fortran arrays call micm%solve(time_step, temperature, pressure, air_density, concentrations, & user_defined_reaction_rates, solver_state, solver_stats, error) ASSERT( error%is_success() ) + ASSERT_EQ(solver_state%get_char_array(), "Converged") + ! solve by passing C pointers + call micm%solve(time_step, c_loc(temperature), c_loc(pressure), c_loc(air_density), & + c_loc(concentrations_c_ptrs), c_loc(user_defined_reaction_rates), & + solver_state, solver_stats, error) + ASSERT( error%is_success() ) ASSERT_EQ(solver_state%get_char_array(), "Converged") + ! check concentrations + do i_cell = 1, NUM_GRID_CELLS + ASSERT_EQ(concentrations(A_index, i_cell), concentrations_c_ptrs(A_index, i_cell)) + ASSERT_EQ(concentrations(B_index, i_cell), concentrations_c_ptrs(B_index, i_cell)) + ASSERT_EQ(concentrations(C_index, i_cell), concentrations_c_ptrs(C_index, i_cell)) + ASSERT_EQ(concentrations(D_index, i_cell), concentrations_c_ptrs(D_index, i_cell)) + ASSERT_EQ(concentrations(E_index, i_cell), concentrations_c_ptrs(E_index, i_cell)) + ASSERT_EQ(concentrations(F_index, i_cell), concentrations_c_ptrs(F_index, i_cell)) + end do + r1%A_ = 0.004 r1%C_ = 50.0 r2%A_ = 0.012 @@ -253,12 +405,12 @@ subroutine test_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accura r2%E_ = 1.0e-6 do i_cell = 1, NUM_GRID_CELLS - initial_A = initial_concentrations((i_cell-1)*NUM_SPECIES+A_index) - initial_C = initial_concentrations((i_cell-1)*NUM_SPECIES+C_index) - initial_D = initial_concentrations((i_cell-1)*NUM_SPECIES+D_index) - initial_F = initial_concentrations((i_cell-1)*NUM_SPECIES+F_index) - k1 = user_defined_reaction_rates((i_cell-1)*NUM_USER_DEFINED_REACTION_RATES+R1_index) - k2 = user_defined_reaction_rates((i_cell-1)*NUM_USER_DEFINED_REACTION_RATES+R2_index) + initial_A = initial_concentrations(A_index, i_cell) + initial_C = initial_concentrations(C_index, i_cell) + initial_D = initial_concentrations(D_index, i_cell) + initial_F = initial_concentrations(F_index, i_cell) + k1 = user_defined_reaction_rates(R1_index, i_cell) + k2 = user_defined_reaction_rates(R2_index, i_cell) k3 = calculate_arrhenius( r1, temperature(i_cell), pressure(i_cell) ) k4 = calculate_arrhenius( r2, temperature(i_cell), pressure(i_cell) ) A = initial_A * exp( -k3 * time_step ) @@ -267,15 +419,15 @@ subroutine test_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test_accura D = initial_D * exp( -k1 * time_step ) E = initial_D * (k1 / (k2 - k1)) * (exp(-k1 * time_step) - exp(-k2 * time_step)) F = initial_F + initial_D * (1.0 + (k1 * exp(-k2 * time_step) - k2 * exp(-k1 * time_step)) / (k2 - k1)) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+A_index), A, test_accuracy) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+B_index), B, test_accuracy) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+C_index), C, test_accuracy) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+D_index), D, test_accuracy) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+E_index), E, test_accuracy) - ASSERT_NEAR(concentrations((i_cell-1)*NUM_SPECIES+F_index), F, test_accuracy) + ASSERT_NEAR(concentrations(A_index, i_cell), A, test_accuracy) + ASSERT_NEAR(concentrations(B_index, i_cell), B, test_accuracy) + ASSERT_NEAR(concentrations(C_index, i_cell), C, test_accuracy) + ASSERT_NEAR(concentrations(D_index, i_cell), D, test_accuracy) + ASSERT_NEAR(concentrations(E_index, i_cell), E, test_accuracy) + ASSERT_NEAR(concentrations(F_index, i_cell), F, test_accuracy) end do - end subroutine test_multiple_grid_cells + end subroutine test_standard_multiple_grid_cells subroutine test_multiple_grid_cell_vector_Rosenbrock() @@ -283,14 +435,14 @@ subroutine test_multiple_grid_cell_vector_Rosenbrock() integer :: num_grid_cells type(error_t) :: error - real(c_double), parameter :: time_step = 200 + real(real64), parameter :: time_step = 200 real, parameter :: test_accuracy = 5.0e-3 - num_grid_cells = 3 + num_grid_cells = MICM_VECTOR_MATRIX_SIZE micm => micm_t( "configs/analytical", Rosenbrock, num_grid_cells, error ) ASSERT( error%is_success() ) - call test_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) + call test_vector_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) deallocate( micm ) @@ -302,14 +454,14 @@ subroutine test_multiple_grid_cell_standard_Rosenbrock() integer :: num_grid_cells type(error_t) :: error - real(c_double), parameter :: time_step = 200 - real, parameter :: test_accuracy = 5.0e-3 + real(real64), parameter :: time_step = 200 + real, parameter :: test_accuracy = 5.0e-3 num_grid_cells = 3 micm => micm_t( "configs/analytical", RosenbrockStandardOrder, num_grid_cells, error ) ASSERT( error%is_success() ) - call test_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) + call test_standard_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) deallocate( micm ) @@ -321,14 +473,14 @@ subroutine test_multiple_grid_cell_vector_BackwardEuler() integer :: num_grid_cells type(error_t) :: error - real(c_double), parameter :: time_step = 10 - real, parameter :: test_accuracy = 0.1 + real(real64), parameter :: time_step = 10 + real, parameter :: test_accuracy = 0.1 - num_grid_cells = 3 + num_grid_cells = MICM_VECTOR_MATRIX_SIZE micm => micm_t( "configs/analytical", BackwardEuler, num_grid_cells, error ) ASSERT( error%is_success() ) - call test_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) + call test_vector_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) deallocate( micm ) @@ -340,14 +492,14 @@ subroutine test_multiple_grid_cell_standard_BackwardEuler() integer :: num_grid_cells type(error_t) :: error - real(c_double), parameter :: time_step = 10 - real, parameter :: test_accuracy = 0.1 + real(real64), parameter :: time_step = 10 + real, parameter :: test_accuracy = 0.1 num_grid_cells = 3 micm => micm_t( "configs/analytical", BackwardEulerStandardOrder, num_grid_cells, error ) ASSERT( error%is_success() ) - call test_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) + call test_standard_multiple_grid_cells( micm, num_grid_cells, time_step, test_accuracy ) deallocate( micm ) diff --git a/fortran/test/fetch_content_integration/test_micm_box_model.F90 b/fortran/test/fetch_content_integration/test_micm_box_model.F90 index e1056884..03a17d51 100644 --- a/fortran/test/fetch_content_integration/test_micm_box_model.F90 +++ b/fortran/test/fetch_content_integration/test_micm_box_model.F90 @@ -3,30 +3,33 @@ program test_micm_box_model use, intrinsic :: iso_c_binding use, intrinsic :: ieee_arithmetic + use iso_fortran_env, only : real64 use musica_util, only: error_t, string_t, mapping_t use musica_micm, only: micm_t, solver_stats_t use musica_micm, only: Rosenbrock, RosenbrockStandardOrder implicit none - call box_model() + call box_model_arrays() + call box_model_c_ptrs() contains - subroutine box_model() + !> Runs a simple box model using the MICM solver and passing in fortran arrays + subroutine box_model_arrays() character(len=256) :: config_path - integer(c_int) :: solver_type - integer(c_int) :: num_grid_cells + integer :: solver_type + integer :: num_grid_cells - real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1 + real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 - real(c_double) :: time_step - real(c_double), target :: temperature(1) - real(c_double), target :: pressure(1) - real(c_double), target :: air_density(1) - real(c_double), target :: concentrations(6) - real(c_double), target :: user_defined_reaction_rates(2) + real(real64) :: time_step + real(real64), target :: temperature(1) + real(real64), target :: pressure(1) + real(real64), target :: air_density(1) + real(real64), target :: concentrations(1,6) + real(real64), target :: user_defined_reaction_rates(1,2) type(string_t) :: solver_state type(solver_stats_t) :: solver_stats @@ -45,8 +48,8 @@ subroutine box_model() pressure(1) = 1.0e5 air_density(:) = pressure(:) / (GAS_CONSTANT * temperature(:)) - concentrations = (/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 /) - user_defined_reaction_rates = (/ 0.001, 0.002 /) + concentrations(1,:) = (/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 /) + user_defined_reaction_rates(1,:) = (/ 0.001, 0.002 /) write(*,*) "Creating MICM solver..." micm => micm_t(config_path, solver_type, num_grid_cells, error) @@ -63,6 +66,61 @@ subroutine box_model() deallocate( micm ) - end subroutine box_model + end subroutine box_model_arrays + + !> Runs a simple box model using the MICM solver and passing in C pointers + subroutine box_model_c_ptrs() + + character(len=256) :: config_path + integer :: solver_type + integer :: num_grid_cells + + real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 + + real(real64) :: time_step + real(real64), target :: temperature(1) + real(real64), target :: pressure(1) + real(real64), target :: air_density(1) + real(real64), target :: concentrations(6) + real(real64), target :: user_defined_reaction_rates(2) + + type(string_t) :: solver_state + type(solver_stats_t) :: solver_stats + type(error_t) :: error + + type(micm_t), pointer :: micm + + integer :: i + + config_path = "configs/analytical" + solver_type = RosenbrockStandardOrder + num_grid_cells = 1 + + time_step = 200 + temperature(1) = 273.0 + pressure(1) = 1.0e5 + air_density(:) = pressure(:) / (GAS_CONSTANT * temperature(:)) + + concentrations = (/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 /) + user_defined_reaction_rates = (/ 0.001, 0.002 /) + + write(*,*) "Creating MICM solver..." + micm => micm_t(config_path, solver_type, num_grid_cells, error) + + do i = 1, micm%species_ordering%size() + print *, "Species Name:", micm%species_ordering%name(i), & + ", Index:", micm%species_ordering%index(i) + end do + + write(*,*) "Solving starts..." + call micm%solve(time_step, c_loc(temperature), c_loc(pressure), & + c_loc(air_density), c_loc(concentrations), & + c_loc(user_defined_reaction_rates), & + solver_state, solver_stats, error) + write(*,*) "After solving, concentrations", concentrations + + deallocate( micm ) + + end subroutine box_model_c_ptrs end program test_micm_box_model diff --git a/fortran/test/unit/util.F90 b/fortran/test/unit/util.F90 index 0fa09513..1efafd28 100644 --- a/fortran/test/unit/util.F90 +++ b/fortran/test/unit/util.F90 @@ -135,6 +135,12 @@ subroutine test_error_t() ASSERT_EQ( b%category(), "bar" ) ASSERT_EQ( b%message(), "foo" ) + a = error_t( 34, "qux", "quux" ) + + ASSERT_EQ( a%code(), 34 ) + ASSERT_EQ( a%category(), "qux" ) + ASSERT_EQ( a%message(), "quux" ) + end subroutine test_error_t !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/fortran/util.F90 b/fortran/util.F90 index 780a3d2c..14f97a5b 100644 --- a/fortran/util.F90 +++ b/fortran/util.F90 @@ -5,6 +5,7 @@ module musica_util use iso_c_binding, only: c_char, c_int, c_ptr, c_size_t, & c_null_ptr, c_f_pointer + use iso_fortran_env, only: real32, real64 implicit none private @@ -15,10 +16,10 @@ module musica_util create_string_c, musica_rk, musica_dk, find_mapping_index !> Single precision kind - integer, parameter :: musica_rk = kind(0.0) + integer, parameter :: musica_rk = real32 !> Double precision kind - integer, parameter :: musica_dk = kind(0.d0) + integer, parameter :: musica_dk = real64 !> Wrapper for a c string type, bind(c) :: string_t_c @@ -64,7 +65,8 @@ module musica_util end type error_t interface error_t - module procedure error_t_constructor + module procedure error_t_constructor_from_error_t_c + module procedure error_t_constructor_from_fortran end interface error_t !> Wrapper for a c configuration @@ -354,7 +356,7 @@ end subroutine configuration_finalize !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Constructor for an error_t object that takes ownership of an error_t_c - function error_t_constructor( c_error ) result( new_error ) + function error_t_constructor_from_error_t_c( c_error ) result( new_error ) type(error_t_c), intent(inout) :: c_error type(error_t) :: new_error @@ -368,7 +370,24 @@ function error_t_constructor( c_error ) result( new_error ) c_error%message_%size_ = 0_c_size_t c_error%code_ = 0_c_int - end function error_t_constructor + end function error_t_constructor_from_error_t_c + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Constructor for an error_t object from fortran data + function error_t_constructor_from_fortran( code, category, message ) & + result( new_error ) + + integer, intent(in) :: code + character(len=*), intent(in) :: category + character(len=*), intent(in) :: message + type(error_t) :: new_error + + new_error%code_ = code + new_error%category_ = category + new_error%message_ = message + + end function error_t_constructor_from_fortran !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/include/musica/micm.hpp b/include/musica/micm.hpp index bb97cc7a..a706c819 100644 --- a/include/musica/micm.hpp +++ b/include/musica/micm.hpp @@ -26,7 +26,7 @@ #include #ifndef MICM_VECTOR_MATRIX_SIZE - #define MICM_VECTOR_MATRIX_SIZE 1 + #define MICM_VECTOR_MATRIX_SIZE 4 #endif namespace musica @@ -41,7 +41,8 @@ namespace musica /// @brief Types of MICM solver enum MICMSolver { - Rosenbrock = 1, // Vector-ordered Rosenbrock solver + UndefinedSolver = 0, // Undefined solver + Rosenbrock, // Vector-ordered Rosenbrock solver RosenbrockStandardOrder, // Standard-ordered Rosenbrock solver BackwardEuler, // Vector-ordered BackwardEuler solver BackwardEulerStandardOrder, // Standard-ordered BackwardEuler solver diff --git a/python/test/test_analytical.py b/python/test/test_analytical.py index fd344143..f5fb7138 100644 --- a/python/test/test_analytical.py +++ b/python/test/test_analytical.py @@ -4,9 +4,8 @@ import random -def TestSingleGridCell(self, solver, places=5): +def TestSingleGridCell(self, solver, time_step, places=5): num_grid_cells = 1 - time_step = 200.0 temperature = 272.5 pressure = 101253.3 GAS_CONSTANT = 8.31446261815324 @@ -93,32 +92,15 @@ def TestSingleGridCell(self, solver, places=5): curr_time += time_step -class TestAnalyticalRosenbrock(unittest.TestCase): - def test_simulation(self): - solver = musica.create_solver( - "configs/analytical", - musica.micmsolver.rosenbrock, - 1) - TestSingleGridCell(self, solver) - - class TestAnalyticalStandardRosenbrock(unittest.TestCase): def test_simulation(self): solver = musica.create_solver( "configs/analytical", musica.micmsolver.rosenbrock_standard_order, 1) - TestSingleGridCell(self, solver) + TestSingleGridCell(self, solver, 200.0, 5) -class TestAnalyticalBackwardEuler(unittest.TestCase): - def test_simulation(self): - solver = musica.create_solver( - "configs/analytical", - musica.micmsolver.backward_euler, - 1) - TestSingleGridCell(self, solver, places=2) - class TestAnalyticalStandardBackwardEuler(unittest.TestCase): def test_simulation(self): @@ -126,11 +108,10 @@ def test_simulation(self): "configs/analytical", musica.micmsolver.backward_euler_standard_order, 1) - TestSingleGridCell(self, solver, places=2) + TestSingleGridCell(self, solver, 10.0, places=2) -def TestMultipleGridCell(self, solver, num_grid_cells, places=5): - time_step = 200.0 +def TestStandardMultipleGridCell(self, solver, num_grid_cells, time_step, places=5): temperatures = [] pressures = [] air_densities = [] @@ -254,13 +235,127 @@ def TestMultipleGridCell(self, solver, num_grid_cells, places=5): curr_time += time_step +def TestVectorMultipleGridCell(self, solver, num_grid_cells, time_step, places=5): + temperatures = [] + pressures = [] + air_densities = [] + concentrations = { + "A": [], + "B": [], + "C": [], + "D": [], + "E": [], + "F": [] + } + rate_constants = { + "USER.reaction 1": [], + "USER.reaction 2": [] + } + for i in range(num_grid_cells): + temperatures.append(275.0 + random.uniform(-10.0, 10.0)) + pressures.append(101253.3 + random.uniform(-500.0, 500.0)) + air_densities.append( + pressures[i] / (8.31446261815324 * temperatures[i])) + concentrations["A"].append(0.75 + random.uniform(-0.05, 0.05)) + concentrations["B"].append(0) + concentrations["C"].append(0.4 + random.uniform(-0.05, 0.05)) + concentrations["D"].append(0.8 + random.uniform(-0.05, 0.05)) + concentrations["E"].append(0) + concentrations["F"].append(0.1 + random.uniform(-0.05, 0.05)) + rate_constants["USER.reaction 1"].append( + 0.001 + random.uniform(-0.0001, 0.0001)) + rate_constants["USER.reaction 2"].append( + 0.002 + random.uniform(-0.0001, 0.0001)) + + rates_ordering = musica.user_defined_reaction_rates(solver) + species_ordering = musica.species_ordering(solver) + + ordered_rate_constants = ( + len(rate_constants.keys()) * num_grid_cells) * [0.0] + for i in range(num_grid_cells): + for key, value in rate_constants.items(): + ordered_rate_constants[i + + rates_ordering[key] * num_grid_cells] = value[i] + + ordered_concentrations = ( + len(concentrations.keys()) * num_grid_cells) * [0.0] + for i in range(num_grid_cells): + for key, value in concentrations.items(): + ordered_concentrations[i + + species_ordering[key] * num_grid_cells] = value[i] + + initial_concentrations = ordered_concentrations + + time_step = 1 + sim_length = 100 + + curr_time = time_step + initial_A = num_grid_cells * [0.0] + initial_C = num_grid_cells * [0.0] + initial_D = num_grid_cells * [0.0] + initial_F = num_grid_cells * [0.0] + for i in range(num_grid_cells): + initial_A[i] = initial_concentrations[i + species_ordering["A"] * num_grid_cells] + initial_C[i] = initial_concentrations[i + species_ordering["C"] * num_grid_cells] + initial_D[i] = initial_concentrations[i + species_ordering["D"] * num_grid_cells] + initial_F[i] = initial_concentrations[i + species_ordering["F"] * num_grid_cells] + + k1 = num_grid_cells * [0.0] + k2 = num_grid_cells * [0.0] + k3 = num_grid_cells * [0.0] + k4 = num_grid_cells * [0.0] + for i in range(num_grid_cells): + k1[i] = ordered_rate_constants[i + rates_ordering["USER.reaction 1"] * num_grid_cells] + k2[i] = ordered_rate_constants[i + rates_ordering["USER.reaction 2"] * num_grid_cells] + k3[i] = 0.004 * np.exp(50.0 / temperatures[i]) + k4[i] = 0.012 * np.exp(75.0 / temperatures[i]) * \ + (temperatures[i] / 50.0)**(-2) * (1.0 + 1.0e-6 * pressures[i]) + + while curr_time <= sim_length: + musica.micm_solve( + solver, + time_step, + temperatures, + pressures, + air_densities, + ordered_concentrations, + ordered_rate_constants) + + for i in range(num_grid_cells): + A_conc = initial_A[i] * np.exp(-(k3[i]) * curr_time) + B_conc = initial_A[i] * (k3[i] / (k4[i] - k3[i])) * \ + (np.exp(-k3[i] * curr_time) - np.exp(-k4[i] * curr_time)) + C_conc = initial_C[i] + initial_A[i] * (1.0 + ( + k3[i] * np.exp(-k4[i] * curr_time) - k4[i] * np.exp(-k3[i] * curr_time)) / (k4[i] - k3[i])) + D_conc = initial_D[i] * np.exp(-(k1[i]) * curr_time) + E_conc = initial_D[i] * (k1[i] / (k2[i] - k1[i])) * \ + (np.exp(-k1[i] * curr_time) - np.exp(-k2[i] * curr_time)) + F_conc = initial_F[i] + initial_D[i] * (1.0 + ( + k1[i] * np.exp(-k2[i] * curr_time) - k2[i] * np.exp(-k1[i] * curr_time)) / (k2[i] - k1[i])) + + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["A"] * num_grid_cells], A_conc, places=places) + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["B"] * num_grid_cells], B_conc, places=places) + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["C"] * num_grid_cells], C_conc, places=places) + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["D"] * num_grid_cells], D_conc, places=places) + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["E"] * num_grid_cells], E_conc, places=places) + self.assertAlmostEqual( + ordered_concentrations[i + species_ordering["F"] * num_grid_cells], F_conc, places=places) + + curr_time += time_step + + class TestAnalyticalRosenbrockMultipleGridCells(unittest.TestCase): def test_simulation(self): solver = musica.create_solver( "configs/analytical", musica.micmsolver.rosenbrock, - 3) - TestMultipleGridCell(self, solver, 3) + 4) + TestVectorMultipleGridCell(self, solver, 4, 200.0, 5) # The number of grid cells must equal the MICM matrix vector dimension class TestAnalyticalStandardRosenbrockMultipleGridCells(unittest.TestCase): @@ -269,7 +364,7 @@ def test_simulation(self): "configs/analytical", musica.micmsolver.rosenbrock_standard_order, 3) - TestMultipleGridCell(self, solver, 3) + TestStandardMultipleGridCell(self, solver, 3, 200.0, 5) class TestAnalyticalBackwardEulerMultipleGridCells(unittest.TestCase): @@ -277,8 +372,8 @@ def test_simulation(self): solver = musica.create_solver( "configs/analytical", musica.micmsolver.backward_euler, - 3) - TestMultipleGridCell(self, solver, 3, places=2) + 4) + TestVectorMultipleGridCell(self, solver, 4, 10.0, places=2) # The number of grid cells must equal the MICM matrix vector dimension class TestAnalyticalStandardBackwardEulerMultipleGridCells(unittest.TestCase): @@ -287,7 +382,7 @@ def test_simulation(self): "configs/analytical", musica.micmsolver.backward_euler_standard_order, 3) - TestMultipleGridCell(self, solver, 3, places=2) + TestStandardMultipleGridCell(self, solver, 3, 10.0, places=2) if __name__ == '__main__': diff --git a/python/test/test_chapman.py b/python/test/test_chapman.py index deb21ead..11038a92 100644 --- a/python/test/test_chapman.py +++ b/python/test/test_chapman.py @@ -13,7 +13,7 @@ def test_micm_solve(self): solver = musica.create_solver( "configs/chapman", - musica.micmsolver.rosenbrock, + musica.micmsolver.rosenbrock_standard_order, num_grid_cells) rate_constant_ordering = musica.user_defined_reaction_rates(solver) diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index ff037d7f..5b46c2d0 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -434,22 +434,18 @@ namespace musica const std::size_t num_species = state.variables_.NumColumns(); const std::size_t num_custom_rate_parameters = state.custom_rate_parameters_.NumColumns(); - int i_species_elem = 0; - int i_param_elem = 0; - for (int i_cell{}; i_cell < num_grid_cells_; ++i_cell) + std::size_t i_cond = 0; + for (auto &cond : state.conditions_) { - state.conditions_[i_cell].temperature_ = temperature[i_cell]; - state.conditions_[i_cell].pressure_ = pressure[i_cell]; - state.conditions_[i_cell].air_density_ = air_density[i_cell]; - for (int i_species{}; i_species < num_species; ++i_species) - { - state.variables_[i_cell][i_species] = concentrations[i_species_elem++]; - } - for (int i_param{}; i_param < num_custom_rate_parameters; ++i_param) - { - state.custom_rate_parameters_[i_cell][i_param] = custom_rate_parameters[i_param_elem++]; - } + cond.temperature_ = temperature[i_cond]; + cond.pressure_ = pressure[i_cond]; + cond.air_density_ = air_density[i_cond++]; } + std::size_t i_species_elem = 0; + for (auto &var : state.variables_.AsVector()) var = concentrations[i_species_elem++]; + + std::size_t i_custom_rate_elem = 0; + for (auto &var : state.custom_rate_parameters_.AsVector()) var = custom_rate_parameters[i_custom_rate_elem++]; solver->CalculateRateConstants(state); auto result = solver->Solve(time_step, state); @@ -467,13 +463,7 @@ namespace musica result.final_time_); i_species_elem = 0; - for (int i_cell{}; i_cell < num_grid_cells_; ++i_cell) - { - for (int i_species{}; i_species < num_species; ++i_species) - { - concentrations[i_species_elem++] = state.variables_[i_cell][i_species]; - } - } + for (auto &var : state.variables_.AsVector()) concentrations[i_species_elem++] = var; DeleteError(error); *error = NoError(); diff --git a/src/test/unit/micm/micm_c_api.cpp b/src/test/unit/micm/micm_c_api.cpp index 8cf1fcdf..bdba8a47 100644 --- a/src/test/unit/micm/micm_c_api.cpp +++ b/src/test/unit/micm/micm_c_api.cpp @@ -268,12 +268,6 @@ void TestSingleGridCell(MICM* micm) DeleteError(&error); } -// Test case for solving system using vector-ordered Rosenbrock solver -TEST_F(MicmCApiTest, SolveUsingVectorOrderedRosenbrock) -{ - TestSingleGridCell(micm); -} - // Test case for solving system using standard-ordered Rosenbrock solver TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock) { @@ -289,21 +283,6 @@ TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock) DeleteError(&error); } -// Test case for solving system using vector-ordered Backward Euler solver -TEST(BackwardEulerVectorOrder, SolveUsingVectordOrderedBackwardEuler) -{ - const char* config_path = "configs/chapman"; - int num_grid_cells = 1; - Error error; - MICM* micm = CreateMicm(config_path, MICMSolver::BackwardEuler, num_grid_cells, &error); - - TestSingleGridCell(micm); - - DeleteMicm(micm, &error); - ASSERT_TRUE(IsSuccess(error)); - DeleteError(&error); -} - // Test case for solving system using standard-ordered Backward Euler solver TEST(BackwardEulerStandardOrder, SolveUsingStandardOrderedBackwardEuler) { @@ -334,14 +313,13 @@ double CalculateArrhenius(const ArrheniusReaction parameters, const double tempe (1.0 + parameters.E_ * pressure); } -// Common test function for solving multiple grid cells -void TestMultipleGridCells(MICM* micm, const size_t num_grid_cells) +// Common test function for solving multiple grid cells with standard-ordered matrices +void TestStandardMultipleGridCells(MICM* micm, const size_t num_grid_cells, const double time_step, const double test_accuracy) { const size_t num_concentrations = 6; const size_t num_user_defined_reaction_rates = 2; constexpr double GAS_CONSTANT = 8.31446261815324; // J mol-1 K-1 - const double time_step = 200.0; // s - + double* temperature = new double[num_grid_cells]; double* pressure = new double[num_grid_cells]; double* air_density = new double[num_grid_cells]; @@ -439,12 +417,131 @@ void TestMultipleGridCells(MICM* micm, const size_t num_grid_cells) double D = initial_D * std::exp(-k1 * time_step); double E = initial_D * (k1 / (k2 - k1)) * (std::exp(-k1 * time_step) - std::exp(-k2 * time_step)); double F = initial_F + initial_D * (1.0 + (k1 * std::exp(-k2 * time_step) - k2 * std::exp(-k1 * time_step)) / (k2 - k1)); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + A_index], A, 5e-3); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + B_index], B, 5e-3); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + C_index], C, 5e-3); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + D_index], D, 5e-3); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + E_index], E, 5e-3); - ASSERT_NEAR(concentrations[i_cell * num_concentrations + F_index], F, 5e-3); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + A_index], A, test_accuracy); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + B_index], B, test_accuracy); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + C_index], C, test_accuracy); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + D_index], D, test_accuracy); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + E_index], E, test_accuracy); + ASSERT_NEAR(concentrations[i_cell * num_concentrations + F_index], F, test_accuracy); + } + delete[] temperature; + delete[] pressure; + delete[] air_density; + delete[] concentrations; + delete[] initial_concentrations; + delete[] user_defined_reaction_rates; +} + +// Common test function for solving multiple grid cells with vectorizable matrices +void TestVectorMultipleGridCells(MICM* micm, const size_t num_grid_cells, const double time_step, const double test_accuracy) +{ + const size_t num_concentrations = 6; + const size_t num_user_defined_reaction_rates = 2; + constexpr double GAS_CONSTANT = 8.31446261815324; // J mol-1 K-1 + + double* temperature = new double[num_grid_cells]; + double* pressure = new double[num_grid_cells]; + double* air_density = new double[num_grid_cells]; + double* concentrations = new double[num_grid_cells * num_concentrations]; + double* initial_concentrations = new double[num_grid_cells * num_concentrations]; + double* user_defined_reaction_rates = new double[num_grid_cells * num_user_defined_reaction_rates]; + + Error error; + + // Get species indices in concentration array + Mappings species_ordering = GetSpeciesOrdering(micm, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(species_ordering.size_, num_concentrations); + std::size_t A_index = FindMappingIndex(species_ordering, "A", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t B_index = FindMappingIndex(species_ordering, "B", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t C_index = FindMappingIndex(species_ordering, "C", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t D_index = FindMappingIndex(species_ordering, "D", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t E_index = FindMappingIndex(species_ordering, "E", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t F_index = FindMappingIndex(species_ordering, "F", &error); + ASSERT_TRUE(IsSuccess(error)); + DeleteMappings(&species_ordering); + + // Get user-defined reaction rates indices in user-defined reaction rates array + Mappings rate_ordering = GetUserDefinedReactionRatesOrdering(micm, &error); + ASSERT_TRUE(IsSuccess(error)); + ASSERT_EQ(rate_ordering.size_, num_user_defined_reaction_rates); + std::size_t R1_index = FindMappingIndex(rate_ordering, "USER.reaction 1", &error); + ASSERT_TRUE(IsSuccess(error)); + std::size_t R2_index = FindMappingIndex(rate_ordering, "USER.reaction 2", &error); + ASSERT_TRUE(IsSuccess(error)); + DeleteMappings(&rate_ordering); + + for (int i = 0; i < num_grid_cells; ++i) + { + temperature[i] = 275.0 + (rand() % 20 - 10); + pressure[i] = 101253.3 + (rand() % 1000 - 500); + air_density[i] = pressure[i] / (GAS_CONSTANT * temperature[i]); + concentrations[i + A_index * num_grid_cells] = 0.75 + (rand() % 10 - 5) * 0.01; + concentrations[i + B_index * num_grid_cells] = 0.0; + concentrations[i + C_index * num_grid_cells] = 0.4 + (rand() % 10 - 5) * 0.01; + concentrations[i + D_index * num_grid_cells] = 0.8 + (rand() % 10 - 5) * 0.01; + concentrations[i + E_index * num_grid_cells] = 0.0; + concentrations[i + F_index * num_grid_cells] = 0.1 + (rand() % 10 - 5) * 0.01; + user_defined_reaction_rates[i + R1_index * num_grid_cells] = 0.001 + (rand() % 10 - 5) * 0.0001; + user_defined_reaction_rates[i + R2_index * num_grid_cells] = 0.002 + (rand() % 10 - 5) * 0.0001; + for (int j = 0; j < num_concentrations; ++j) + { + initial_concentrations[i + j * num_grid_cells] = concentrations[i + j * num_grid_cells]; + } + } + + String solver_state; + SolverResultStats solver_stats; + MicmSolve( + micm, + time_step, + temperature, + pressure, + air_density, + concentrations, + user_defined_reaction_rates, + &solver_state, + &solver_stats, + &error); + ASSERT_TRUE(IsSuccess(error)); + DeleteError(&error); + + // Add assertions to check the solved concentrations + ArrheniusReaction arr1; + arr1.A_ = 0.004; + arr1.C_ = 50.0; + ArrheniusReaction arr2{ 0.012, -2, 75, 50, 1.0e-6 }; + + ASSERT_STREQ(solver_state.value_, "Converged"); + DeleteString(&solver_state); + + for (int i_cell = 0; i_cell < num_grid_cells; ++i_cell) + { + double initial_A = initial_concentrations[i_cell + A_index * num_grid_cells]; + double initial_C = initial_concentrations[i_cell + C_index * num_grid_cells]; + double initial_D = initial_concentrations[i_cell + D_index * num_grid_cells]; + double initial_F = initial_concentrations[i_cell + F_index * num_grid_cells]; + double k1 = user_defined_reaction_rates[i_cell + R1_index * num_grid_cells]; + double k2 = user_defined_reaction_rates[i_cell + R2_index * num_grid_cells]; + double k3 = CalculateArrhenius(arr1, temperature[i_cell], pressure[i_cell]); + double k4 = CalculateArrhenius(arr2, temperature[i_cell], pressure[i_cell]); + double A = initial_A * std::exp(-k3 * time_step); + double B = initial_A * (k3 / (k4 - k3)) * (std::exp(-k3 * time_step) - std::exp(-k4 * time_step)); + double C = initial_C + initial_A * (1.0 + (k3 * std::exp(-k4 * time_step) - k4 * std::exp(-k3 * time_step)) / (k4 - k3)); + double D = initial_D * std::exp(-k1 * time_step); + double E = initial_D * (k1 / (k2 - k1)) * (std::exp(-k1 * time_step) - std::exp(-k2 * time_step)); + double F = initial_F + initial_D * (1.0 + (k1 * std::exp(-k2 * time_step) - k2 * std::exp(-k1 * time_step)) / (k2 - k1)); + ASSERT_NEAR(concentrations[i_cell + A_index * num_grid_cells], A, test_accuracy); + ASSERT_NEAR(concentrations[i_cell + B_index * num_grid_cells], B, test_accuracy); + ASSERT_NEAR(concentrations[i_cell + C_index * num_grid_cells], C, test_accuracy); + ASSERT_NEAR(concentrations[i_cell + D_index * num_grid_cells], D, test_accuracy); + ASSERT_NEAR(concentrations[i_cell + E_index * num_grid_cells], E, test_accuracy); + ASSERT_NEAR(concentrations[i_cell + F_index * num_grid_cells], F, test_accuracy); } delete[] temperature; delete[] pressure; @@ -457,14 +554,16 @@ void TestMultipleGridCells(MICM* micm, const size_t num_grid_cells) // Test case for solving multiple grid cells using vector-ordered Rosenbrock solver TEST_F(MicmCApiTest, SolveMultipleGridCellsUsingVectorOrderedRosenbrock) { - constexpr size_t num_grid_cells = 3; + constexpr size_t num_grid_cells = MICM_VECTOR_MATRIX_SIZE; + constexpr double time_step = 200.0; + constexpr double test_accuracy = 5.0e-3; const char* config_path = "configs/analytical"; Error error; DeleteMicm(micm, &error); ASSERT_TRUE(IsSuccess(error)); micm = CreateMicm(config_path, MICMSolver::Rosenbrock, num_grid_cells, &error); ASSERT_TRUE(IsSuccess(error)); - TestMultipleGridCells(micm, num_grid_cells); + TestVectorMultipleGridCells(micm, num_grid_cells, time_step, test_accuracy); DeleteError(&error); } @@ -472,13 +571,47 @@ TEST_F(MicmCApiTest, SolveMultipleGridCellsUsingVectorOrderedRosenbrock) TEST_F(MicmCApiTest, SolveMultipleGridCellsUsingStandardOrderedRosenbrock) { constexpr size_t num_grid_cells = 3; + constexpr double time_step = 200.0; + constexpr double test_accuracy = 5.0e-3; const char* config_path = "configs/analytical"; Error error; DeleteMicm(micm, &error); ASSERT_TRUE(IsSuccess(error)); micm = CreateMicm(config_path, MICMSolver::RosenbrockStandardOrder, num_grid_cells, &error); ASSERT_TRUE(IsSuccess(error)); - TestMultipleGridCells(micm, num_grid_cells); + TestStandardMultipleGridCells(micm, num_grid_cells, time_step, test_accuracy); + DeleteError(&error); +} + +// Test case for solving multiple grid cells using vector-ordered Backward Euler solver +TEST_F(MicmCApiTest, SolveMultipleGridCellsUsingVectorOrderedBackwardEuler) +{ + constexpr size_t num_grid_cells = MICM_VECTOR_MATRIX_SIZE; + constexpr double time_step = 10.0; + constexpr double test_accuracy = 0.1; + const char* config_path = "configs/analytical"; + Error error; + DeleteMicm(micm, &error); + ASSERT_TRUE(IsSuccess(error)); + micm = CreateMicm(config_path, MICMSolver::BackwardEuler, num_grid_cells, &error); + ASSERT_TRUE(IsSuccess(error)); + TestVectorMultipleGridCells(micm, num_grid_cells, time_step, test_accuracy); + DeleteError(&error); +} + +// Test case for solving multiple grid cells using standard-ordered Backward Euler solver +TEST_F(MicmCApiTest, SolveMultipleGridCellsUsingStandardOrderedBackwardEuler) +{ + constexpr size_t num_grid_cells = 3; + constexpr double time_step = 10.0; + constexpr double test_accuracy = 0.1; + const char* config_path = "configs/analytical"; + Error error; + DeleteMicm(micm, &error); + ASSERT_TRUE(IsSuccess(error)); + micm = CreateMicm(config_path, MICMSolver::BackwardEulerStandardOrder, num_grid_cells, &error); + ASSERT_TRUE(IsSuccess(error)); + TestStandardMultipleGridCells(micm, num_grid_cells, time_step, test_accuracy); DeleteError(&error); } From e4c978662079f1f66f7e2e70206ba40a1bf6426e Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Tue, 5 Nov 2024 10:33:02 -0800 Subject: [PATCH 09/13] Auto-format code using Clang-Format (#241) Co-authored-by: GitHub Actions --- src/micm/micm.cpp | 9 ++++++--- src/test/unit/micm/micm_c_api.cpp | 8 ++++++-- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index 5b46c2d0..d10656b6 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -442,10 +442,12 @@ namespace musica cond.air_density_ = air_density[i_cond++]; } std::size_t i_species_elem = 0; - for (auto &var : state.variables_.AsVector()) var = concentrations[i_species_elem++]; + for (auto &var : state.variables_.AsVector()) + var = concentrations[i_species_elem++]; std::size_t i_custom_rate_elem = 0; - for (auto &var : state.custom_rate_parameters_.AsVector()) var = custom_rate_parameters[i_custom_rate_elem++]; + for (auto &var : state.custom_rate_parameters_.AsVector()) + var = custom_rate_parameters[i_custom_rate_elem++]; solver->CalculateRateConstants(state); auto result = solver->Solve(time_step, state); @@ -463,7 +465,8 @@ namespace musica result.final_time_); i_species_elem = 0; - for (auto &var : state.variables_.AsVector()) concentrations[i_species_elem++] = var; + for (auto &var : state.variables_.AsVector()) + concentrations[i_species_elem++] = var; DeleteError(error); *error = NoError(); diff --git a/src/test/unit/micm/micm_c_api.cpp b/src/test/unit/micm/micm_c_api.cpp index bdba8a47..6e96d6bf 100644 --- a/src/test/unit/micm/micm_c_api.cpp +++ b/src/test/unit/micm/micm_c_api.cpp @@ -314,12 +314,16 @@ double CalculateArrhenius(const ArrheniusReaction parameters, const double tempe } // Common test function for solving multiple grid cells with standard-ordered matrices -void TestStandardMultipleGridCells(MICM* micm, const size_t num_grid_cells, const double time_step, const double test_accuracy) +void TestStandardMultipleGridCells( + MICM* micm, + const size_t num_grid_cells, + const double time_step, + const double test_accuracy) { const size_t num_concentrations = 6; const size_t num_user_defined_reaction_rates = 2; constexpr double GAS_CONSTANT = 8.31446261815324; // J mol-1 K-1 - + double* temperature = new double[num_grid_cells]; double* pressure = new double[num_grid_cells]; double* air_density = new double[num_grid_cells]; From 3f1718e2cd405cac3a74b13fa0169860f271f51e Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Tue, 5 Nov 2024 10:33:26 -0800 Subject: [PATCH 10/13] Auto-format code using Clang-Format (#240) Co-authored-by: GitHub Actions --- python/test/test_analytical.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/test/test_analytical.py b/python/test/test_analytical.py index f5fb7138..334bc492 100644 --- a/python/test/test_analytical.py +++ b/python/test/test_analytical.py @@ -101,7 +101,6 @@ def test_simulation(self): TestSingleGridCell(self, solver, 200.0, 5) - class TestAnalyticalStandardBackwardEuler(unittest.TestCase): def test_simulation(self): solver = musica.create_solver( @@ -355,7 +354,8 @@ def test_simulation(self): "configs/analytical", musica.micmsolver.rosenbrock, 4) - TestVectorMultipleGridCell(self, solver, 4, 200.0, 5) # The number of grid cells must equal the MICM matrix vector dimension + # The number of grid cells must equal the MICM matrix vector dimension + TestVectorMultipleGridCell(self, solver, 4, 200.0, 5) class TestAnalyticalStandardRosenbrockMultipleGridCells(unittest.TestCase): @@ -373,7 +373,8 @@ def test_simulation(self): "configs/analytical", musica.micmsolver.backward_euler, 4) - TestVectorMultipleGridCell(self, solver, 4, 10.0, places=2) # The number of grid cells must equal the MICM matrix vector dimension + # The number of grid cells must equal the MICM matrix vector dimension + TestVectorMultipleGridCell(self, solver, 4, 10.0, places=2) class TestAnalyticalStandardBackwardEulerMultipleGridCells(unittest.TestCase): From 9993aa333be7c8e339675df8f0b037edd04177b6 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Wed, 6 Nov 2024 10:39:55 -0800 Subject: [PATCH 11/13] Add `contiguous` attribute for arrays passed to c functions (#242) * add contiguous attribute for arrays passed to c functions * add tests with array slices and fix intents * revert to previous intel image to avoid gfortran ICE --- docker/Dockerfile.fortran-intel | 2 +- fortran/micm.F90 | 20 +- .../test_micm_api.F90 | 102 ++-- .../test_tuvx_api.F90 | 500 +++++++++--------- fortran/tuvx/grid.F90 | 26 +- fortran/tuvx/profile.F90 | 38 +- fortran/tuvx/radiator.F90 | 36 +- fortran/tuvx/tuvx.F90 | 12 +- fortran/util.F90 | 6 +- 9 files changed, 380 insertions(+), 362 deletions(-) diff --git a/docker/Dockerfile.fortran-intel b/docker/Dockerfile.fortran-intel index a3f3ec94..da120605 100644 --- a/docker/Dockerfile.fortran-intel +++ b/docker/Dockerfile.fortran-intel @@ -1,5 +1,5 @@ # versions and sizes from here: https://hub.docker.com/r/intel/oneapi-hpckit/tags -FROM intel/oneapi-hpckit:latest +FROM intel/oneapi-hpckit:2024.0.1-devel-ubuntu22.04 # Based off of this: https://dgpu-docs.intel.com/driver/installation.html#repository-public-key-used-for-package-and-repository-signing # however those docs (at the time of this writing are incorrect) and this is the correct url diff --git a/fortran/micm.F90 b/fortran/micm.F90 index 92383920..e55ab0d3 100644 --- a/fortran/micm.F90 +++ b/fortran/micm.F90 @@ -252,16 +252,16 @@ subroutine solve_arrays(this, time_step, temperature, pressure, air_density, & 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 - class(micm_t), intent(in) :: this - real(real64), intent(in) :: time_step - 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(:,:) - type(string_t), intent(out) :: solver_state - type(solver_stats_t), intent(out) :: solver_stats - type(error_t), intent(out) :: error + class(micm_t), intent(in) :: this + real(real64), intent(in) :: time_step + real(real64), target, contiguous, intent(in) :: temperature(:) + real(real64), target, contiguous, intent(in) :: pressure(:) + real(real64), target, contiguous, intent(in) :: air_density(:) + real(real64), target, contiguous, intent(inout) :: concentrations(:,:) + real(real64), target, contiguous, 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 diff --git a/fortran/test/fetch_content_integration/test_micm_api.F90 b/fortran/test/fetch_content_integration/test_micm_api.F90 index 7168f853..87d1a909 100644 --- a/fortran/test/fetch_content_integration/test_micm_api.F90 +++ b/fortran/test/fetch_content_integration/test_micm_api.F90 @@ -182,13 +182,18 @@ subroutine test_vector_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test integer, parameter :: NUM_SPECIES = 6 integer, parameter :: NUM_USER_DEFINED_REACTION_RATES = 2 - real(real64), target :: temperature(NUM_GRID_CELLS) - real(real64), target :: pressure(NUM_GRID_CELLS) - real(real64), target :: air_density(NUM_GRID_CELLS) - real(real64), target :: concentrations(NUM_GRID_CELLS,NUM_SPECIES) + ! set up arrays to pass to MICM as slices to ensure contiguous memory is passed to c functions + real(real64), target :: temperature(2,NUM_GRID_CELLS) + real(real64), target :: temperature_c_ptrs(NUM_GRID_CELLS) + real(real64), target :: pressure(2,NUM_GRID_CELLS) + real(real64), target :: pressure_c_ptrs(NUM_GRID_CELLS) + real(real64), target :: air_density(3,NUM_GRID_CELLS) + real(real64), target :: air_density_c_ptrs(NUM_GRID_CELLS) + real(real64), target :: concentrations(4,NUM_GRID_CELLS,NUM_SPECIES) real(real64), target :: concentrations_c_ptrs(NUM_GRID_CELLS,NUM_SPECIES) - real(real64), target :: initial_concentrations(NUM_GRID_CELLS,NUM_SPECIES) - real(real64), target :: user_defined_reaction_rates(NUM_GRID_CELLS,NUM_USER_DEFINED_REACTION_RATES) + real(real64), target :: initial_concentrations(4,NUM_GRID_CELLS,NUM_SPECIES) + real(real64), target :: user_defined_reaction_rates(3,NUM_GRID_CELLS,NUM_USER_DEFINED_REACTION_RATES) + real(real64), target :: user_defined_reaction_rates_c_ptrs(NUM_GRID_CELLS,NUM_USER_DEFINED_REACTION_RATES) type(string_t) :: solver_state type(solver_stats_t) :: solver_stats integer :: solver_type @@ -221,51 +226,62 @@ subroutine test_vector_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test R2_index = micm%user_defined_reaction_rates%index( "USER.reaction 2", error ) ASSERT( error%is_success() ) + temperature(:,:) = 1.0e300_real64 + pressure(:,:) = 1.0e300_real64 + air_density(:,:) = 1.0e300_real64 + concentrations(:,:,:) = 1.0e300_real64 + user_defined_reaction_rates(:,:,:) = 1.0e300_real64 do i_cell = 1, NUM_GRID_CELLS call random_number( temp ) - temperature(i_cell) = 265.0 + temp * 20.0 + temperature(2,i_cell) = 265.0 + temp * 20.0 call random_number( temp ) - pressure(i_cell) = 100753.3 + temp * 1000.0 - air_density(i_cell) = pressure(i_cell) / ( GAS_CONSTANT * temperature(i_cell) ) + pressure(2,i_cell) = 100753.3 + temp * 1000.0 + air_density(2,i_cell) = pressure(2,i_cell) / ( GAS_CONSTANT * temperature(2,i_cell) ) call random_number( temp ) - concentrations(i_cell,A_index) = 0.7 + temp * 0.1 - concentrations(i_cell,B_index) = 0.0 + concentrations(2,i_cell,A_index) = 0.7 + temp * 0.1 + concentrations(2,i_cell,B_index) = 0.0 call random_number( temp ) - concentrations(i_cell,C_index) = 0.35 + temp * 0.1 + concentrations(2,i_cell,C_index) = 0.35 + temp * 0.1 call random_number( temp ) - concentrations(i_cell,D_index) = 0.75 + temp * 0.1 - concentrations(i_cell,E_index) = 0.0 + concentrations(2,i_cell,D_index) = 0.75 + temp * 0.1 + concentrations(2,i_cell,E_index) = 0.0 call random_number( temp ) - concentrations(i_cell,F_index) = 0.05 + temp * 0.1 + concentrations(2,i_cell,F_index) = 0.05 + temp * 0.1 call random_number( temp ) - user_defined_reaction_rates(i_cell,R1_index) = 0.0005 + temp * 0.0001 + user_defined_reaction_rates(2,i_cell,R1_index) = 0.0005 + temp * 0.0001 call random_number( temp ) - user_defined_reaction_rates(i_cell,R2_index) = 0.0015 + temp * 0.0001 + user_defined_reaction_rates(2,i_cell,R2_index) = 0.0015 + temp * 0.0001 end do - initial_concentrations(:,:) = concentrations(:,:) - concentrations_c_ptrs(:,:) = concentrations(:,:) + initial_concentrations(:,:,:) = concentrations(:,:,:) + concentrations_c_ptrs(:,:) = concentrations(2,:,:) + user_defined_reaction_rates_c_ptrs(:,:) = user_defined_reaction_rates(2,:,:) + temperature_c_ptrs(:) = temperature(2,:) + pressure_c_ptrs(:) = pressure(2,:) + air_density_c_ptrs(:) = air_density(2,:) ! solve by passing fortran arrays - call micm%solve(time_step, temperature, pressure, air_density, concentrations, & - user_defined_reaction_rates, solver_state, solver_stats, error) + call micm%solve(time_step, temperature(2,:), pressure(2,:), air_density(2,:), & + concentrations(2,:,:), user_defined_reaction_rates(2,:,:), & + solver_state, solver_stats, error) ASSERT( error%is_success() ) ASSERT_EQ(solver_state%get_char_array(), "Converged") ! solve by passing C pointers - call micm%solve(time_step, c_loc(temperature), c_loc(pressure), c_loc(air_density), & - c_loc(concentrations_c_ptrs), c_loc(user_defined_reaction_rates), & + call micm%solve(time_step, c_loc(temperature_c_ptrs), c_loc(pressure_c_ptrs), & + c_loc(air_density_c_ptrs), c_loc(concentrations_c_ptrs), & + c_loc(user_defined_reaction_rates_c_ptrs), & solver_state, solver_stats, error) ASSERT( error%is_success() ) ASSERT_EQ(solver_state%get_char_array(), "Converged") ! check concentrations do i_cell = 1, NUM_GRID_CELLS - ASSERT_EQ(concentrations(i_cell,A_index), concentrations_c_ptrs(i_cell,A_index)) - ASSERT_EQ(concentrations(i_cell,B_index), concentrations_c_ptrs(i_cell,B_index)) - ASSERT_EQ(concentrations(i_cell,C_index), concentrations_c_ptrs(i_cell,C_index)) - ASSERT_EQ(concentrations(i_cell,D_index), concentrations_c_ptrs(i_cell,D_index)) - ASSERT_EQ(concentrations(i_cell,E_index), concentrations_c_ptrs(i_cell,E_index)) - ASSERT_EQ(concentrations(i_cell,F_index), concentrations_c_ptrs(i_cell,F_index)) + ASSERT_EQ(concentrations(2,i_cell,A_index), concentrations_c_ptrs(i_cell,A_index)) + ASSERT_EQ(concentrations(2,i_cell,B_index), concentrations_c_ptrs(i_cell,B_index)) + ASSERT_EQ(concentrations(2,i_cell,C_index), concentrations_c_ptrs(i_cell,C_index)) + ASSERT_EQ(concentrations(2,i_cell,D_index), concentrations_c_ptrs(i_cell,D_index)) + ASSERT_EQ(concentrations(2,i_cell,E_index), concentrations_c_ptrs(i_cell,E_index)) + ASSERT_EQ(concentrations(2,i_cell,F_index), concentrations_c_ptrs(i_cell,F_index)) end do r1%A_ = 0.004 @@ -277,26 +293,26 @@ subroutine test_vector_multiple_grid_cells(micm, NUM_GRID_CELLS, time_step, test r2%E_ = 1.0e-6 do i_cell = 1, NUM_GRID_CELLS - initial_A = initial_concentrations(i_cell,A_index) - initial_C = initial_concentrations(i_cell,C_index) - initial_D = initial_concentrations(i_cell,D_index) - initial_F = initial_concentrations(i_cell,F_index) - k1 = user_defined_reaction_rates(i_cell,R1_index) - k2 = user_defined_reaction_rates(i_cell,R2_index) - k3 = calculate_arrhenius( r1, temperature(i_cell), pressure(i_cell) ) - k4 = calculate_arrhenius( r2, temperature(i_cell), pressure(i_cell) ) + initial_A = initial_concentrations(2,i_cell,A_index) + initial_C = initial_concentrations(2,i_cell,C_index) + initial_D = initial_concentrations(2,i_cell,D_index) + initial_F = initial_concentrations(2,i_cell,F_index) + k1 = user_defined_reaction_rates(2,i_cell,R1_index) + k2 = user_defined_reaction_rates(2,i_cell,R2_index) + k3 = calculate_arrhenius( r1, temperature(2,i_cell), pressure(2,i_cell) ) + k4 = calculate_arrhenius( r2, temperature(2,i_cell), pressure(2,i_cell) ) A = initial_A * exp( -k3 * time_step ) B = initial_A * (k3 / (k4 - k3)) * (exp(-k3 * time_step) - exp(-k4 * time_step)) C = initial_C + initial_A * (1.0 + (k3 * exp(-k4 * time_step) - k4 * exp(-k3 * time_step)) / (k4 - k3)) D = initial_D * exp( -k1 * time_step ) E = initial_D * (k1 / (k2 - k1)) * (exp(-k1 * time_step) - exp(-k2 * time_step)) F = initial_F + initial_D * (1.0 + (k1 * exp(-k2 * time_step) - k2 * exp(-k1 * time_step)) / (k2 - k1)) - ASSERT_NEAR(concentrations(i_cell,A_index), A, test_accuracy) - ASSERT_NEAR(concentrations(i_cell,B_index), B, test_accuracy) - ASSERT_NEAR(concentrations(i_cell,C_index), C, test_accuracy) - ASSERT_NEAR(concentrations(i_cell,D_index), D, test_accuracy) - ASSERT_NEAR(concentrations(i_cell,E_index), E, test_accuracy) - ASSERT_NEAR(concentrations(i_cell,F_index), F, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,A_index), A, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,B_index), B, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,C_index), C, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,D_index), D, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,E_index), E, test_accuracy) + ASSERT_NEAR(concentrations(2,i_cell,F_index), F, test_accuracy) end do end subroutine test_vector_multiple_grid_cells diff --git a/fortran/test/fetch_content_integration/test_tuvx_api.F90 b/fortran/test/fetch_content_integration/test_tuvx_api.F90 index a8e91432..a5e660c7 100644 --- a/fortran/test/fetch_content_integration/test_tuvx_api.F90 +++ b/fortran/test/fetch_content_integration/test_tuvx_api.F90 @@ -78,36 +78,38 @@ end subroutine test_tuvx_api_invalid_config subroutine test_tuvx_solve() - type(tuvx_t), pointer :: tuvx - type(error_t) :: error - character(len=256) :: config_path - type(grid_map_t), pointer :: grids, grids_from_host - type(grid_t), pointer :: grid, height_grid, wavelength_grid - type(profile_map_t), pointer :: profiles, profiles_from_host - type(profile_t), pointer :: profile, profile_copy - type(radiator_map_t), pointer :: radiators, radiators_from_host - type(radiator_t), pointer :: radiator, radiator_copy - real*8, dimension(5), target :: edges, edge_values, temp_edge - real*8, dimension(4), target :: midpoints, midpoint_values, layer_densities, temp_midpoint - real*8 :: temp_real - integer :: num_vertical_layers, num_wavelength_bins - real*8, dimension(3,2), target :: optical_depths, temp_od - real*8, dimension(3,2), target :: single_scattering_albedos, temp_ssa - real*8, dimension(3,2,1), target :: asymmetry_factors, temp_asym - - edges = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /) - midpoints = (/ 15.0, 25.0, 35.0, 45.0 /) - edge_values = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) - midpoint_values = (/ 15.0, 25.0, 35.0, 45.0 /) - layer_densities = (/ 2.0, 4.0, 1.0, 7.0 /) + type(tuvx_t), pointer :: tuvx + type(error_t) :: error + character(len=256) :: config_path + type(grid_map_t), pointer :: grids, grids_from_host + type(grid_t), pointer :: grid, height_grid, wavelength_grid + type(profile_map_t), pointer :: profiles, profiles_from_host + type(profile_t), pointer :: profile, profile_copy + type(radiator_map_t), pointer :: radiators, radiators_from_host + type(radiator_t), pointer :: radiator, radiator_copy + ! set up arrays with extra dimensions to test whether arrays passed to + ! c functions are contiguous + real*8, dimension(3,5), target :: edges, edge_values, temp_edge + real*8, dimension(2,4), target :: midpoints, midpoint_values, layer_densities, temp_midpoint + real*8 :: temp_real + integer :: num_vertical_layers, num_wavelength_bins + real*8, dimension(4,3,2), target :: optical_depths, temp_od + real*8, dimension(3,3,2), target :: single_scattering_albedos, temp_ssa + real*8, dimension(2,3,2,1), target :: asymmetry_factors, temp_asym + + edges(2,:) = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /) + midpoints(2,:) = (/ 15.0, 25.0, 35.0, 45.0 /) + edge_values(2,:) = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) + midpoint_values(2,:) = (/ 15.0, 25.0, 35.0, 45.0 /) + layer_densities(2,:) = (/ 2.0, 4.0, 1.0, 7.0 /) num_vertical_layers = 3 num_wavelength_bins = 2 - optical_depths(:,1) = (/ 30.0, 20.0, 10.0 /) - optical_depths(:,2) = (/ 70.0, 80.0, 90.0 /) - single_scattering_albedos(:,1) = (/ 300.0, 200.0, 100.0 /) - single_scattering_albedos(:,2) = (/ 700.0, 800.0, 900.0 /) - asymmetry_factors(:,1,1) = (/ 3.0, 2.0, 1.0 /) - asymmetry_factors(:,2,1) = (/ 7.0, 8.0, 9.0 /) + optical_depths(2,:,1) = (/ 30.0, 20.0, 10.0 /) + optical_depths(2,:,2) = (/ 70.0, 80.0, 90.0 /) + single_scattering_albedos(2,:,1) = (/ 300.0, 200.0, 100.0 /) + single_scattering_albedos(2,:,2) = (/ 700.0, 800.0, 900.0 /) + asymmetry_factors(2,:,1,1) = (/ 3.0, 2.0, 1.0 /) + asymmetry_factors(2,:,2,1) = (/ 7.0, 8.0, 9.0 /) config_path = "examples/ts1_tsmlt.json" @@ -137,132 +139,132 @@ subroutine test_tuvx_solve() ASSERT_EQ( grid%number_of_sections( error ), 4 ) ASSERT( error%is_success() ) - call grid%set_edges( edges, error ) + call grid%set_edges( edges(2,:), error ) ASSERT( error%is_success() ) - call grid%get_edges( temp_edge, error ) + call grid%get_edges( temp_edge(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_edge(1), 1.0 ) - ASSERT_EQ( temp_edge(2), 2.0 ) - ASSERT_EQ( temp_edge(3), 3.0 ) - ASSERT_EQ( temp_edge(4), 4.0 ) - ASSERT_EQ( temp_edge(5), 5.0 ) + ASSERT_EQ( temp_edge(2,1), 1.0 ) + ASSERT_EQ( temp_edge(2,2), 2.0 ) + ASSERT_EQ( temp_edge(2,3), 3.0 ) + ASSERT_EQ( temp_edge(2,4), 4.0 ) + ASSERT_EQ( temp_edge(2,5), 5.0 ) - edges = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) + edges(2,:) = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) - call grid%set_edges( edges, error ) + call grid%set_edges( edges(2,:), error ) ASSERT( error%is_success() ) - call grid%set_midpoints( midpoints, error ) + call grid%set_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - call grid%get_edges( temp_edge, error ) + call grid%get_edges( temp_edge(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_edge(1), 10.0 ) - ASSERT_EQ( temp_edge(2), 20.0 ) - ASSERT_EQ( temp_edge(3), 30.0 ) - ASSERT_EQ( temp_edge(4), 40.0 ) - ASSERT_EQ( temp_edge(5), 50.0 ) + ASSERT_EQ( temp_edge(2,1), 10.0 ) + ASSERT_EQ( temp_edge(2,2), 20.0 ) + ASSERT_EQ( temp_edge(2,3), 30.0 ) + ASSERT_EQ( temp_edge(2,4), 40.0 ) + ASSERT_EQ( temp_edge(2,5), 50.0 ) - call grid%get_midpoints( temp_midpoint, error ) + call grid%get_midpoints( temp_midpoint(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_midpoint(1), 15.0 ) - ASSERT_EQ( temp_midpoint(2), 25.0 ) - ASSERT_EQ( temp_midpoint(3), 35.0 ) - ASSERT_EQ( temp_midpoint(4), 45.0 ) + ASSERT_EQ( temp_midpoint(2,1), 15.0 ) + ASSERT_EQ( temp_midpoint(2,2), 25.0 ) + ASSERT_EQ( temp_midpoint(2,3), 35.0 ) + ASSERT_EQ( temp_midpoint(2,4), 45.0 ) call grids%add( grid, error ) - edges(:) = 0.0 - midpoints(:) = 0.0 + edges(2,:) = 0.0 + midpoints(2,:) = 0.0 - call grid%get_edges( edges, error ) + call grid%get_edges( edges(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( edges(1), 10.0 ) - ASSERT_EQ( edges(2), 20.0 ) - ASSERT_EQ( edges(3), 30.0 ) - ASSERT_EQ( edges(4), 40.0 ) - ASSERT_EQ( edges(5), 50.0 ) + ASSERT_EQ( edges(2,1), 10.0 ) + ASSERT_EQ( edges(2,2), 20.0 ) + ASSERT_EQ( edges(2,3), 30.0 ) + ASSERT_EQ( edges(2,4), 40.0 ) + ASSERT_EQ( edges(2,5), 50.0 ) - call grid%get_midpoints( midpoints, error ) + call grid%get_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( midpoints(1), 15.0 ) - ASSERT_EQ( midpoints(2), 25.0 ) - ASSERT_EQ( midpoints(3), 35.0 ) - ASSERT_EQ( midpoints(4), 45.0 ) + ASSERT_EQ( midpoints(2,1), 15.0 ) + ASSERT_EQ( midpoints(2,2), 25.0 ) + ASSERT_EQ( midpoints(2,3), 35.0 ) + ASSERT_EQ( midpoints(2,4), 45.0 ) - edges = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /) - midpoints = (/ 1.5, 2.5, 3.5, 4.5 /) + edges(2,:) = (/ 1.0, 2.0, 3.0, 4.0, 5.0 /) + midpoints(2,:) = (/ 1.5, 2.5, 3.5, 4.5 /) - call grid%set_edges( edges, error ) + call grid%set_edges( edges(2,:), error ) ASSERT( error%is_success() ) - call grid%set_midpoints( midpoints, error ) + call grid%set_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - edges(:) = 0.0 - midpoints(:) = 0.0 + edges(2,:) = 0.0 + midpoints(2,:) = 0.0 - call grid%get_edges( edges, error ) + call grid%get_edges( edges(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( edges(1), 1.0 ) - ASSERT_EQ( edges(2), 2.0 ) - ASSERT_EQ( edges(3), 3.0 ) - ASSERT_EQ( edges(4), 4.0 ) - ASSERT_EQ( edges(5), 5.0 ) + ASSERT_EQ( edges(2,1), 1.0 ) + ASSERT_EQ( edges(2,2), 2.0 ) + ASSERT_EQ( edges(2,3), 3.0 ) + ASSERT_EQ( edges(2,4), 4.0 ) + ASSERT_EQ( edges(2,5), 5.0 ) - call grid%get_midpoints( midpoints, error ) + call grid%get_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( midpoints(1), 1.5 ) - ASSERT_EQ( midpoints(2), 2.5 ) - ASSERT_EQ( midpoints(3), 3.5 ) - ASSERT_EQ( midpoints(4), 4.5 ) + ASSERT_EQ( midpoints(2,1), 1.5 ) + ASSERT_EQ( midpoints(2,2), 2.5 ) + ASSERT_EQ( midpoints(2,3), 3.5 ) + ASSERT_EQ( midpoints(2,4), 4.5 ) deallocate( grid ) grid => grids%get( "foo", "bars", error ) ASSERT( error%is_success() ) - edges(:) = 0.0 - midpoints(:) = 0.0 + edges(2,:) = 0.0 + midpoints(2,:) = 0.0 - call grid%get_edges( edges, error ) + call grid%get_edges( edges(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( edges(1), 1.0 ) - ASSERT_EQ( edges(2), 2.0 ) - ASSERT_EQ( edges(3), 3.0 ) - ASSERT_EQ( edges(4), 4.0 ) - ASSERT_EQ( edges(5), 5.0 ) + ASSERT_EQ( edges(2,1), 1.0 ) + ASSERT_EQ( edges(2,2), 2.0 ) + ASSERT_EQ( edges(2,3), 3.0 ) + ASSERT_EQ( edges(2,4), 4.0 ) + ASSERT_EQ( edges(2,5), 5.0 ) - call grid%get_midpoints( midpoints, error ) + call grid%get_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( midpoints(1), 1.5 ) - ASSERT_EQ( midpoints(2), 2.5 ) - ASSERT_EQ( midpoints(3), 3.5 ) - ASSERT_EQ( midpoints(4), 4.5 ) + ASSERT_EQ( midpoints(2,1), 1.5 ) + ASSERT_EQ( midpoints(2,2), 2.5 ) + ASSERT_EQ( midpoints(2,3), 3.5 ) + ASSERT_EQ( midpoints(2,4), 4.5 ) - edges = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) - midpoints = (/ 15.0, 25.0, 35.0, 45.0 /) + edges(2,:) = (/ 10.0, 20.0, 30.0, 40.0, 50.0 /) + midpoints(2,:) = (/ 15.0, 25.0, 35.0, 45.0 /) - call grid%set_edges( edges, error ) - call grid%set_midpoints( midpoints, error ) + call grid%set_edges( edges(2,:), error ) + call grid%set_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - edges(:) = 0.0 - midpoints(:) = 0.0 + edges(2,:) = 0.0 + midpoints(2,:) = 0.0 - call grid%get_edges( edges, error ) + call grid%get_edges( edges(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( edges(1), 10.0 ) - ASSERT_EQ( edges(2), 20.0 ) - ASSERT_EQ( edges(3), 30.0 ) - ASSERT_EQ( edges(4), 40.0 ) - ASSERT_EQ( edges(5), 50.0 ) + ASSERT_EQ( edges(2,1), 10.0 ) + ASSERT_EQ( edges(2,2), 20.0 ) + ASSERT_EQ( edges(2,3), 30.0 ) + ASSERT_EQ( edges(2,4), 40.0 ) + ASSERT_EQ( edges(2,5), 50.0 ) - call grid%get_midpoints( midpoints, error ) + call grid%get_midpoints( midpoints(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( midpoints(1), 15.0 ) - ASSERT_EQ( midpoints(2), 25.0 ) - ASSERT_EQ( midpoints(3), 35.0 ) - ASSERT_EQ( midpoints(4), 45.0 ) + ASSERT_EQ( midpoints(2,1), 15.0 ) + ASSERT_EQ( midpoints(2,2), 25.0 ) + ASSERT_EQ( midpoints(2,3), 35.0 ) + ASSERT_EQ( midpoints(2,4), 45.0 ) profiles => tuvx%get_profiles( error ) ASSERT( error%is_success() ) @@ -278,36 +280,36 @@ subroutine test_tuvx_solve() profile => profile_t( "baz", "qux", grid, error ) ASSERT( error%is_success() ) - call profile%set_edge_values( edge_values, error ) + call profile%set_edge_values( edge_values(2,:), error ) ASSERT( error%is_success() ) - call profile%get_edge_values( temp_edge, error ) + call profile%get_edge_values( temp_edge(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_edge(1), 10.0 ) - ASSERT_EQ( temp_edge(2), 20.0 ) - ASSERT_EQ( temp_edge(3), 30.0 ) - ASSERT_EQ( temp_edge(4), 40.0 ) - ASSERT_EQ( temp_edge(5), 50.0 ) + ASSERT_EQ( temp_edge(2,1), 10.0 ) + ASSERT_EQ( temp_edge(2,2), 20.0 ) + ASSERT_EQ( temp_edge(2,3), 30.0 ) + ASSERT_EQ( temp_edge(2,4), 40.0 ) + ASSERT_EQ( temp_edge(2,5), 50.0 ) - call profile%set_midpoint_values( midpoint_values, error ) + call profile%set_midpoint_values( midpoint_values(2,:), error ) ASSERT( error%is_success() ) - call profile%get_midpoint_values( temp_midpoint, error ) + call profile%get_midpoint_values( temp_midpoint(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_midpoint(1), 15.0 ) - ASSERT_EQ( temp_midpoint(2), 25.0 ) - ASSERT_EQ( temp_midpoint(3), 35.0 ) - ASSERT_EQ( temp_midpoint(4), 45.0 ) + ASSERT_EQ( temp_midpoint(2,1), 15.0 ) + ASSERT_EQ( temp_midpoint(2,2), 25.0 ) + ASSERT_EQ( temp_midpoint(2,3), 35.0 ) + ASSERT_EQ( temp_midpoint(2,4), 45.0 ) - call profile%set_layer_densities( layer_densities, error ) + call profile%set_layer_densities( layer_densities(2,:), error ) ASSERT( error%is_success() ) - call profile%get_layer_densities( temp_midpoint, error ) + call profile%get_layer_densities( temp_midpoint(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_midpoint(1), 2.0 ) - ASSERT_EQ( temp_midpoint(2), 4.0 ) - ASSERT_EQ( temp_midpoint(3), 1.0 ) - ASSERT_EQ( temp_midpoint(4), 7.0 ) + ASSERT_EQ( temp_midpoint(2,1), 2.0 ) + ASSERT_EQ( temp_midpoint(2,2), 4.0 ) + ASSERT_EQ( temp_midpoint(2,3), 1.0 ) + ASSERT_EQ( temp_midpoint(2,4), 7.0 ) call profile%set_exo_layer_density( 1.0d0, error ) ASSERT( error%is_success() ) @@ -316,12 +318,12 @@ subroutine test_tuvx_solve() ASSERT( error%is_success() ) ASSERT_EQ( temp_real, 1.0 ) - call profile%get_layer_densities( temp_midpoint, error ) + call profile%get_layer_densities( temp_midpoint(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_midpoint(1), 2.0 ) - ASSERT_EQ( temp_midpoint(2), 4.0 ) - ASSERT_EQ( temp_midpoint(3), 1.0 ) - ASSERT_EQ( temp_midpoint(4), 7.0 + 1.0 ) + ASSERT_EQ( temp_midpoint(2,1), 2.0 ) + ASSERT_EQ( temp_midpoint(2,2), 4.0 ) + ASSERT_EQ( temp_midpoint(2,3), 1.0 ) + ASSERT_EQ( temp_midpoint(2,4), 7.0 + 1.0 ) call profile%calculate_exo_layer_density( 10.0d0, error ) ASSERT( error%is_success() ) @@ -331,34 +333,34 @@ subroutine test_tuvx_solve() ! Revisit this after non-SI units are converted in the TUV-x internal functions ASSERT_EQ( temp_real, 10.0 * 7.0 * 100.0 ) - call profile%get_layer_densities( temp_midpoint, error ) + call profile%get_layer_densities( temp_midpoint(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_midpoint(1), 2.0 ) - ASSERT_EQ( temp_midpoint(2), 4.0 ) - ASSERT_EQ( temp_midpoint(3), 1.0 ) - ASSERT_EQ( temp_midpoint(4), 7.0 + 10.0 * 7.0 * 100.0 ) + ASSERT_EQ( temp_midpoint(2,1), 2.0 ) + ASSERT_EQ( temp_midpoint(2,2), 4.0 ) + ASSERT_EQ( temp_midpoint(2,3), 1.0 ) + ASSERT_EQ( temp_midpoint(2,4), 7.0 + 10.0 * 7.0 * 100.0 ) call profiles%add( profile, error ) profile_copy => profiles%get( "baz", "qux", error ) - call profile_copy%get_edge_values( temp_edge, error ) + call profile_copy%get_edge_values( temp_edge(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_edge(1), 10.0 ) - ASSERT_EQ( temp_edge(2), 20.0 ) - ASSERT_EQ( temp_edge(3), 30.0 ) - ASSERT_EQ( temp_edge(4), 40.0 ) - ASSERT_EQ( temp_edge(5), 50.0 ) + ASSERT_EQ( temp_edge(2,1), 10.0 ) + ASSERT_EQ( temp_edge(2,2), 20.0 ) + ASSERT_EQ( temp_edge(2,3), 30.0 ) + ASSERT_EQ( temp_edge(2,4), 40.0 ) + ASSERT_EQ( temp_edge(2,5), 50.0 ) - edge_values = (/ 32.0, 34.0, 36.0, 38.0, 40.0 /) - call profile_copy%set_edge_values( edge_values, error ) + edge_values(2,:) = (/ 32.0, 34.0, 36.0, 38.0, 40.0 /) + call profile_copy%set_edge_values( edge_values(2,:), error ) - call profile%get_edge_values( temp_edge, error ) + call profile%get_edge_values( temp_edge(2,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_edge(1), 32.0 ) - ASSERT_EQ( temp_edge(2), 34.0 ) - ASSERT_EQ( temp_edge(3), 36.0 ) - ASSERT_EQ( temp_edge(4), 38.0 ) - ASSERT_EQ( temp_edge(5), 40.0 ) + ASSERT_EQ( temp_edge(2,1), 32.0 ) + ASSERT_EQ( temp_edge(2,2), 34.0 ) + ASSERT_EQ( temp_edge(2,3), 36.0 ) + ASSERT_EQ( temp_edge(2,4), 38.0 ) + ASSERT_EQ( temp_edge(2,5), 40.0 ) radiators => tuvx%get_radiators( error ) ASSERT( error%is_success() ) @@ -376,118 +378,118 @@ subroutine test_tuvx_solve() radiator => radiator_t( "foo_radiator", height_grid, wavelength_grid, error ) ASSERT( error%is_success() ) - call radiator%set_optical_depths( optical_depths, error ) + call radiator%set_optical_depths( optical_depths(2,:,:), error ) ASSERT( error%is_success() ) - call radiator%get_optical_depths( temp_od, error ) + call radiator%get_optical_depths( temp_od(2,:,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_od(1,1), 30.0 ) - ASSERT_EQ( temp_od(2,1), 20.0 ) - ASSERT_EQ( temp_od(3,1), 10.0 ) - ASSERT_EQ( temp_od(1,2), 70.0 ) - ASSERT_EQ( temp_od(2,2), 80.0 ) - ASSERT_EQ( temp_od(3,2), 90.0 ) + ASSERT_EQ( temp_od(2,1,1), 30.0 ) + ASSERT_EQ( temp_od(2,2,1), 20.0 ) + ASSERT_EQ( temp_od(2,3,1), 10.0 ) + ASSERT_EQ( temp_od(2,1,2), 70.0 ) + ASSERT_EQ( temp_od(2,2,2), 80.0 ) + ASSERT_EQ( temp_od(2,3,2), 90.0 ) - call radiator%set_single_scattering_albedos( single_scattering_albedos, error ) + call radiator%set_single_scattering_albedos( single_scattering_albedos(2,:,:), error ) ASSERT( error%is_success() ) - call radiator%get_single_scattering_albedos( temp_ssa, error ) + call radiator%get_single_scattering_albedos( temp_ssa(2,:,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_ssa(1,1), 300.0 ) - ASSERT_EQ( temp_ssa(2,1), 200.0 ) - ASSERT_EQ( temp_ssa(3,1), 100.0 ) - ASSERT_EQ( temp_ssa(1,2), 700.0 ) - ASSERT_EQ( temp_ssa(2,2), 800.0 ) - ASSERT_EQ( temp_ssa(3,2), 900.0 ) + ASSERT_EQ( temp_ssa(2,1,1), 300.0 ) + ASSERT_EQ( temp_ssa(2,2,1), 200.0 ) + ASSERT_EQ( temp_ssa(2,3,1), 100.0 ) + ASSERT_EQ( temp_ssa(2,1,2), 700.0 ) + ASSERT_EQ( temp_ssa(2,2,2), 800.0 ) + ASSERT_EQ( temp_ssa(2,3,2), 900.0 ) - call radiator%set_asymmetry_factors( asymmetry_factors, error ) + call radiator%set_asymmetry_factors( asymmetry_factors(2,:,:,:), error ) ASSERT( error%is_success() ) - call radiator%get_asymmetry_factors( temp_asym, error ) + call radiator%get_asymmetry_factors( temp_asym(2,:,:,:), error ) ASSERT( error%is_success() ) - ASSERT_EQ( temp_asym(1,1,1), 3.0 ) - ASSERT_EQ( temp_asym(2,1,1), 2.0 ) - ASSERT_EQ( temp_asym(3,1,1), 1.0 ) - ASSERT_EQ( temp_asym(1,2,1), 7.0 ) - ASSERT_EQ( temp_asym(2,2,1), 8.0 ) - ASSERT_EQ( temp_asym(3,2,1), 9.0 ) + ASSERT_EQ( temp_asym(2,1,1,1), 3.0 ) + ASSERT_EQ( temp_asym(2,2,1,1), 2.0 ) + ASSERT_EQ( temp_asym(2,3,1,1), 1.0 ) + ASSERT_EQ( temp_asym(2,1,2,1), 7.0 ) + ASSERT_EQ( temp_asym(2,2,2,1), 8.0 ) + ASSERT_EQ( temp_asym(2,3,2,1), 9.0 ) ! call radiators%add( radiator, error ) radiator_copy => radiators%get( "foo_radiator", error ) - optical_depths(:,:) = 0.0 - single_scattering_albedos(:,:) = 0.0 - asymmetry_factors(:,:,:) = 0.0 - - call radiator_copy%get_optical_depths( optical_depths, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( optical_depths(1,1), 30.0 ) - ASSERT_EQ( optical_depths(2,1), 20.0 ) - ASSERT_EQ( optical_depths(3,1), 10.0 ) - ASSERT_EQ( optical_depths(1,2), 70.0 ) - ASSERT_EQ( optical_depths(2,2), 80.0 ) - ASSERT_EQ( optical_depths(3,2), 90.0 ) - - call radiator_copy%get_single_scattering_albedos( single_scattering_albedos, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( single_scattering_albedos(1,1), 300.0 ) - ASSERT_EQ( single_scattering_albedos(2,1), 200.0 ) - ASSERT_EQ( single_scattering_albedos(3,1), 100.0 ) - ASSERT_EQ( single_scattering_albedos(1,2), 700.0 ) - ASSERT_EQ( single_scattering_albedos(2,2), 800.0 ) - ASSERT_EQ( single_scattering_albedos(3,2), 900.0 ) - - call radiator_copy%get_asymmetry_factors( asymmetry_factors, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( asymmetry_factors(1,1,1), 3.0 ) - ASSERT_EQ( asymmetry_factors(2,1,1), 2.0 ) - ASSERT_EQ( asymmetry_factors(3,1,1), 1.0 ) - ASSERT_EQ( asymmetry_factors(1,2,1), 7.0 ) - ASSERT_EQ( asymmetry_factors(2,2,1), 8.0 ) - ASSERT_EQ( asymmetry_factors(3,2,1), 9.0 ) - - optical_depths(:,1) = (/ 90.0, 80.0, 70.0 /) - optical_depths(:,2) = (/ 75.0, 85.0, 95.0 /) - single_scattering_albedos(:,1) = (/ 900.0, 800.0, 700.0 /) - single_scattering_albedos(:,2) = (/ 750.0, 850.0, 950.0 /) - asymmetry_factors(:,1,1) = (/ 9.0, 8.0, 7.0 /) - asymmetry_factors(:,2,1) = (/ 5.0, 4.0, 3.0 /) - - call radiator_copy%set_optical_depths( optical_depths, error ) - call radiator_copy%set_single_scattering_albedos( single_scattering_albedos, error ) - call radiator_copy%set_asymmetry_factors( asymmetry_factors, error ) - ASSERT( error%is_success() ) - - optical_depths(:,:) = 0.0 - single_scattering_albedos(:,:) = 0.0 - asymmetry_factors(:,:,:) = 0.0 - - call radiator%get_optical_depths( optical_depths, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( optical_depths(1,1), 90.0 ) - ASSERT_EQ( optical_depths(2,1), 80.0 ) - ASSERT_EQ( optical_depths(3,1), 70.0 ) - ASSERT_EQ( optical_depths(1,2), 75.0 ) - ASSERT_EQ( optical_depths(2,2), 85.0 ) - ASSERT_EQ( optical_depths(3,2), 95.0 ) - - call radiator%get_single_scattering_albedos( single_scattering_albedos, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( single_scattering_albedos(1,1), 900.0 ) - ASSERT_EQ( single_scattering_albedos(2,1), 800.0 ) - ASSERT_EQ( single_scattering_albedos(3,1), 700.0 ) - ASSERT_EQ( single_scattering_albedos(1,2), 750.0 ) - ASSERT_EQ( single_scattering_albedos(2,2), 850.0 ) - ASSERT_EQ( single_scattering_albedos(3,2), 950.0 ) - - call radiator%get_asymmetry_factors( asymmetry_factors, error ) - ASSERT( error%is_success() ) - ASSERT_EQ( asymmetry_factors(1,1,1), 9.0 ) - ASSERT_EQ( asymmetry_factors(2,1,1), 8.0 ) - ASSERT_EQ( asymmetry_factors(3,1,1), 7.0 ) - ASSERT_EQ( asymmetry_factors(1,2,1), 5.0 ) - ASSERT_EQ( asymmetry_factors(2,2,1), 4.0 ) - ASSERT_EQ( asymmetry_factors(3,2,1), 3.0 ) + optical_depths(2,:,:) = 0.0 + single_scattering_albedos(2,:,:) = 0.0 + asymmetry_factors(2,:,:,:) = 0.0 + + call radiator_copy%get_optical_depths( optical_depths(2,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( optical_depths(2,1,1), 30.0 ) + ASSERT_EQ( optical_depths(2,2,1), 20.0 ) + ASSERT_EQ( optical_depths(2,3,1), 10.0 ) + ASSERT_EQ( optical_depths(2,1,2), 70.0 ) + ASSERT_EQ( optical_depths(2,2,2), 80.0 ) + ASSERT_EQ( optical_depths(2,3,2), 90.0 ) + + call radiator_copy%get_single_scattering_albedos( single_scattering_albedos(2,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( single_scattering_albedos(2,1,1), 300.0 ) + ASSERT_EQ( single_scattering_albedos(2,2,1), 200.0 ) + ASSERT_EQ( single_scattering_albedos(2,3,1), 100.0 ) + ASSERT_EQ( single_scattering_albedos(2,1,2), 700.0 ) + ASSERT_EQ( single_scattering_albedos(2,2,2), 800.0 ) + ASSERT_EQ( single_scattering_albedos(2,3,2), 900.0 ) + + call radiator_copy%get_asymmetry_factors( asymmetry_factors(2,:,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( asymmetry_factors(2,1,1,1), 3.0 ) + ASSERT_EQ( asymmetry_factors(2,2,1,1), 2.0 ) + ASSERT_EQ( asymmetry_factors(2,3,1,1), 1.0 ) + ASSERT_EQ( asymmetry_factors(2,1,2,1), 7.0 ) + ASSERT_EQ( asymmetry_factors(2,2,2,1), 8.0 ) + ASSERT_EQ( asymmetry_factors(2,3,2,1), 9.0 ) + + optical_depths(2,:,1) = (/ 90.0, 80.0, 70.0 /) + optical_depths(2,:,2) = (/ 75.0, 85.0, 95.0 /) + single_scattering_albedos(2,:,1) = (/ 900.0, 800.0, 700.0 /) + single_scattering_albedos(2,:,2) = (/ 750.0, 850.0, 950.0 /) + asymmetry_factors(2,:,1,1) = (/ 9.0, 8.0, 7.0 /) + asymmetry_factors(2,:,2,1) = (/ 5.0, 4.0, 3.0 /) + + call radiator_copy%set_optical_depths( optical_depths(2,:,:), error ) + call radiator_copy%set_single_scattering_albedos( single_scattering_albedos(2,:,:), error ) + call radiator_copy%set_asymmetry_factors( asymmetry_factors(2,:,:,:), error ) + ASSERT( error%is_success() ) + + optical_depths(:,:,:) = 0.0 + single_scattering_albedos(:,:,:) = 0.0 + asymmetry_factors(:,:,:,:) = 0.0 + + call radiator%get_optical_depths( optical_depths(2,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( optical_depths(2,1,1), 90.0 ) + ASSERT_EQ( optical_depths(2,2,1), 80.0 ) + ASSERT_EQ( optical_depths(2,3,1), 70.0 ) + ASSERT_EQ( optical_depths(2,1,2), 75.0 ) + ASSERT_EQ( optical_depths(2,2,2), 85.0 ) + ASSERT_EQ( optical_depths(2,3,2), 95.0 ) + + call radiator%get_single_scattering_albedos( single_scattering_albedos(2,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( single_scattering_albedos(2,1,1), 900.0 ) + ASSERT_EQ( single_scattering_albedos(2,2,1), 800.0 ) + ASSERT_EQ( single_scattering_albedos(2,3,1), 700.0 ) + ASSERT_EQ( single_scattering_albedos(2,1,2), 750.0 ) + ASSERT_EQ( single_scattering_albedos(2,2,2), 850.0 ) + ASSERT_EQ( single_scattering_albedos(2,3,2), 950.0 ) + + call radiator%get_asymmetry_factors( asymmetry_factors(2,:,:,:), error ) + ASSERT( error%is_success() ) + ASSERT_EQ( asymmetry_factors(2,1,1,1), 9.0 ) + ASSERT_EQ( asymmetry_factors(2,2,1,1), 8.0 ) + ASSERT_EQ( asymmetry_factors(2,3,1,1), 7.0 ) + ASSERT_EQ( asymmetry_factors(2,1,2,1), 5.0 ) + ASSERT_EQ( asymmetry_factors(2,2,2,1), 4.0 ) + ASSERT_EQ( asymmetry_factors(2,3,2,1), 3.0 ) deallocate( grid ) deallocate( grids ) diff --git a/fortran/tuvx/grid.F90 b/fortran/tuvx/grid.F90 index db9b9563..2311ff08 100644 --- a/fortran/tuvx/grid.F90 +++ b/fortran/tuvx/grid.F90 @@ -157,7 +157,7 @@ integer function number_of_sections(this, error) result( n_sections ) use musica_util, only: error_t, error_t_c ! Arguments - class(grid_t), intent(inout) :: this + class(grid_t), intent(in) :: this type(error_t), intent(inout) :: error ! Local variables @@ -175,9 +175,9 @@ subroutine set_edges(this, edges, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(grid_t), intent(inout) :: this - real(dk), target, dimension(:), intent(in) :: edges - type(error_t), intent(inout) :: error + class(grid_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: edges(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -197,9 +197,9 @@ subroutine get_edges(this, edges, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(grid_t), intent(inout) :: this - real(dk), target, dimension(:), intent(inout) :: edges - type(error_t), intent(inout) :: error + class(grid_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: edges(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -219,9 +219,9 @@ subroutine set_midpoints(this, midpoints, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(grid_t), intent(inout) :: this - real(dk), target, dimension(:), intent(in) :: midpoints - type(error_t), intent(inout) :: error + class(grid_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: midpoints(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -241,9 +241,9 @@ subroutine get_midpoints(this, midpoints, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(grid_t), intent(inout) :: this - real(dk), target, dimension(:), intent(inout) :: midpoints - type(error_t), intent(inout) :: error + class(grid_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: midpoints(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c diff --git a/fortran/tuvx/profile.F90 b/fortran/tuvx/profile.F90 index 757a08fc..15bf6bb3 100644 --- a/fortran/tuvx/profile.F90 +++ b/fortran/tuvx/profile.F90 @@ -205,9 +205,9 @@ subroutine set_edge_values(this, edge_values, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(in) :: edge_values - type(error_t), intent(inout) :: error + class(profile_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: edge_values(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -228,9 +228,9 @@ subroutine get_edge_values(this, edge_values, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(inout) :: edge_values - type(error_t), intent(inout) :: error + class(profile_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: edge_values(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -251,9 +251,9 @@ subroutine set_midpoint_values(this, midpoint_values, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(in) :: midpoint_values - type(error_t), intent(inout) :: error + class(profile_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: midpoint_values(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -274,9 +274,9 @@ subroutine get_midpoint_values(this, midpoint_values, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(inout) :: midpoint_values - type(error_t), intent(inout) :: error + class(profile_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: midpoint_values(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -297,9 +297,9 @@ subroutine set_layer_densities(this, layer_densities, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(in) :: layer_densities - type(error_t), intent(inout) :: error + class(profile_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: layer_densities(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -320,9 +320,9 @@ subroutine get_layer_densities(this, layer_densities, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this - real(dk), target, dimension(:), intent(inout) :: layer_densities - type(error_t), intent(inout) :: error + class(profile_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: layer_densities(:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -383,7 +383,7 @@ function get_exo_layer_density(this, error) result(exo_layer_density) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(profile_t), intent(inout) :: this + class(profile_t), intent(in) :: this type(error_t), intent(inout) :: error ! Return value diff --git a/fortran/tuvx/radiator.F90 b/fortran/tuvx/radiator.F90 index 60a79ee6..9ebd0b27 100644 --- a/fortran/tuvx/radiator.F90 +++ b/fortran/tuvx/radiator.F90 @@ -182,9 +182,9 @@ subroutine set_optical_depths(this, optical_depths, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:), intent(in) :: optical_depths - type(error_t), intent(inout) :: error + class(radiator_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: optical_depths(:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -207,9 +207,9 @@ subroutine get_optical_depths(this, optical_depths, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:), intent(in) :: optical_depths - type(error_t), intent(inout) :: error + class(radiator_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: optical_depths(:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -233,9 +233,9 @@ subroutine set_single_scattering_albedos(this, single_scattering_albedos, & use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:), intent(in) :: single_scattering_albedos - type(error_t), intent(inout) :: error + class(radiator_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: single_scattering_albedos(:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -260,9 +260,9 @@ subroutine get_single_scattering_albedos(this, single_scattering_albedos, & use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:), intent(in) :: single_scattering_albedos - type(error_t), intent(inout) :: error + class(radiator_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: single_scattering_albedos(:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -286,9 +286,9 @@ subroutine set_asymmetry_factors(this, asymmetry_factors, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:,:), intent(in) :: asymmetry_factors - type(error_t), intent(inout) :: error + class(radiator_t), intent(inout) :: this + real(dk), target, contiguous, intent(in) :: asymmetry_factors(:,:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c @@ -313,9 +313,9 @@ subroutine get_asymmetry_factors(this, asymmetry_factors, error) use musica_util, only: error_t, error_t_c, dk => musica_dk ! Arguments - class(radiator_t), intent(inout) :: this - real(dk), target, dimension(:,:,:), intent(in) :: asymmetry_factors - type(error_t), intent(inout) :: error + class(radiator_t), intent(in) :: this + real(dk), target, contiguous, intent(out) :: asymmetry_factors(:,:,:) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c diff --git a/fortran/tuvx/tuvx.F90 b/fortran/tuvx/tuvx.F90 index d5c2ca2f..d17a5126 100644 --- a/fortran/tuvx/tuvx.F90 +++ b/fortran/tuvx/tuvx.F90 @@ -285,12 +285,12 @@ subroutine run(this, solar_zenith_angle, earth_sun_distance, & 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 + 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, contiguous, intent(inout) :: photolysis_rate_constants(:,:) ! s-1 (layer, reaction) + real(kind=dk), target, contiguous, intent(inout) :: heating_rates(:,:) ! K s-1 (layer, reaction) + type(error_t), intent(inout) :: error ! Local variables type(error_t_c) :: error_c diff --git a/fortran/util.F90 b/fortran/util.F90 index 14f97a5b..fda31594 100644 --- a/fortran/util.F90 +++ b/fortran/util.F90 @@ -747,9 +747,9 @@ subroutine copy_data( this, source, target ) use iso_c_binding, only: c_loc - class(index_mappings_t), intent(inout) :: this - real(kind=musica_dk), target, intent(in) :: source(:) - real(kind=musica_dk), target, intent(in) :: target(:) + class(index_mappings_t), intent(inout) :: this + real(kind=musica_dk), target, contiguous, intent(in) :: source(:) + real(kind=musica_dk), target, contiguous, intent(in) :: target(:) call copy_data_c( this%mappings_c_, c_loc( source ), c_loc( target ) ) From 326b5119768d5be9654baf96ae3bd6a1b757fdc8 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Wed, 6 Nov 2024 14:01:26 -0800 Subject: [PATCH 12/13] fit copy_data argument intent, add tests, update TUV-x commit (#243) --- cmake/dependencies.cmake | 2 +- fortran/test/unit/util.F90 | 16 ++++++++-------- fortran/util.F90 | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index 6a3e3b65..a9f4d2ba 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -86,7 +86,7 @@ if (MUSICA_ENABLE_TUVX AND MUSICA_BUILD_C_CXX_INTERFACE) set(TUVX_INSTALL_INCLUDE_DIR ${MUSICA_INSTALL_INCLUDE_DIR} CACHE STRING "" FORCE) set_git_default(TUVX_GIT_REPOSITORY https://github.com/NCAR/tuv-x.git) - set_git_default(TUVX_GIT_TAG v0.10.0) + set_git_default(TUVX_GIT_TAG fbe0f8aa73f6630d230c6463b603d6ba64c65dcf) FetchContent_Declare(tuvx GIT_REPOSITORY ${TUVX_GIT_REPOSITORY} diff --git a/fortran/test/unit/util.F90 b/fortran/test/unit/util.F90 index 1efafd28..05a527b6 100644 --- a/fortran/test/unit/util.F90 +++ b/fortran/test/unit/util.F90 @@ -280,7 +280,7 @@ subroutine build_and_check_index_mapping_t(config) type(mappings_t), pointer :: source_map, target_map type(index_mappings_t), pointer :: index_mappings type(error_t) :: error - real(dk), allocatable :: source_data(:), target_data(:) + real(dk) :: source_data(2,5), target_data(3,4) allocate( f_map( 2 ) ) map => mapping_t( "Test", 2 ) @@ -302,14 +302,14 @@ subroutine build_and_check_index_mapping_t(config) index_mappings => index_mappings_t( config, source_map, target_map, error ) ASSERT( error%is_success() ) - source_data = (/ 1.0_dk, 2.0_dk, 3.0_dk, 4.0_dk, 5.0_dk /) - target_data = (/ 10.0_dk, 20.0_dk, 30.0_dk, 40.0_dk /) + source_data(2,:) = (/ 1.0_dk, 2.0_dk, 3.0_dk, 4.0_dk, 5.0_dk /) + target_data(2,:) = (/ 10.0_dk, 20.0_dk, 30.0_dk, 40.0_dk /) - call index_mappings%copy_data( source_data, target_data ) - ASSERT_EQ( target_data( 1 ), 5.0_dk * 0.82_dk ) - ASSERT_EQ( target_data( 2 ), 20.0_dk ) - ASSERT_EQ( target_data( 3 ), 2.0_dk ) - ASSERT_EQ( target_data( 4 ), 40.0_dk ) + call index_mappings%copy_data( source_data(2,:), target_data(2,:) ) + ASSERT_EQ( target_data( 2, 1 ), 5.0_dk * 0.82_dk ) + ASSERT_EQ( target_data( 2, 2 ), 20.0_dk ) + ASSERT_EQ( target_data( 2, 3 ), 2.0_dk ) + ASSERT_EQ( target_data( 2, 4 ), 40.0_dk ) deallocate( index_mappings ) deallocate( source_map ) deallocate( target_map ) diff --git a/fortran/util.F90 b/fortran/util.F90 index fda31594..76e4f278 100644 --- a/fortran/util.F90 +++ b/fortran/util.F90 @@ -749,7 +749,7 @@ subroutine copy_data( this, source, target ) class(index_mappings_t), intent(inout) :: this real(kind=musica_dk), target, contiguous, intent(in) :: source(:) - real(kind=musica_dk), target, contiguous, intent(in) :: target(:) + real(kind=musica_dk), target, contiguous, intent(inout) :: target(:) call copy_data_c( this%mappings_c_, c_loc( source ), c_loc( target ) ) From 1debffba28d0601cbd646e475cd6751d06e5a8b3 Mon Sep 17 00:00:00 2001 From: Jiwon Gim <55209567+boulderdaze@users.noreply.github.com> Date: Thu, 19 Dec 2024 15:08:57 -0700 Subject: [PATCH 13/13] Update version (#246) * update version * fix tuvx tag * correctly creating backward euler * removing 3.7 --------- Co-authored-by: Kyle Shores --- .github/workflows/pip.yml | 2 +- CITATION.cff | 2 +- CMakeLists.txt | 2 +- cmake/dependencies.cmake | 4 ++-- src/micm/micm.cpp | 8 ++------ 5 files changed, 7 insertions(+), 11 deletions(-) diff --git a/.github/workflows/pip.yml b/.github/workflows/pip.yml index 724bbc0a..782e4e80 100644 --- a/.github/workflows/pip.yml +++ b/.github/workflows/pip.yml @@ -13,7 +13,7 @@ jobs: fail-fast: false matrix: platform: [windows-latest, macos-13, ubuntu-latest] - python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] # python versions: https://devguide.python.org/versions/ runs-on: ${{ matrix.platform }} diff --git a/CITATION.cff b/CITATION.cff index 65f9cbd8..0c2c1a53 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -64,4 +64,4 @@ number: 10 page: "E1743 - E1760" doi: "10.1175/BAMS-D-19-0331.1" url: "https://journals.ametsoc.org/view/journals/bams/101/10/bamsD190331.xml" -version: 0.8.0 +version: 0.9.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index a0fb629a..bb898bcf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.21) # must be on the same line so that pyproject.toml can correctly identify the version -project(musica-distribution VERSION 0.8.1) +project(musica-distribution VERSION 0.9.0) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}/cmake) set(CMAKE_USER_MAKE_RULES_OVERRIDE ${CMAKE_MODULE_PATH}/SetDefaults.cmake) diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake index a9f4d2ba..b6de32d6 100644 --- a/cmake/dependencies.cmake +++ b/cmake/dependencies.cmake @@ -62,7 +62,7 @@ endif() if (MUSICA_ENABLE_MICM AND MUSICA_BUILD_C_CXX_INTERFACE) set_git_default(MICM_GIT_REPOSITORY https://github.com/NCAR/micm.git) - set_git_default(MICM_GIT_TAG b3c462a) + set_git_default(MICM_GIT_TAG v.3.7.0) FetchContent_Declare(micm GIT_REPOSITORY ${MICM_GIT_REPOSITORY} @@ -86,7 +86,7 @@ if (MUSICA_ENABLE_TUVX AND MUSICA_BUILD_C_CXX_INTERFACE) set(TUVX_INSTALL_INCLUDE_DIR ${MUSICA_INSTALL_INCLUDE_DIR} CACHE STRING "" FORCE) set_git_default(TUVX_GIT_REPOSITORY https://github.com/NCAR/tuv-x.git) - set_git_default(TUVX_GIT_TAG fbe0f8aa73f6630d230c6463b603d6ba64c65dcf) + set_git_default(TUVX_GIT_TAG v0.10.1) FetchContent_Declare(tuvx GIT_REPOSITORY ${TUVX_GIT_REPOSITORY} diff --git a/src/micm/micm.cpp b/src/micm/micm.cpp index d10656b6..d4bde009 100644 --- a/src/micm/micm.cpp +++ b/src/micm/micm.cpp @@ -359,15 +359,11 @@ namespace musica solver_parameters_ = std::make_unique(solver_config.GetSolverParams()); auto solver = std::make_unique( - micm::SolverBuilder< + micm::CpuSolverBuilder< micm::BackwardEulerSolverParameters, micm::VectorMatrix, micm::SparseMatrix>, - micm::ProcessSet, - micm::LinearSolver< - micm::SparseMatrix>, - micm::LuDecomposition>, - VectorState>(micm::BackwardEulerSolverParameters()) + micm::LuDecompositionDoolittle>(micm::BackwardEulerSolverParameters()) .SetSystem(solver_parameters_->system_) .SetReactions(solver_parameters_->processes_) .SetNumberOfGridCells(num_grid_cells_)