From f7a69ef440aabaa74eb1786830250d50f1782b6e Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 5 Nov 2024 10:31:14 -0800 Subject: [PATCH] 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); }