Python error checking (#189)

* adding error checking, some debug statements for the python wrapper

* adding error checking to python, updating micm tag

* getting tests to pass

* correcting a fortran test

* newest tuvx commit

* I don't know what's wrong with nvidia

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.7.0)
+project(musica-distribution VERSION 0.7.1)
 
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH};${PROJECT_SOURCE_DIR}/cmake)
 set(CMAKE_USER_MAKE_RULES_OVERRIDE ${CMAKE_MODULE_PATH}/SetDefaults.cmake)

cmake/dependencies.cmake

@@ -51,7 +51,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 v3.5.0)
+  set_git_default(MICM_GIT_TAG 6b1c58a9be14095e3f3c6df403c91c4e800f23de)
 
   FetchContent_Declare(micm
       GIT_REPOSITORY ${MICM_GIT_REPOSITORY}

configs/chapman/species.json

@@ -3,7 +3,7 @@
         {
             "name": "M",
             "type": "CHEM_SPEC",
-            "tracer type": "CONSTANT"
+            "tracer type": "THIRD_BODY"
         },
         {
             "name": "O2",

fortran/test/fetch_content_integration/test_micm_api.F90

@@ -30,7 +30,7 @@ subroutine test_api()
     real(c_double)                :: pressure
     real(c_double)                :: air_density
     integer(c_int)                :: num_concentrations, num_user_defined_reaction_rates
-    real(c_double), dimension(5)  :: concentrations 
+    real(c_double), dimension(4)  :: concentrations 
     real(c_double), dimension(3)  :: user_defined_reaction_rates 
     character(len=256)            :: config_path
     integer(c_int)                :: solver_type
@@ -52,8 +52,8 @@ subroutine test_api()
     temperature = 272.5
     pressure = 101253.4
     air_density = pressure / ( GAS_CONSTANT * temperature )
-    num_concentrations = 5
-    concentrations = (/ 0.75, 0.4, 0.8, 0.01, 0.02 /)
+    num_concentrations = 4
+    concentrations = (/ 0.4, 0.8, 0.01, 0.02 /)
     num_user_defined_reaction_rates = 3
     user_defined_reaction_rates = (/ 0.1, 0.2, 0.3 /)
 
@@ -109,21 +109,16 @@ subroutine test_api()
     ASSERT( logical( bool_value ) )
 
     string_value = micm%get_species_property_string( "O3", "missing property", error )
-    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, \
-                      MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
+    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
     double_value = micm%get_species_property_double( "O3", "missing property", error )
-    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, \
-                      MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
+    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
     int_value = micm%get_species_property_int( "O3", "missing property", error )
-    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, \
-                      MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
+    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
     bool_value = micm%get_species_property_bool( "O3", "missing property", error )
-    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, \
-                      MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
+    ASSERT( error%is_error( MICM_ERROR_CATEGORY_SPECIES, MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND ) )
     deallocate( micm )
     micm => micm_t( "configs/invalid", solver_type, num_grid_cells, error )
-    ASSERT( error%is_error( MICM_ERROR_CATEGORY_CONFIGURATION, \
-                      MICM_CONFIGURATION_ERROR_CODE_INVALID_FILE_PATH ) )
+    ASSERT( error%is_error( MICM_ERROR_CATEGORY_CONFIGURATION, MICM_CONFIGURATION_ERROR_CODE_INVALID_FILE_PATH ) )
     ASSERT( .not. associated( micm ) )
 
     write(*,*) "[test micm fort api] Finished."

include/musica/micm.hpp

@@ -216,6 +216,7 @@ namespace musica
     using RosenbrockVectorType = typename micm::RosenbrockSolverParameters::
         template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixVector, micm::LuDecomposition>>;
     using Rosenbrock = micm::Solver<RosenbrockVectorType, micm::State<DenseMatrixVector, SparseMatrixVector>>;
+    using VectorState = micm::State<DenseMatrixVector, SparseMatrixVector>;
     std::unique_ptr<Rosenbrock> rosenbrock_;
 
     /// @brief Standard-ordered Rosenbrock solver type
@@ -224,6 +225,7 @@ namespace musica
     using RosenbrockStandardType = typename micm::RosenbrockSolverParameters::
         template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixStandard, micm::LuDecomposition>>;
     using RosenbrockStandard = micm::Solver<RosenbrockStandardType, micm::State<DenseMatrixStandard, SparseMatrixStandard>>;
+    using StandardState = micm::State<DenseMatrixStandard, SparseMatrixStandard>;
     std::unique_ptr<RosenbrockStandard> rosenbrock_standard_;
 
    private:

python/test/test_chapman.py

@@ -10,7 +10,7 @@ def test_micm_solve(self):
         pressure = 101253.3
         GAS_CONSTANT = 8.31446261815324
         air_density = pressure / (GAS_CONSTANT * temperature)
-        concentrations = [0.75, 0.4, 0.8, 0.01, 0.02]
+        concentrations = [0.4, 0.8, 0.01, 0.02]
 
         solver = musica.create_solver(
             "configs/chapman",
@@ -41,16 +41,15 @@ def test_micm_solve(self):
 
         self.assertEqual(
             ordering, {
-                'M': 0, 'O': 2, 'O1D': 1, 'O2': 3, 'O3': 4})
+                'O': 1, 'O1D': 0, 'O2': 2, 'O3': 3})
         self.assertEqual(
             rate_constant_ordering, {
                 'PHOTO.R1': 0, 'PHOTO.R3': 1, 'PHOTO.R5': 2})
 
-        self.assertEqual(concentrations[0], 0.75)
-        self.assertNotEqual(concentrations[1], 0.4)
-        self.assertNotEqual(concentrations[2], 0.8)
-        self.assertNotEqual(concentrations[3], 0.01)
-        self.assertNotEqual(concentrations[4], 0.02)
+        self.assertNotEqual(concentrations[0], 0.4)
+        self.assertNotEqual(concentrations[1], 0.8)
+        self.assertNotEqual(concentrations[2], 0.01)
+        self.assertNotEqual(concentrations[3], 0.02)
 
 
 if __name__ == '__main__':

python/wrapper.cpp

@@ -78,6 +78,12 @@ PYBIND11_MODULE(musica, m)
             &solver_state,
             &solver_stats,
             &error);
+        if (!musica::IsSuccess(error))
+        {
+          std::string message = "Error solving system: " + std::string(error.message_.value_);
+          DeleteError(&error);
+          throw std::runtime_error(message);
+        }
 
         // Update the concentrations list after solving
         for (std::size_t i = 0; i < concentrations_cpp.size(); ++i)

src/micm/micm.cpp

@@ -248,7 +248,8 @@ namespace musica
               micm::ProcessSet,
               micm::LinearSolver<
                   micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>,
-                  micm::LuDecomposition>>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
+                  micm::LuDecomposition>,
+              VectorState>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters())
               .SetSystem(solver_parameters_->system_)
               .SetReactions(solver_parameters_->processes_)
               .SetNumberOfGridCells(num_grid_cells_)

src/test/unit/micm/micm_c_api.cpp

@@ -97,7 +97,7 @@ TEST_F(MicmCApiTest, GetSpeciesOrdering)
   Mapping* species_ordering = GetSpeciesOrdering(micm, &array_size, &error);
   ASSERT_TRUE(IsSuccess(error));
   DeleteError(&error);
-  ASSERT_EQ(array_size, 5);
+  ASSERT_EQ(array_size, 4);
   bool found = false;
   for (std::size_t i = 0; i < array_size; i++)
   {
@@ -130,16 +130,6 @@ TEST_F(MicmCApiTest, GetSpeciesOrdering)
   ASSERT_TRUE(found);
   found = false;
   for (std::size_t i = 0; i < array_size; i++)
-  {
-    if (strcmp(species_ordering[i].name_.value_, "M") == 0)
-    {
-      found = true;
-      break;
-    }
-  }
-  ASSERT_TRUE(found);
-  found = false;
-  for (std::size_t i = 0; i < array_size; i++)
   {
     if (strcmp(species_ordering[i].name_.value_, "O1D") == 0)
     {
@@ -201,8 +191,8 @@ TEST_F(MicmCApiTest, SolveUsingVectorOrderedRosenbrock)
   double pressure = 101253.3;
   constexpr double GAS_CONSTANT = 8.31446261815324;  // J mol-1 K-1
   double air_density = pressure / (GAS_CONSTANT * temperature);
-  int num_concentrations = 5;
-  double concentrations[] = { 0.75, 0.4, 0.8, 0.01, 0.02 };
+  int num_concentrations = 4;
+  double concentrations[] = { 0.4, 0.8, 0.01, 0.02 };
   std::size_t num_user_defined_reaction_rates = 3;
   double user_defined_reaction_rates[] = { 0.1, 0.2, 0.3 };
   String solver_state;
@@ -234,11 +224,10 @@ TEST_F(MicmCApiTest, SolveUsingVectorOrderedRosenbrock)
   ASSERT_TRUE(IsSuccess(error));
 
   // Add assertions to check the solved concentrations
-  ASSERT_EQ(concentrations[0], 0.75);
-  ASSERT_NE(concentrations[1], 0.4);
-  ASSERT_NE(concentrations[2], 0.8);
-  ASSERT_NE(concentrations[3], 0.01);
-  ASSERT_NE(concentrations[4], 0.02);
+  ASSERT_NE(concentrations[0], 0.4);
+  ASSERT_NE(concentrations[1], 0.8);
+  ASSERT_NE(concentrations[2], 0.01);
+  ASSERT_NE(concentrations[3], 0.02);
 
   std::cout << "Solver state: " << solver_state.value_ << std::endl;
   std::cout << "Function Calls: " << solver_stats.function_calls_ << std::endl;
@@ -269,8 +258,8 @@ TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock)
   double pressure = 101253.3;
   constexpr double GAS_CONSTANT = 8.31446261815324;  // J mol-1 K-1
   double air_density = pressure / (GAS_CONSTANT * temperature);
-  int num_concentrations = 5;
-  double concentrations[] = { 0.75, 0.4, 0.8, 0.01, 0.02 };
+  int num_concentrations = 4;
+  double concentrations[] = { 0.4, 0.8, 0.01, 0.02 };
   std::size_t num_user_defined_reaction_rates = 3;
   double user_defined_reaction_rates[] = { 0.1, 0.2, 0.3 };
   String solver_state;
@@ -300,12 +289,10 @@ TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock)
       &error);
   ASSERT_TRUE(IsSuccess(error));
 
-  // Add assertions to check the solved concentrations
-  ASSERT_EQ(concentrations[0], 0.75);
-  ASSERT_NE(concentrations[1], 0.4);
-  ASSERT_NE(concentrations[2], 0.8);
-  ASSERT_NE(concentrations[3], 0.01);
-  ASSERT_NE(concentrations[4], 0.02);
+  ASSERT_NE(concentrations[0], 0.4);
+  ASSERT_NE(concentrations[1], 0.8);
+  ASSERT_NE(concentrations[2], 0.01);
+  ASSERT_NE(concentrations[3], 0.02);
 
   std::cout << "Solver state: " << solver_state.value_ << std::endl;
   std::cout << "Function Calls: " << solver_stats.function_calls_ << std::endl;