Merge pull request #314 from NCAR/develop-245-terminator

Develop 245 terminator

include/micm/process/process.hpp

@@ -36,13 +36,11 @@ namespace micm
     /// @param processes The set of processes being solved
     /// @param state The solver state to update
     template<template<class> class MatrixPolicy>
-    requires(!VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
-        const std::vector<Process>& processes,
-        State<MatrixPolicy>& state);
+      requires(!VectorizableDense<MatrixPolicy<double>>)
+    static void UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state);
     template<template<class> class MatrixPolicy>
-    requires(VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
-        const std::vector<Process>& processes,
-        State<MatrixPolicy>& state);
+      requires(VectorizableDense<MatrixPolicy<double>>)
+    static void UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state);
 
     friend class ProcessBuilder;
     static ProcessBuilder create();
@@ -92,9 +90,8 @@ namespace micm
   };
 
   template<template<class> class MatrixPolicy>
-  requires(!VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
-      const std::vector<Process>& processes,
-      State<MatrixPolicy>& state)
+    requires(!VectorizableDense<MatrixPolicy<double>>)
+  void Process::UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state)
   {
     for (std::size_t i{}; i < state.custom_rate_parameters_.size(); ++i)
     {
@@ -103,17 +100,20 @@ namespace micm
       std::size_t i_rate_constant = 0;
       for (auto& process : processes)
       {
+        double fixed_reactants = 1.0;
+        for (auto& reactant : process.reactants_)
+          if (reactant.IsParameterized())
+            fixed_reactants *= reactant.parameterize_(state.conditions_[i]);
         state.rate_constants_[i][(i_rate_constant++)] =
-            process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter);
+            process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter) * fixed_reactants;
         custom_parameters_iter += process.rate_constant_->SizeCustomParameters();
       }
     }
   }
 
   template<template<class> class MatrixPolicy>
-  requires(VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
-      const std::vector<Process>& processes,
-      State<MatrixPolicy>& state)
+    requires(VectorizableDense<MatrixPolicy<double>>)
+  void Process::UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state)
   {
     const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector();
     auto& v_rate_constants = state.rate_constants_.AsVector();
@@ -133,8 +133,13 @@ namespace micm
             params[i_param] = v_custom_parameters[offset_params + i_param * L + i_cell];
           }
           std::vector<double>::const_iterator custom_parameters_iter = params.begin();
+          double fixed_reactants = 1.0;
+          for (auto& reactant : process.reactants_)
+            if (reactant.IsParameterized())
+              fixed_reactants *= reactant.parameterize_(state.conditions_[i_group * L + i_cell]);
           v_rate_constants[offset_rc + i_cell] =
-              process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter);
+              process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter) *
+              fixed_reactants;
         }
         offset_params += params.size() * L;
         offset_rc += L;

include/micm/process/process_set.hpp

@@ -47,12 +47,13 @@ namespace micm
     /// @param state_variables Current state variable values (grid cell, state variable)
     /// @param forcing Forcing terms for each state variable (grid cell, state variable)
     template<template<class> typename MatrixPolicy>
-    requires(!VectorizableDense<MatrixPolicy<double>>) void AddForcingTerms(
+      requires(!VectorizableDense<MatrixPolicy<double>>)
+    void AddForcingTerms(
         const MatrixPolicy<double>& rate_constants,
         const MatrixPolicy<double>& state_variables,
         MatrixPolicy<double>& forcing) const;
     template<template<class> typename MatrixPolicy>
-    requires VectorizableDense<MatrixPolicy<double>>
+      requires VectorizableDense<MatrixPolicy<double>>
     void AddForcingTerms(
         const MatrixPolicy<double>& rate_constants,
         const MatrixPolicy<double>& state_variables,
@@ -63,12 +64,14 @@ namespace micm
     /// @param state_variables Current state variable values (grid cell, state variable)
     /// @param jacobian Jacobian matrix for the system (grid cell, dependent variable, independent variable)
     template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
-    requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>) void AddJacobianTerms(
+      requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>)
+    void AddJacobianTerms(
         const MatrixPolicy<double>& rate_constants,
         const MatrixPolicy<double>& state_variables,
         SparseMatrixPolicy<double>& jacobian) const;
     template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
-    requires(VectorizableDense<MatrixPolicy<double>>&& VectorizableSparse<SparseMatrixPolicy<double>>) void AddJacobianTerms(
+      requires(VectorizableDense<MatrixPolicy<double>> && VectorizableSparse<SparseMatrixPolicy<double>>)
+    void AddJacobianTerms(
         const MatrixPolicy<double>& rate_constants,
         const MatrixPolicy<double>& state_variables,
         SparseMatrixPolicy<double>& jacobian) const;
@@ -84,17 +87,25 @@ namespace micm
   {
     for (auto& process : processes)
     {
-      number_of_reactants_.push_back(process.reactants_.size());
-      number_of_products_.push_back(process.products_.size());
+      std::size_t number_of_reactants = 0;
+      std::size_t number_of_products = 0;
       for (auto& reactant : process.reactants_)
       {
+        if (reactant.IsParameterized())
+          continue;  // Skip reactants that are parameterizations
         reactant_ids_.push_back(state.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_));
         yields_.push_back(product.second);
+        ++number_of_products;
       }
+      number_of_reactants_.push_back(number_of_reactants);
+      number_of_products_.push_back(number_of_products);
     }
   };
 
@@ -147,7 +158,8 @@ namespace micm
   }
 
   template<template<class> typename MatrixPolicy>
-  requires(!VectorizableDense<MatrixPolicy<double>>) inline void ProcessSet::AddForcingTerms(
+    requires(!VectorizableDense<MatrixPolicy<double>>)
+  inline void ProcessSet::AddForcingTerms(
       const MatrixPolicy<double>& rate_constants,
       const MatrixPolicy<double>& state_variables,
       MatrixPolicy<double>& forcing) const
@@ -183,7 +195,7 @@ namespace micm
   };
 
   template<template<class> typename MatrixPolicy>
-  requires VectorizableDense<MatrixPolicy<double>>
+    requires VectorizableDense<MatrixPolicy<double>>
   inline void ProcessSet::AddForcingTerms(
       const MatrixPolicy<double>& rate_constants,
       const MatrixPolicy<double>& state_variables,
@@ -224,12 +236,11 @@ namespace micm
   }
 
   template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
-  requires(
-      !VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>) inline void ProcessSet::
-      AddJacobianTerms(
-          const MatrixPolicy<double>& rate_constants,
-          const MatrixPolicy<double>& state_variables,
-          SparseMatrixPolicy<double>& jacobian) const
+    requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>)
+  inline void ProcessSet::AddJacobianTerms(
+      const MatrixPolicy<double>& rate_constants,
+      const MatrixPolicy<double>& state_variables,
+      SparseMatrixPolicy<double>& jacobian) const
   {
     auto cell_jacobian = jacobian.AsVector().begin();
 
@@ -271,11 +282,11 @@ namespace micm
   }
 
   template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
-  requires(VectorizableDense<MatrixPolicy<double>>&& VectorizableSparse<SparseMatrixPolicy<double>>) inline void ProcessSet::
-      AddJacobianTerms(
-          const MatrixPolicy<double>& rate_constants,
-          const MatrixPolicy<double>& state_variables,
-          SparseMatrixPolicy<double>& jacobian) const
+    requires(VectorizableDense<MatrixPolicy<double>> && VectorizableSparse<SparseMatrixPolicy<double>>)
+  inline void ProcessSet::AddJacobianTerms(
+      const MatrixPolicy<double>& rate_constants,
+      const MatrixPolicy<double>& state_variables,
+      SparseMatrixPolicy<double>& jacobian) const
   {
     const auto& v_rate_constants = rate_constants.AsVector();
     const auto& v_state_variables = state_variables.AsVector();

include/micm/solver/rosenbrock.hpp

@@ -50,9 +50,11 @@ namespace micm
     double rejection_factor_decrease_{ 0.1 };                     // used to decrease the step after 2 successive rejections
     double safety_factor_{ 0.9 };                                 // safety factor in new step size computation
 
-    double h_min_{ 0 };        // step size min
-    double h_max_{ 0.5 };      // step size max
-    double h_start_{ 0.005 };  // step size start
+    double h_min_{ 0.0 };  // step size min [s]
+    double h_max_{
+      0.0
+    };  // step size max [s] (if zero or greater than the solver time-step, the time-step passed to the solver will be used)
+    double h_start_{ 0.0 };  // step size start [s] (if zero, 1.0e-6 will be used, if greater than h_max, h_max will be used)
 
     // Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
     // or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
@@ -173,7 +175,8 @@ namespace micm
       uint64_t decompositions{};
       /// @brief The number of linear solvers
       uint64_t solves{};
-      /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix is never singular
+      /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix
+      /// is never singular
       uint64_t singular{};
       /// @brief The cumulative amount of time spent calculating the forcing function
       std::chrono::duration<double, std::nano> total_forcing_time{};

include/micm/solver/rosenbrock.inl

@@ -1,17 +1,20 @@
 // Copyright (C) 2023 National Center for Atmospheric Research
 // SPDX-License-Identifier: Apache-2.0
 
-#define TIMED_METHOD(assigned_increment, time_it, method, ...) \
-    { \
-      if constexpr (time_it) { \
-          auto start = std::chrono::high_resolution_clock::now(); \
-          method(__VA_ARGS__); \
-          auto end = std::chrono::high_resolution_clock::now(); \
-          assigned_increment += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); \
-      } else { \
-          method(__VA_ARGS__); \
-      } \
-    } \
+#define TIMED_METHOD(assigned_increment, time_it, method, ...)                                 \
+  {                                                                                            \
+    if constexpr (time_it)                                                                     \
+    {                                                                                          \
+      auto start = std::chrono::high_resolution_clock::now();                                  \
+      method(__VA_ARGS__);                                                                     \
+      auto end = std::chrono::high_resolution_clock::now();                                    \
+      assigned_increment += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); \
+    }                                                                                          \
+    else                                                                                       \
+    {                                                                                          \
+      method(__VA_ARGS__);                                                                     \
+    }                                                                                          \
+  }
 
 namespace micm
 {
@@ -517,9 +520,9 @@ namespace micm
     MatrixPolicy<double> forcing(Y.size(), Y[0].size(), 0.0);
     MatrixPolicy<double> temp(Y.size(), Y[0].size(), 0.0);
     std::vector<MatrixPolicy<double>> K{};
-
-    parameters_.h_max_ = time_step;
-    parameters_.h_start_ = std::max(parameters_.h_min_, delta_min_);
+    const double h_max = parameters_.h_max_ == 0.0 ? time_step : std::min(time_step, parameters_.h_max_);
+    const double h_start =
+        parameters_.h_start_ == 0.0 ? std::max(parameters_.h_min_, delta_min_) : std::min(h_max, parameters_.h_start_);
 
     stats_.Reset();
     UpdateState(state);
@@ -528,12 +531,14 @@ namespace micm
       K.push_back(MatrixPolicy<double>(Y.size(), Y[0].size(), 0.0));
 
     double present_time = 0.0;
-    double H =
-        std::min(std::max(std::abs(parameters_.h_min_), std::abs(parameters_.h_start_)), std::abs(parameters_.h_max_));
+    double H = std::min(std::max(std::abs(parameters_.h_min_), std::abs(h_start)), std::abs(h_max));
 
     if (std::abs(H) <= 10 * parameters_.round_off_)
       H = delta_min_;
 
+    // TODO: the logic above this point should be moved to the constructor and should return an error
+    //       if the parameters are invalid (e.g., h_min > h_max or h_start > h_max)
+
     bool reject_last_h = false;
     bool reject_more_h = false;
 
@@ -639,7 +644,7 @@ namespace micm
           stats_.accepted += 1;
           present_time = present_time + H;
           Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end());
-          Hnew = std::max(parameters_.h_min_, std::min(Hnew, parameters_.h_max_));
+          Hnew = std::max(parameters_.h_min_, std::min(Hnew, h_max));
           if (reject_last_h)
           {
             // No step size increase after a rejected step

include/micm/system/phase.hpp

@@ -49,5 +49,21 @@ namespace micm
       species_ = other.species_;
       return *this;
     }
+
+    /// @brief Returns the number of non-parameterized species
+    std::size_t StateSize() const
+    {
+      return std::count_if(species_.begin(), species_.end(), [](const Species& s) { return !s.IsParameterized(); });
+    }
+
+    /// @brief Returns a set of unique names for each non-parameterized species
+    std::vector<std::string> UniqueNames() const
+    {
+      std::vector<std::string> names{};
+      for (const auto& species : species_)
+        if (!species.IsParameterized())
+          names.push_back(species.name_);
+      return names;
+    }
   };
 }  // namespace micm

include/micm/system/species.hpp

@@ -3,9 +3,12 @@
 // SPDX-License-Identifier: Apache-2.0
 #pragma once
 
+#include <functional>
 #include <string>
 #include <vector>
 
+#include "conditions.hpp"
+
 namespace micm
 {
 
@@ -19,6 +22,12 @@ namespace micm
     /// @brief A list of properties of this species
     std::map<std::string, double> properties_;
 
+    /// @brief A function that if provided will be used to parameterize
+    ///        the concentration of this species during solving.
+    ///        Species with this function defined will be excluded from
+    ///        the solver state.
+    std::function<double(const Conditions)> parameterize_{ nullptr };
+
     /// @brief Copy assignment
     /// @param other species to copy
     Species& operator=(const Species& other);
@@ -35,6 +44,12 @@ namespace micm
     /// @param name The name of the species
     /// @param properties The properties of the species
     Species(const std::string& name, const std::map<std::string, double>& properties);
+
+    /// @brief Returns whether a species is parameterized
+    bool IsParameterized() const;
+
+    /// @brief Return a Species instance parameterized on air density
+    static Species ThirdBody();
   };
 
   inline Species& Species::operator=(const Species& other)
@@ -46,13 +61,15 @@ namespace micm
 
     name_ = other.name_;
     properties_ = other.properties_;  // This performs a shallow copy
+    parameterize_ = other.parameterize_;
 
     return *this;
   }
 
   inline Species::Species(const Species& other)
       : name_(other.name_),
-        properties_(other.properties_){};
+        properties_(other.properties_),
+        parameterize_(other.parameterize_){};
 
   inline Species::Species(const std::string& name)
       : name_(name){};
@@ -61,4 +78,15 @@ namespace micm
       : name_(name),
         properties_(properties){};
 
+  inline bool Species::IsParameterized() const
+  {
+    return parameterize_ != nullptr;
+  }
+
+  inline Species Species::ThirdBody()
+  {
+    Species third_body{ "M" };
+    third_body.parameterize_ = [](const Conditions& c) { return c.air_density_; };
+    return third_body;
+  }
 }  // namespace micm