279 make solve function for ode solvers thread safe

.github/workflows/ubuntu.yml

@@ -48,7 +48,7 @@ jobs:
         run: cmake -S . -B build -D CMAKE_BUILD_TYPE=${{ matrix.build_type }}
 
       - name: Build
-        run: cmake --build build --parallel 10
+        run: cmake --build build --verbose
 
       - name: Run tests
         run: |

CMakeLists.txt

@@ -51,6 +51,12 @@ endif()
 
 add_compile_definitions(DEFAULT_VECTOR_SIZE=${DEFAULT_VECTOR_MATRIX_SIZE})
 
+# on ubuntu with clang, an incorrect version of the c++ standard library was being linked
+if (${CMAKE_HOST_SYSTEM_NAME} MATCHES "Linux" AND "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+    # If the compiler is Clang, use libc++
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++")
+endif()
+
 ################################################################################
 # Dependencies
 

cmake/test_util.cmake

@@ -25,14 +25,6 @@ function(create_standard_test)
     target_link_libraries(test_${TEST_NAME} PUBLIC ${library})
   endforeach()
 
-  if(ENABLE_JSON)
-    target_link_libraries(test_${TEST_NAME} PRIVATE nlohmann_json::nlohmann_json)
-  endif()
-
-  if(ENABLE_LLVM)
-    target_link_libraries(test_${TEST_NAME} PRIVATE ${llvm_libs})
-  endif()
-
   if(NOT DEFINED TEST_WORKING_DIRECTORY)
     set(TEST_WORKING_DIRECTORY "${CMAKE_BINARY_DIR}")
   endif()

docker/Dockerfile

@@ -20,7 +20,7 @@ RUN mkdir /build \
         -D CMAKE_BUILD_TYPE=debug \
         -D ENABLE_CLANG_TIDY:BOOL=FALSE \
         ../micm \
-      && make install -j 8
+      && make install -j
 
 # now test if we can use the installed files
 RUN cd /micm/test/integration/cmake/find_package \

docs/source/conf.py

@@ -66,6 +66,8 @@
         "json_url": "https://ncar.github.io/micm/_static/switcher.json",
         "version_match": release,
     },
+   "pygment_light_style": "tango",
+   "pygment_dark_style": "monokai"
 }
 
 html_css_files = [

docs/source/user_guide/user_defined_rate_constant_tutorial.rst

@@ -194,7 +194,7 @@ Finally, set and upate the rate constants as needed:
              // so we need to track how much time the solver was able to integrate for and continue
              // solving until we finish
              double elapsed_solve_time = 0;
-          +  state.SetCustomRateParameter("PHOTO.my photolysis rate", photo_rate);
+             +state.SetCustomRateParameter("PHOTO.my photolysis rate", photo_rate);
 
              while (elapsed_solve_time < time_step)
              {
@@ -205,7 +205,7 @@ Finally, set and upate the rate constants as needed:
              }
 
              print_state(time_step * (i + 1), state);
-          +   photo_rate *= 1.5;
+             +photo_rate *= 1.5;
            }
 
 And this is final output. Notice that the concentration of G ends up much higher than in 

include/micm/process/arrhenius_rate_constant.hpp

@@ -47,6 +47,11 @@ namespace micm
     /// @return A rate constant based off of the conditions in the system
     double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters) const override;
 
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions) const override;
+
     double calculate(const double& temperature, const double& pressure) const;
   };
 
@@ -72,6 +77,11 @@ namespace micm
     return calculate(conditions.temperature_, conditions.pressure_);
   }
 
+  inline double ArrheniusRateConstant::calculate(const Conditions& conditions) const
+  {
+    return calculate(conditions.temperature_, conditions.pressure_);
+  }
+
   inline double ArrheniusRateConstant::calculate(const double& temperature, const double& pressure) const
   {
     return parameters_.A_ * std::exp(parameters_.C_ / temperature) * std::pow(temperature / parameters_.D_, parameters_.B_) *

include/micm/process/branched_rate_constant.hpp

@@ -54,6 +54,11 @@ namespace micm
     /// @return A rate constant based off of the conditions in the system
     double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters) const override;
 
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions) const override;
+
     /// @brief Calculate the rate constant
     /// @param temperature Temperature in [K]
     /// @param air_number_density Number density in [mol m-3]
@@ -92,6 +97,11 @@ namespace micm
     return calculate(conditions.temperature_, conditions.air_density_);
   }
 
+  inline double BranchedRateConstant::calculate(const Conditions& conditions) const
+  {
+    return calculate(conditions.temperature_, conditions.air_density_);
+  }
+
   inline double BranchedRateConstant::calculate(const double& temperature, const double& air_number_density) const
   {
     double pre = parameters_.X_ * std::exp(-parameters_.Y_ / temperature);

include/micm/process/cuda_process_set.hpp

@@ -5,12 +5,8 @@
 
 #include <micm/process/process_set.hpp>
 #include <micm/util/cuda_param.hpp>
+#include <micm/process/cuda_process_set.cuh>
 
-#ifdef USE_CUDA
-#  include <micm/process/cuda_process_set.cuh>
-#endif
-
-#ifdef USE_CUDA
 namespace micm
 {
   /// @brief A GPU-based implementation of ProcessSet
@@ -20,8 +16,7 @@ namespace micm
     /// @brief Create a process set calculator for a given set of processes
     /// @param processes Processes to create calculator for
     /// @param state Solver state
-    template<template<class> class MatrixPolicy>
-    CudaProcessSet(const std::vector<Process>& processes, const State<MatrixPolicy>& state);
+    CudaProcessSet(const std::vector<Process>& processes, const std::map<std::string, std::size_t>& variable_map);
 
     template<template<class> typename MatrixPolicy>
     requires VectorizableDense<MatrixPolicy<double>> std::chrono::nanoseconds AddForcingTerms(
@@ -39,9 +34,8 @@ namespace micm
     const;
   };
 
-  template<template<class> class MatrixPolicy>
-  inline CudaProcessSet::CudaProcessSet(const std::vector<Process>& processes, const State<MatrixPolicy>& state)
-      : ProcessSet(processes, state)
+  inline CudaProcessSet::CudaProcessSet(const std::vector<Process>& processes, const std::map<std::string, std::size_t>& variable_map)
+      : ProcessSet(processes, variable_map)
   {
   }
 
@@ -73,6 +67,7 @@ namespace micm
     std::chrono::nanoseconds kernel_duration = micm::cuda::AddForcingTermsKernelDriver(matrix, processSet);
     return kernel_duration;  // time performance of kernel function
   }
+
   template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
   requires VectorizableDense<MatrixPolicy<double>> && VectorizableSparse<SparseMatrixPolicy<double>>
   inline std::chrono::nanoseconds CudaProcessSet::AddJacobianTerms(
@@ -104,5 +99,4 @@ namespace micm
     std::chrono::nanoseconds kernel_duration = micm::cuda::AddJacobianTermsKernelDriver(matrix, sparseMatrix, processSet);
     return kernel_duration;  // time performance of kernel function
   }
-}  // namespace micm
-#endif
+}  // namespace micm

include/micm/process/jit_process_set.hpp

@@ -33,11 +33,10 @@ namespace micm
     /// @param compiler JIT compiler
     /// @param processes Processes to create calculator for
     /// @param state Solver state
-    template<template<class> class MatrixPolicy>
     JitProcessSet(
         std::shared_ptr<JitCompiler> compiler,
         const std::vector<Process> &processes,
-        const State<MatrixPolicy> &state);
+        const std::map<std::string, std::size_t>& variable_map);
 
     ~JitProcessSet();
 
@@ -103,20 +102,15 @@ namespace micm
   }
 
   template<std::size_t L>
-  template<template<class> class MatrixPolicy>
   inline JitProcessSet<L>::JitProcessSet(
       std::shared_ptr<JitCompiler> compiler,
       const std::vector<Process> &processes,
-      const State<MatrixPolicy> &state)
-      : ProcessSet(processes, state),
+      const std::map<std::string, std::size_t>& variable_map)
+      : ProcessSet(processes, variable_map),
         compiler_(compiler)
   {
     forcing_function_ = NULL;
     jacobian_function_ = NULL;
-    if (state.variables_.size() != L || state.variables_.GroupVectorSize() != L)
-    {
-      throw std::runtime_error("Invalid state for JitProcessSet. Check the the VectorMatrix template parameters.");
-    }
     this->GenerateForcingFunction();
   }
 

include/micm/process/process.hpp

@@ -42,14 +42,14 @@ namespace micm
     /// @brief Update the solver state rate constants
     /// @param processes The set of processes being solved
     /// @param state The solver state to update
-    template<template<class> class MatrixPolicy>
+    template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
     requires(!VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
         const std::vector<Process>& processes,
-        State<MatrixPolicy>& state);
-    template<template<class> class MatrixPolicy>
+        State<MatrixPolicy, SparseMatrixPolicy>& state);
+    template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
     requires(VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
         const std::vector<Process>& processes,
-        State<MatrixPolicy>& state);
+        State<MatrixPolicy, SparseMatrixPolicy>& state);
 
     friend class ProcessBuilder;
     static ProcessBuilder create();
@@ -105,10 +105,10 @@ namespace micm
     ProcessBuilder& phase(const Phase& phase);
   };
 
-  template<template<class> class MatrixPolicy>
+  template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
   requires(!VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
       const std::vector<Process>& processes,
-      State<MatrixPolicy>& state)
+      State<MatrixPolicy, SparseMatrixPolicy>& state)
   {
     for (std::size_t i{}; i < state.custom_rate_parameters_.size(); ++i)
     {
@@ -128,10 +128,10 @@ namespace micm
     }
   }
 
-  template<template<class> class MatrixPolicy>
+  template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
   requires(VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
       const std::vector<Process>& processes,
-      State<MatrixPolicy>& state)
+      State<MatrixPolicy, SparseMatrixPolicy>& state)
   {
     const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector();
     auto& v_rate_constants = state.rate_constants_.AsVector();

include/micm/process/process_set.hpp

@@ -29,9 +29,8 @@ namespace micm
 
     /// @brief Create a process set calculator for a given set of processes
     /// @param processes Processes to create calculator for
-    /// @param state Solver state
-    template<template<class> class MatrixPolicy>
-    ProcessSet(const std::vector<Process>& processes, const State<MatrixPolicy>& state);
+    /// @param StateParameters Solver state
+    ProcessSet(const std::vector<Process>& processes, const std::map<std::string, std::size_t>& variable_map);
 
     /// @brief Return the full set of non-zero Jacobian elements for the set of processes
     /// @return Jacobian elements as a set of index pairs
@@ -74,8 +73,7 @@ namespace micm
         SparseMatrixPolicy<double>& jacobian) const;
   };
 
-  template<template<class> class MatrixPolicy>
-  inline ProcessSet::ProcessSet(const std::vector<Process>& processes, const State<MatrixPolicy>& state)
+  inline ProcessSet::ProcessSet(const std::vector<Process>& processes, const std::map<std::string, std::size_t>& variable_map)
       : number_of_reactants_(),
         reactant_ids_(),
         number_of_products_(),
@@ -90,14 +88,14 @@ namespace micm
       {
         if (reactant.IsParameterized())
           continue;  // Skip reactants that are parameterizations
-        reactant_ids_.push_back(state.variable_map_.at(reactant.name_));
+        reactant_ids_.push_back(variable_map.at(reactant.name_));
         ++number_of_reactants;
       }
       for (auto& product : process.products_)
       {
         if (product.first.IsParameterized())
           continue;  // Skip products that are parameterizations
-        product_ids_.push_back(state.variable_map_.at(product.first.name_));
+        product_ids_.push_back(variable_map.at(product.first.name_));
         yields_.push_back(product.second);
         ++number_of_products;
       }

include/micm/process/rate_constant.hpp

@@ -36,6 +36,15 @@ namespace micm
       return 0;
     }
 
+    /// @brief Calculate the rate constant for a set of conditions
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @param custom_parameters User defined rate constant parameters
+    /// @return The reaction rate constant
+    virtual double calculate(const Conditions& conditions) const
+    {
+      return 0;
+    }
+
     /// @brief Calculate the rate constant for a set of conditions
     /// @param conditions The current environmental conditions of the chemical system
     /// @param custom_parameters User defined rate constant parameters

include/micm/process/surface_rate_constant.hpp

@@ -56,6 +56,11 @@ namespace micm
     /// @param custom_parameters User-defined rate constant parameters
     /// @return A rate constant based off of the conditions in the system
     double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters) const override;
+
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions) const override;
   };
 
   inline SurfaceRateConstant::SurfaceRateConstant(const SurfaceRateConstantParameters& parameters)
@@ -70,6 +75,12 @@ namespace micm
     return std::unique_ptr<RateConstant>{ new SurfaceRateConstant{ *this } };
   }
 
+  inline double SurfaceRateConstant::calculate(
+      const Conditions& conditions) const
+  {
+    throw std::runtime_error("Surface rate constants must be supplied with a radius and number density using the alternative calculate function");
+  }
+
   inline double SurfaceRateConstant::calculate(
       const Conditions& conditions,
       std::vector<double>::const_iterator custom_parameters) const

include/micm/process/ternary_chemical_activation_rate_constant.hpp

@@ -52,6 +52,11 @@ namespace micm
     /// @return A rate constant based off of the conditions in the system
     double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters) const override;
 
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions) const override;
+
     /// @brief Calculate the rate constant
     /// @param temperature Temperature in [K]
     /// @param air_number_density Number density in [mol m-3]
@@ -75,6 +80,12 @@ namespace micm
     return std::unique_ptr<RateConstant>{ new TernaryChemicalActivationRateConstant{ *this } };
   }
 
+  inline double TernaryChemicalActivationRateConstant::calculate(const Conditions& conditions) const
+  {
+    double val = calculate(conditions.temperature_, conditions.air_density_);
+    return val;
+  }
+
   inline double TernaryChemicalActivationRateConstant::calculate(
       const Conditions& conditions,
       std::vector<double>::const_iterator custom_parameters) const

include/micm/process/troe_rate_constant.hpp

@@ -52,6 +52,11 @@ namespace micm
     /// @return A rate constant based off of the conditions in the system
     double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters) const override;
 
+    /// @brief Calculate the rate constant
+    /// @param conditions The current environmental conditions of the chemical system
+    /// @return A rate constant based off of the conditions in the system
+    double calculate(const Conditions& conditions) const override;
+
     /// @brief Calculate the rate constant
     /// @param temperature Temperature in [K]
     /// @param air_number_density Number density in [mol m-3]
@@ -74,6 +79,11 @@ namespace micm
     return std::unique_ptr<RateConstant>{ new TroeRateConstant{ *this } };
   }
 
+  inline double TroeRateConstant::calculate(const Conditions& conditions) const
+  {
+    return calculate(conditions.temperature_, conditions.air_density_);
+  }
+
   inline double TroeRateConstant::calculate(
       const Conditions& conditions,
       std::vector<double>::const_iterator custom_parameters) const